220 CHARACTER(len=255),
PARAMETER:: subcategory=
"grid_transform_class"
224 double precision :: boxdx
225 double precision :: boxdy
226 double precision :: radius
231 DOUBLE PRECISION :: percentile
240 END TYPE interval_info
272 TYPE(vol7d_level) :: input_levtype
273 TYPE(vol7d_var) :: input_coordvar
274 TYPE(vol7d_level) :: output_levtype
284 CHARACTER(len=80) :: trans_type
285 CHARACTER(len=80) :: sub_type
287 TYPE(rect_ind) :: rect_ind
288 TYPE(rect_coo) :: rect_coo
289 TYPE(area_info) :: area_info
290 TYPE(arrayof_georef_coord_array) :: poly
291 TYPE(stat_info) :: stat_info
292 TYPE(interval_info) :: interval_info
293 TYPE(box_info) :: box_info
294 TYPE(vertint) :: vertint
295 INTEGER :: time_definition
296 INTEGER :: category = 0
308 TYPE(transform_def),
PUBLIC :: trans
309 INTEGER :: innx = imiss
310 INTEGER :: inny = imiss
311 INTEGER :: innz = imiss
312 INTEGER :: outnx = imiss
313 INTEGER :: outny = imiss
314 INTEGER :: outnz = imiss
315 INTEGER :: levshift = imiss
316 INTEGER :: levused = imiss
317 INTEGER :: iniox,inioy,infox,infoy,outinx,outiny,outfnx,outfny
318 INTEGER,
POINTER :: inter_index_x(:,:) => null()
319 INTEGER,
POINTER :: inter_index_y(:,:) => null()
320 INTEGER,
POINTER :: inter_index_z(:) => null()
321 INTEGER,
POINTER :: point_index(:,:) => null()
322 DOUBLE PRECISION,
POINTER :: inter_x(:,:) => null()
323 DOUBLE PRECISION,
POINTER :: inter_y(:,:) => null()
324 DOUBLE PRECISION,
POINTER :: inter_xp(:,:) => null()
325 DOUBLE PRECISION,
POINTER :: inter_yp(:,:) => null()
326 DOUBLE PRECISION,
POINTER :: inter_zp(:) => null()
327 DOUBLE PRECISION,
POINTER :: vcoord_in(:) => null()
328 DOUBLE PRECISION,
POINTER :: vcoord_out(:) => null()
329 LOGICAL,
POINTER :: point_mask(:,:) => null()
330 LOGICAL,
POINTER :: stencil(:,:) => null()
331 REAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
332 REAL,
ALLOCATABLE :: val_mask(:,:)
335 LOGICAL :: recur = .false.
336 LOGICAL :: dolog = .false.
340 TYPE(vol7d_level),
POINTER :: output_level_auto(:) => null()
341 INTEGER :: category = 0
342 LOGICAL :: valid = .false.
343 PROCEDURE(basic_find_index),
NOPASS,
POINTER :: find_index => basic_find_index
349 MODULE PROCEDURE transform_init, grid_transform_levtype_levtype_init, &
350 grid_transform_init, &
351 grid_transform_grid_vol7d_init, grid_transform_vol7d_grid_init, &
352 grid_transform_vol7d_vol7d_init
357 MODULE PROCEDURE transform_delete, grid_transform_delete
362 MODULE PROCEDURE transform_get_val, grid_transform_get_val
368 MODULE PROCEDURE grid_transform_compute, grid_transform_v7d_grid_compute
374 MODULE PROCEDURE grid_transform_c_e
380 PUBLIC interval_info, interval_info_new, interval_info_valid, basic_find_index
385 FUNCTION interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
RESULT(this)
386 REAL,
INTENT(in),
OPTIONAL :: interv_gt
387 REAL,
INTENT(in),
OPTIONAL :: interv_ge
388 REAL,
INTENT(in),
OPTIONAL :: interv_lt
389 REAL,
INTENT(in),
OPTIONAL :: interv_le
391 TYPE(interval_info) :: this
393 IF (
PRESENT(interv_gt))
THEN
394 IF (
c_e(interv_gt))
THEN
399 IF (
PRESENT(interv_ge))
THEN
400 IF (
c_e(interv_ge))
THEN
405 IF (
PRESENT(interv_lt))
THEN
406 IF (
c_e(interv_lt))
THEN
411 IF (
PRESENT(interv_le))
THEN
412 IF (
c_e(interv_le))
THEN
418 END FUNCTION interval_info_new
425 ELEMENTAL FUNCTION interval_info_valid(this, val)
426 TYPE(interval_info),
INTENT(in) :: this
427 REAL,
INTENT(in) :: val
429 REAL :: interval_info_valid
431 interval_info_valid = 1.0
433 IF (
c_e(this%gt))
THEN
434 IF (val < this%gt) interval_info_valid = 0.0
435 IF (.NOT.this%ge)
THEN
436 IF (val == this%gt) interval_info_valid = 0.0
439 IF (
c_e(this%lt))
THEN
440 IF (val > this%lt) interval_info_valid = 0.0
441 IF (.NOT.this%le)
THEN
442 IF (val == this%lt) interval_info_valid = 0.0
446 END FUNCTION interval_info_valid
454 SUBROUTINE transform_init(this, trans_type, sub_type, &
455 ix, iy, fx, fy, ilon, ilat, flon, flat, &
456 npx, npy, boxdx, boxdy, radius, poly, percentile, &
457 interv_gt, interv_ge, interv_lt, interv_le, &
458 extrap, time_definition, &
459 input_levtype, input_coordvar, output_levtype, categoryappend)
460 TYPE(transform_def),
INTENT(out) :: this
461 CHARACTER(len=*) :: trans_type
462 CHARACTER(len=*) :: sub_type
463 INTEGER,
INTENT(in),
OPTIONAL :: ix
464 INTEGER,
INTENT(in),
OPTIONAL :: iy
465 INTEGER,
INTENT(in),
OPTIONAL :: fx
466 INTEGER,
INTENT(in),
OPTIONAL :: fy
467 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilon
468 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilat
469 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flon
470 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flat
471 INTEGER,
INTENT(IN),
OPTIONAL :: npx
472 INTEGER,
INTENT(IN),
OPTIONAL :: npy
473 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdx
474 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdy
475 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: radius
476 TYPE(arrayof_georef_coord_array),
OPTIONAL :: poly
477 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: percentile
478 REAL,
INTENT(in),
OPTIONAL :: interv_gt
479 REAL,
INTENT(in),
OPTIONAL :: interv_ge
480 REAL,
INTENT(in),
OPTIONAL :: interv_lt
481 REAL,
INTENT(in),
OPTIONAL :: interv_le
482 LOGICAL,
INTENT(IN),
OPTIONAL :: extrap
483 INTEGER,
INTENT(IN),
OPTIONAL :: time_definition
484 TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: input_levtype
485 TYPE(vol7d_var),
INTENT(IN),
OPTIONAL :: input_coordvar
486 TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: output_levtype
487 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
489 character(len=512) :: a_name
491 IF (
PRESENT(categoryappend))
THEN
492 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
493 trim(categoryappend))
495 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
497 this%category=l4f_category_get(a_name)
499 this%trans_type = trans_type
500 this%sub_type = sub_type
502 CALL optio(extrap,this%extrap)
504 call optio(ix,this%rect_ind%ix)
505 call optio(iy,this%rect_ind%iy)
506 call optio(fx,this%rect_ind%fx)
507 call optio(fy,this%rect_ind%fy)
509 call optio(ilon,this%rect_coo%ilon)
510 call optio(ilat,this%rect_coo%ilat)
511 call optio(flon,this%rect_coo%flon)
512 call optio(flat,this%rect_coo%flat)
514 CALL optio(boxdx,this%area_info%boxdx)
515 CALL optio(boxdy,this%area_info%boxdy)
516 CALL optio(radius,this%area_info%radius)
517 IF (
PRESENT(poly)) this%poly = poly
518 CALL optio(percentile,this%stat_info%percentile)
520 this%interval_info = interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
522 CALL optio(npx,this%box_info%npx)
523 CALL optio(npy,this%box_info%npy)
525 IF (
PRESENT(input_levtype))
THEN
526 this%vertint%input_levtype = input_levtype
528 this%vertint%input_levtype = vol7d_level_miss
530 IF (
PRESENT(input_coordvar))
THEN
531 this%vertint%input_coordvar = input_coordvar
533 this%vertint%input_coordvar = vol7d_var_miss
535 IF (
PRESENT(output_levtype))
THEN
536 this%vertint%output_levtype = output_levtype
538 this%vertint%output_levtype = vol7d_level_miss
541 call optio(time_definition,this%time_definition)
542 if (
c_e(this%time_definition) .and. &
543 (this%time_definition < 0 .OR. this%time_definition > 1))
THEN
544 call l4f_category_log(this%category,l4f_error,
"Error in time_definition: "//
to_char(this%time_definition))
545 call raise_fatal_error()
549 IF (this%trans_type ==
'zoom')
THEN
551 IF (this%sub_type ==
'coord' .OR. this%sub_type ==
'projcoord')
THEN
553 if (
c_e(this%rect_coo%ilon) .and.
c_e(this%rect_coo%ilat) .and. &
554 c_e(this%rect_coo%flon) .and.
c_e(this%rect_coo%flat))
then
557 if ( this%rect_coo%ilon > this%rect_coo%flon .or. &
558 this%rect_coo%ilat > this%rect_coo%flat )
then
560 call l4f_category_log(this%category,l4f_error, &
561 "invalid zoom coordinates: ")
562 call l4f_category_log(this%category,l4f_error, &
563 trim(
to_char(this%rect_coo%ilon))//
'/'// &
564 trim(
to_char(this%rect_coo%flon)))
565 call l4f_category_log(this%category,l4f_error, &
566 trim(
to_char(this%rect_coo%ilat))//
'/'// &
568 call raise_fatal_error()
573 call l4f_category_log(this%category,l4f_error,
"zoom: coord parameters missing")
574 call raise_fatal_error()
578 else if (this%sub_type ==
'coordbb')
then
580 if (
c_e(this%rect_coo%ilon) .and.
c_e(this%rect_coo%ilat) .and. &
581 c_e(this%rect_coo%flon) .and.
c_e(this%rect_coo%flat))
then
584 call l4f_category_log(this%category,l4f_error,
"zoom: coordbb parameters missing")
585 call raise_fatal_error()
589 else if (this%sub_type ==
'index')
then
591 IF (
c_e(this%rect_ind%ix) .AND.
c_e(this%rect_ind%iy) .AND. &
592 c_e(this%rect_ind%fx) .AND.
c_e(this%rect_ind%fy))
THEN
595 IF (this%rect_ind%ix > this%rect_ind%fx .OR. &
596 this%rect_ind%iy > this%rect_ind%fy)
THEN
598 CALL l4f_category_log(this%category,l4f_error,
'invalid zoom indices: ')
599 CALL l4f_category_log(this%category,l4f_error, &
600 trim(
to_char(this%rect_ind%ix))//
'/'// &
601 trim(
to_char(this%rect_ind%fx)))
602 CALL l4f_category_log(this%category,l4f_error, &
603 trim(
to_char(this%rect_ind%iy))//
'/'// &
604 trim(
to_char(this%rect_ind%fy)))
606 CALL raise_fatal_error()
611 CALL l4f_category_log(this%category,l4f_error,&
612 'zoom: index parameters ix, iy, fx, fy not provided')
613 CALL raise_fatal_error()
618 CALL sub_type_error()
622 ELSE IF (this%trans_type ==
'inter')
THEN
624 IF (this%sub_type ==
'near' .OR. this%sub_type ==
'bilin' .OR. &
625 this%sub_type ==
'linear' .OR. this%sub_type ==
'shapiro_near')
THEN
628 CALL sub_type_error()
632 ELSE IF (this%trans_type ==
'boxregrid' .OR. this%trans_type ==
'boxinter' .OR. &
633 this%trans_type ==
'polyinter' .OR. this%trans_type ==
'maskinter' .OR. &
634 this%trans_type ==
'stencilinter')
THEN
636 IF (this%trans_type ==
'boxregrid')
THEN
637 IF (
c_e(this%box_info%npx) .AND.
c_e(this%box_info%npy))
THEN
638 IF (this%box_info%npx <= 0 .OR. this%box_info%npy <= 0 )
THEN
639 CALL l4f_category_log(this%category,l4f_error,
'boxregrid: invalid parameters '//&
640 trim(
to_char(this%box_info%npx))//
' '//trim(
to_char(this%box_info%npy)))
641 CALL raise_fatal_error()
644 CALL l4f_category_log(this%category,l4f_error, &
645 'boxregrid: parameters npx, npy missing')
646 CALL raise_fatal_error()
650 IF (this%trans_type ==
'polyinter')
THEN
651 IF (this%poly%arraysize <= 0)
THEN
652 CALL l4f_category_log(this%category,l4f_error, &
653 "polyinter: poly parameter missing or empty")
654 CALL raise_fatal_error()
658 IF (this%trans_type ==
'stencilinter')
THEN
659 IF (.NOT.
c_e(this%area_info%radius))
THEN
660 CALL l4f_category_log(this%category,l4f_error, &
661 "stencilinter: radius parameter missing")
662 CALL raise_fatal_error()
666 IF (this%sub_type ==
'average' .OR. this%sub_type ==
'stddev' &
667 .OR. this%sub_type ==
'stddevnm1')
THEN
668 this%stat_info%percentile = rmiss
669 ELSE IF (this%sub_type ==
'max')
THEN
670 this%stat_info%percentile = 101.
671 ELSE IF (this%sub_type ==
'min')
THEN
672 this%stat_info%percentile = -1.
673 ELSE IF (this%sub_type ==
'percentile')
THEN
674 IF (.NOT.
c_e(this%stat_info%percentile))
THEN
675 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
676 ':percentile: percentile value not provided')
677 CALL raise_fatal_error()
678 ELSE IF (this%stat_info%percentile >= 100.)
THEN
679 this%sub_type =
'max'
680 ELSE IF (this%stat_info%percentile <= 0.)
THEN
681 this%sub_type =
'min'
683 ELSE IF (this%sub_type ==
'frequency')
THEN
684 IF (.NOT.
c_e(this%interval_info%gt) .AND. .NOT.
c_e(this%interval_info%gt))
THEN
685 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
686 ':frequency: lower and/or upper limit not provided')
687 CALL raise_fatal_error()
690 CALL sub_type_error()
694 ELSE IF (this%trans_type ==
'maskgen')
THEN
696 IF (this%sub_type ==
'poly')
THEN
698 IF (this%poly%arraysize <= 0)
THEN
699 CALL l4f_category_log(this%category,l4f_error,
"maskgen:poly poly parameter missing or empty")
700 CALL raise_fatal_error()
703 ELSE IF (this%sub_type ==
'grid')
THEN
707 CALL sub_type_error()
711 ELSE IF (this%trans_type ==
'vertint')
THEN
713 IF (this%vertint%input_levtype == vol7d_level_miss)
THEN
714 CALL l4f_category_log(this%category,l4f_error, &
715 'vertint parameter input_levtype not provided')
716 CALL raise_fatal_error()
719 IF (this%vertint%output_levtype == vol7d_level_miss)
THEN
720 CALL l4f_category_log(this%category,l4f_error, &
721 'vertint parameter output_levtype not provided')
722 CALL raise_fatal_error()
725 IF (this%sub_type ==
'linear' .OR. this%sub_type ==
'linearsparse')
THEN
728 CALL sub_type_error()
732 ELSE IF (this%trans_type ==
'metamorphosis')
THEN
734 IF (this%sub_type ==
'all')
THEN
736 ELSE IF (this%sub_type ==
'coordbb')
THEN
738 IF (
c_e(this%rect_coo%ilon) .AND.
c_e(this%rect_coo%ilat) .AND. &
739 c_e(this%rect_coo%flon) .AND.
c_e(this%rect_coo%flat))
THEN
742 CALL l4f_category_log(this%category,l4f_error,
"metamorphosis: coordbb parameters missing")
743 CALL raise_fatal_error()
747 ELSE IF (this%sub_type ==
'poly')
THEN
749 IF (this%poly%arraysize <= 0)
THEN
750 CALL l4f_category_log(this%category,l4f_error,
"metamorphosis:poly: poly parameter missing or empty")
751 CALL raise_fatal_error()
754 ELSE IF (this%sub_type ==
'mask' .OR. this%sub_type ==
'maskvalid' .OR. &
755 this%sub_type ==
'maskinvalid' .OR. this%sub_type ==
'setinvalidto' .OR. &
756 this%sub_type ==
'settoinvalid' .OR. this%sub_type ==
'lemaskinvalid' .OR. &
757 this%sub_type ==
'ltmaskinvalid' .OR. this%sub_type ==
'gemaskinvalid' .OR. &
758 this%sub_type ==
'gtmaskinvalid')
THEN
761 CALL sub_type_error()
766 CALL trans_type_error()
772 SUBROUTINE sub_type_error()
774 CALL l4f_category_log(this%category, l4f_error, trim(this%trans_type) &
775 //
': sub_type '//trim(this%sub_type)//
' is not defined')
776 CALL raise_fatal_error()
778 END SUBROUTINE sub_type_error
780 SUBROUTINE trans_type_error()
782 CALL l4f_category_log(this%category, l4f_error,
'trans_type '//this%trans_type &
784 CALL raise_fatal_error()
786 END SUBROUTINE trans_type_error
789 END SUBROUTINE transform_init
795 SUBROUTINE transform_delete(this)
798 this%trans_type=cmiss
801 this%rect_ind%ix=imiss
802 this%rect_ind%iy=imiss
803 this%rect_ind%fx=imiss
804 this%rect_ind%fy=imiss
806 this%rect_coo%ilon=dmiss
807 this%rect_coo%ilat=dmiss
808 this%rect_coo%flon=dmiss
809 this%rect_coo%flat=dmiss
811 this%box_info%npx=imiss
812 this%box_info%npy=imiss
817 CALL l4f_category_delete(this%category)
819 END SUBROUTINE transform_delete
823 SUBROUTINE transform_get_val(this, time_definition, trans_type, sub_type, &
824 input_levtype, output_levtype)
826 INTEGER,
INTENT(out),
OPTIONAL :: time_definition
827 CHARACTER(len=*),
INTENT(out),
OPTIONAL :: trans_type
828 CHARACTER(len=*),
INTENT(out),
OPTIONAL :: sub_type
829 TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: input_levtype
831 TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: output_levtype
834 IF (
PRESENT(time_definition)) time_definition=this%time_definition
835 IF (
PRESENT(trans_type)) trans_type = this%trans_type
836 IF (
PRESENT(sub_type)) sub_type = this%sub_type
837 IF (
PRESENT(input_levtype)) input_levtype = this%vertint%input_levtype
838 IF (
PRESENT(output_levtype)) output_levtype = this%vertint%output_levtype
841 END SUBROUTINE transform_get_val
887 SUBROUTINE grid_transform_levtype_levtype_init(this, trans, lev_in, lev_out, &
888 coord_3d_in, categoryappend)
891 TYPE(vol7d_level),
INTENT(in) :: lev_in(:)
892 TYPE(vol7d_level),
INTENT(in) :: lev_out(:)
893 REAL,
INTENT(inout),
OPTIONAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
894 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
896 DOUBLE PRECISION :: coord_in(SIZE(lev_in))
897 DOUBLE PRECISION,
ALLOCATABLE :: coord_out(:)
898 LOGICAL :: mask_in(SIZE(lev_in))
899 LOGICAL,
ALLOCATABLE :: mask_out(:)
901 INTEGER :: i, j, icache, inused, istart, iend, ostart, oend
904 CALL grid_transform_init_common(this, trans, categoryappend)
906 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vertint")
909 IF (this%trans%trans_type ==
'vertint')
THEN
911 IF (
c_e(trans%vertint%input_levtype%level2) .AND. &
912 trans%vertint%input_levtype%level1 /= trans%vertint%input_levtype%level2)
THEN
913 CALL l4f_category_log(this%category, l4f_error, &
914 'vertint: input upper and lower surface must be of the same type, '// &
915 t2c(trans%vertint%input_levtype%level1)//
'/='// &
916 t2c(trans%vertint%input_levtype%level2))
920 IF (
c_e(trans%vertint%output_levtype%level2) .AND. &
921 trans%vertint%output_levtype%level1 /= trans%vertint%output_levtype%level2)
THEN
922 CALL l4f_category_log(this%category, l4f_error, &
923 'vertint: output upper and lower surface must be of the same type'// &
924 t2c(trans%vertint%output_levtype%level1)//
'/='// &
925 t2c(trans%vertint%output_levtype%level2))
930 mask_in(:) = (lev_in(:)%level1 == trans%vertint%input_levtype%level1) .AND. &
931 (lev_in(:)%level2 == trans%vertint%input_levtype%level2)
932 CALL make_vert_coord(lev_in, mask_in, coord_in, dolog)
933 this%innz =
SIZE(lev_in)
934 istart = firsttrue(mask_in)
935 iend = lasttrue(mask_in)
936 inused = iend - istart + 1
937 IF (inused /= count(mask_in))
THEN
938 CALL l4f_category_log(this%category, l4f_error, &
939 'grid_transform_levtype_levtype_init: input levels badly sorted '//&
940 t2c(inused)//
'/'//t2c(count(mask_in)))
944 this%levshift = istart-1
945 this%levused = inused
947 IF (trans%vertint%input_levtype%level1 /= trans%vertint%output_levtype%level1)
THEN
949 CALL l4f_category_log(this%category, l4f_debug, &
950 'vertint: different input and output level types '// &
951 t2c(trans%vertint%input_levtype%level1)//
' '// &
952 t2c(trans%vertint%output_levtype%level1))
955 ALLOCATE(mask_out(
SIZE(lev_out)), this%vcoord_out(
SIZE(lev_out)))
956 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
957 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
958 CALL make_vert_coord(lev_out, mask_out, this%vcoord_out, dolog)
959 this%outnz =
SIZE(mask_out)
962 IF (.NOT.
PRESENT(coord_3d_in))
THEN
963 CALL l4f_category_log(this%category, l4f_warn, &
964 'vertint: different input and output level types &
965 &and no coord_3d_in, expecting vert. coord. in volume')
968 IF (
SIZE(coord_3d_in,3) /= inused)
THEN
969 CALL l4f_category_log(this%category, l4f_error, &
970 'vertint: vertical size of coord_3d_in (vertical coordinate) &
971 &different from number of input levels suitable for interpolation')
972 CALL l4f_category_log(this%category, l4f_error, &
973 'coord_3d_in: '//t2c(
SIZE(coord_3d_in,3))// &
974 ', input levels for interpolation: '//t2c(inused))
979 CALL move_alloc(coord_3d_in, this%coord_3d_in)
981 WHERE(
c_e(this%coord_3d_in) .AND. this%coord_3d_in > 0.0)
982 this%coord_3d_in = log(this%coord_3d_in)
984 this%coord_3d_in = rmiss
995 CALL l4f_category_log(this%category, l4f_debug, &
996 'vertint: equal input and output level types '// &
997 t2c(trans%vertint%input_levtype%level1))
1000 IF (
SIZE(lev_out) > 0)
THEN
1001 ALLOCATE(mask_out(
SIZE(lev_out)), coord_out(
SIZE(lev_out)))
1002 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
1003 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
1004 CALL make_vert_coord(lev_out, mask_out, coord_out, dolog)
1007 IF (
c_e(trans%vertint%input_levtype%level2) .AND. &
1008 .NOT.
c_e(trans%vertint%output_levtype%level2))
THEN
1009 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
1010 trans%vertint%output_levtype%level1 == 150)
THEN
1011 ALLOCATE(this%output_level_auto(inused-1))
1012 CALL l4f_category_log(this%category,l4f_info, &
1013 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
1014 //
'/'//t2c(iend-istart)//
' output levels (f->h)')
1015 DO i = istart, iend - 1
1016 CALL init(this%output_level_auto(i-istart+1), &
1017 trans%vertint%input_levtype%level1, lev_in(i)%l2)
1020 CALL l4f_category_log(this%category, l4f_error, &
1021 'grid_transform_levtype_levtype_init: automatic generation of output levels &
1022 &available only for hybrid levels')
1026 ELSE IF (.NOT.
c_e(trans%vertint%input_levtype%level2) .AND. &
1027 c_e(trans%vertint%output_levtype%level2))
THEN
1028 ALLOCATE(this%output_level_auto(inused-1))
1029 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
1030 trans%vertint%output_levtype%level1 == 150)
THEN
1031 CALL l4f_category_log(this%category,l4f_info, &
1032 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
1033 //
'/'//t2c(iend-istart)//
' output levels (h->f)')
1034 DO i = istart, iend - 1
1035 CALL init(this%output_level_auto(i-istart+1), trans%vertint%input_levtype%level1, &
1036 lev_in(i)%l1, trans%vertint%input_levtype%level1, &
1040 CALL l4f_category_log(this%category, l4f_error, &
1041 'grid_transform_levtype_levtype_init: automatic generation of output levels &
1042 &available only for hybrid levels')
1047 CALL l4f_category_log(this%category, l4f_error, &
1048 'grid_transform_levtype_levtype_init: strange situation'// &
1049 to_char(
c_e(trans%vertint%input_levtype%level2))//
' '// &
1050 to_char(
c_e(trans%vertint%output_levtype%level2)))
1054 ALLOCATE(coord_out(inused-1), mask_out(inused-1))
1055 mask_out(:) = .true.
1056 CALL make_vert_coord(this%output_level_auto, mask_out, coord_out, dolog)
1059 this%outnz =
SIZE(mask_out)
1060 ostart = firsttrue(mask_out)
1061 oend = lasttrue(mask_out)
1064 IF (istart == 0)
THEN
1065 CALL l4f_category_log(this%category, l4f_warn, &
1066 'grid_transform_levtype_levtype_init: &
1067 &input contains no vertical levels of type ('// &
1068 trim(
to_char(trans%vertint%input_levtype%level1))//
','// &
1069 trim(
to_char(trans%vertint%input_levtype%level2))// &
1070 ') suitable for interpolation')
1073 ELSE IF (istart == iend)
THEN
1074 CALL l4f_category_log(this%category, l4f_warn, &
1075 'grid_transform_levtype_levtype_init: &
1076 &input contains only 1 vertical level of type ('// &
1077 trim(
to_char(trans%vertint%input_levtype%level1))//
','// &
1078 trim(
to_char(trans%vertint%input_levtype%level2))// &
1079 ') suitable for interpolation')
1081 IF (ostart == 0)
THEN
1082 CALL l4f_category_log(this%category, l4f_warn, &
1083 'grid_transform_levtype_levtype_init: &
1084 &output contains no vertical levels of type ('// &
1085 trim(
to_char(trans%vertint%output_levtype%level1))//
','// &
1086 trim(
to_char(trans%vertint%output_levtype%level2))// &
1087 ') suitable for interpolation')
1093 IF (this%trans%sub_type ==
'linear')
THEN
1095 ALLOCATE(this%inter_index_z(this%outnz), this%inter_zp(this%outnz))
1096 this%inter_index_z(:) = imiss
1097 this%inter_zp(:) = dmiss
1098 IF (this%trans%extrap .AND. istart > 0)
THEN
1101 this%inter_index_z(:) = istart
1102 this%inter_zp(:) = 1.0d0
1107 outlev:
DO j = ostart, oend
1108 inlev:
DO i = icache, iend
1109 IF (coord_in(i) >= coord_out(j))
THEN
1110 IF (coord_out(j) >= coord_in(i-1))
THEN
1111 this%inter_index_z(j) = i - 1
1112 this%inter_zp(j) = (coord_out(j)-coord_in(i-1)) / &
1113 (coord_in(i)-coord_in(i-1))
1120 IF (this%trans%extrap .AND. iend > 1)
THEN
1121 this%inter_index_z(j) = iend - 1
1122 this%inter_zp(j) = 0.0d0
1126 DEALLOCATE(coord_out, mask_out)
1129 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN
1131 ALLOCATE(this%vcoord_in(this%levused), this%vcoord_out(
SIZE(coord_out)))
1132 this%vcoord_in(:) = coord_in(this%levshift+1:this%levshift+this%levused)
1133 this%vcoord_out(:) = coord_out(:)
1134 DEALLOCATE(coord_out, mask_out)
1145 END SUBROUTINE grid_transform_levtype_levtype_init
1150 SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1151 TYPE(vol7d_level),
INTENT(in) :: lev(:)
1152 LOGICAL,
INTENT(inout) :: mask(:)
1153 DOUBLE PRECISION,
INTENT(out) :: coord(:)
1154 LOGICAL,
INTENT(out) :: dolog
1157 DOUBLE PRECISION :: fact
1164 IF (any(lev(k)%level1 == height_level))
THEN
1166 ELSE IF (any(lev(k)%level1 == thermo_level))
THEN
1168 ELSE IF (any(lev(k)%level1 == sigma_level))
THEN
1174 IF (
c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2)
THEN
1175 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108)
THEN
1176 WHERE(mask(:) .AND. lev(:)%l1 > 0 .AND. lev(:)%l2 > 0)
1177 coord(:) = (log(dble(lev(:)%l1)*fact) + log(dble(lev(:)%l2)*fact))*0.5d0
1182 coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0
1186 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108)
THEN
1187 WHERE(mask(:) .AND. lev(:)%l1 > 0)
1188 coord(:) = log(dble(lev(:)%l1)*fact)
1193 coord(:) = lev(:)%l1*fact
1199 mask(:) = mask(:) .AND.
c_e(coord(:))
1201 END SUBROUTINE make_vert_coord
1221 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1227 REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1228 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1229 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1231 INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1232 xnmin, xnmax, ynmin, ynmax
1233 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1234 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1235 TYPE(geo_proj) :: proj_in, proj_out
1237 LOGICAL,
ALLOCATABLE :: point_mask(:,:)
1241 CALL grid_transform_init_common(this, trans, categoryappend)
1243 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-vg6d")
1247 CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1248 CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1250 IF (this%trans%trans_type ==
'zoom')
THEN
1252 IF (this%trans%sub_type ==
'coord')
THEN
1254 CALL griddim_zoom_coord(in, &
1255 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1256 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1257 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1258 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1260 ELSE IF (this%trans%sub_type ==
'projcoord')
THEN
1262 CALL griddim_zoom_projcoord(in, &
1263 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1264 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1265 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1266 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1268 ELSE IF (this%trans%sub_type ==
'coordbb')
THEN
1273 CALL get_val(lin, nx=nx, ny=ny)
1275 ALLOCATE(point_mask(nx,ny))
1276 point_mask(:,:) = .false.
1282 IF (lin%dim%lon(i,j) > this%trans%rect_coo%ilon .AND. &
1283 lin%dim%lon(i,j) < this%trans%rect_coo%flon .AND. &
1284 lin%dim%lat(i,j) > this%trans%rect_coo%ilat .AND. &
1285 lin%dim%lat(i,j) < this%trans%rect_coo%flat)
THEN
1286 point_mask(i,j) = .true.
1293 IF (any(point_mask(i,:)))
EXIT
1295 this%trans%rect_ind%ix = i
1296 DO i = nx, this%trans%rect_ind%ix, -1
1297 IF (any(point_mask(i,:)))
EXIT
1299 this%trans%rect_ind%fx = i
1302 IF (any(point_mask(:,j)))
EXIT
1304 this%trans%rect_ind%iy = j
1305 DO j = ny, this%trans%rect_ind%iy, -1
1306 IF (any(point_mask(:,j)))
EXIT
1308 this%trans%rect_ind%fy = j
1310 DEALLOCATE(point_mask)
1312 IF (this%trans%rect_ind%ix > this%trans%rect_ind%fx .OR. &
1313 this%trans%rect_ind%iy > this%trans%rect_ind%fy)
THEN
1315 CALL l4f_category_log(this%category,l4f_error, &
1316 "zoom coordbb: no points inside bounding box "//&
1317 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
1318 trim(
to_char(this%trans%rect_coo%flon))//
","// &
1319 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
1320 trim(
to_char(this%trans%rect_coo%flat)))
1329 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1330 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1333 this%iniox = min(max(this%trans%rect_ind%ix,1),nx)
1334 this%inioy = min(max(this%trans%rect_ind%iy,1),ny)
1335 this%infox = max(min(this%trans%rect_ind%fx,nx),1)
1336 this%infoy = max(min(this%trans%rect_ind%fy,ny),1)
1338 this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx)
1339 this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny)
1340 this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1
1341 this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1
1343 xmin=xmin+steplon*(this%trans%rect_ind%ix-1)
1344 ymin=ymin+steplat*(this%trans%rect_ind%iy-1)
1345 xmax=xmax+steplon*(this%trans%rect_ind%fx-nx)
1346 ymax=ymax+steplat*(this%trans%rect_ind%fy-ny)
1350 CALL dealloc(out%dim)
1352 out%dim%nx = this%trans%rect_ind%fx - this%trans%rect_ind%ix + 1
1353 out%dim%ny = this%trans%rect_ind%fy - this%trans%rect_ind%iy + 1
1354 this%outnx = out%dim%nx
1355 this%outny = out%dim%ny
1359 CALL set_val(out, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1363 ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN
1365 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1366 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1372 xmin_new = xmin + (this%trans%box_info%npx - 1)*0.5d0*steplon
1373 ymin_new = ymin + (this%trans%box_info%npy - 1)*0.5d0*steplat
1377 CALL dealloc(out%dim)
1379 out%dim%nx = nx/this%trans%box_info%npx
1380 out%dim%ny = ny/this%trans%box_info%npy
1381 this%outnx = out%dim%nx
1382 this%outny = out%dim%ny
1383 steplon = steplon*this%trans%box_info%npx
1384 steplat = steplat*this%trans%box_info%npy
1386 CALL set_val(out, xmin=xmin_new, ymin=ymin_new, &
1387 xmax=xmin_new + dble(out%dim%nx-1)*steplon, dx=steplon, &
1388 ymax=ymin_new + dble(out%dim%ny-1)*steplat, dy=steplat)
1392 ELSE IF (this%trans%trans_type ==
'inter')
THEN
1394 CALL outgrid_setup()
1396 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin'&
1397 .OR. this%trans%sub_type ==
'shapiro_near')
THEN
1399 CALL get_val(in, nx=this%innx, ny=this%inny, &
1400 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1401 CALL get_val(out, nx=this%outnx, ny=this%outny)
1403 ALLOCATE(this%inter_index_x(this%outnx,this%outny), &
1404 this%inter_index_y(this%outnx,this%outny))
1405 CALL copy(out, lout)
1408 IF (this%trans%sub_type ==
'bilin')
THEN
1409 CALL this%find_index(in, .false., &
1410 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1411 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1412 this%inter_index_x, this%inter_index_y)
1414 ALLOCATE(this%inter_x(this%innx,this%inny), &
1415 this%inter_y(this%innx,this%inny))
1416 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1417 this%inter_yp(this%outnx,this%outny))
1420 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1422 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1425 CALL this%find_index(in, .true., &
1426 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1427 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1428 this%inter_index_x, this%inter_index_y)
1436 ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
1438 CALL outgrid_setup()
1439 CALL get_val(in, nx=this%innx, ny=this%inny)
1440 CALL get_val(out, nx=this%outnx, ny=this%outny, &
1441 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1444 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
1445 CALL get_val(out, dx=this%trans%area_info%boxdx)
1446 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
1447 CALL get_val(out, dx=this%trans%area_info%boxdy)
1449 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1450 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1452 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1453 this%inter_index_y(this%innx,this%inny))
1459 CALL this%find_index(out, .true., &
1460 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1461 lin%dim%lon, lin%dim%lat, .false., &
1462 this%inter_index_x, this%inter_index_y)
1467 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
1469 CALL outgrid_setup()
1471 CALL get_val(in, nx=this%innx, ny=this%inny, &
1472 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1473 CALL get_val(out, nx=this%outnx, ny=this%outny)
1475 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1476 this%inter_index_y(this%outnx,this%outny))
1477 CALL copy(out, lout)
1479 CALL this%find_index(in, .true., &
1480 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1481 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1482 this%inter_index_x, this%inter_index_y)
1485 nr = int(this%trans%area_info%radius)
1488 r2 = this%trans%area_info%radius**2
1489 ALLOCATE(this%stencil(n,n))
1490 this%stencil(:,:) = .true.
1493 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1500 xnmax = this%innx - nr
1502 ynmax = this%inny - nr
1503 DO iy = 1, this%outny
1504 DO ix = 1, this%outnx
1505 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1506 this%inter_index_x(ix,iy) > xnmax .OR. &
1507 this%inter_index_y(ix,iy) < ynmin .OR. &
1508 this%inter_index_y(ix,iy) > ynmax)
THEN
1509 this%inter_index_x(ix,iy) = imiss
1510 this%inter_index_y(ix,iy) = imiss
1516 CALL l4f_category_log(this%category, l4f_debug, &
1517 'stencilinter: stencil size '//t2c(n*n))
1518 CALL l4f_category_log(this%category, l4f_debug, &
1519 'stencilinter: stencil points '//t2c(count(this%stencil)))
1525 ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
1527 IF (this%trans%sub_type ==
'poly')
THEN
1530 CALL get_val(in, nx=this%innx, ny=this%inny)
1531 this%outnx = this%innx
1532 this%outny = this%inny
1535 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1536 this%inter_index_y(this%innx,this%inny))
1537 this%inter_index_x(:,:) = imiss
1538 this%inter_index_y(:,:) = 1
1547 inside_x:
DO i = 1, this%innx
1548 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1550 DO n = nprev, this%trans%poly%arraysize
1551 IF (
inside(point, this%trans%poly%array(n)))
THEN
1552 this%inter_index_x(i,j) = n
1557 DO n = nprev-1, 1, -1
1558 IF (
inside(point, this%trans%poly%array(n)))
THEN
1559 this%inter_index_x(i,j) = n
1570 ELSE IF (this%trans%sub_type ==
'grid')
THEN
1573 CALL copy(out, lout)
1576 CALL get_val(in, nx=this%innx, ny=this%inny)
1577 this%outnx = this%innx
1578 this%outny = this%inny
1579 CALL get_val(lout, nx=nx, ny=ny, &
1580 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1583 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1584 this%inter_index_y(this%innx,this%inny))
1590 CALL this%find_index(lout, .true., &
1591 nx, ny, xmin, xmax, ymin, ymax, &
1592 out%dim%lon, out%dim%lat, .false., &
1593 this%inter_index_x, this%inter_index_y)
1595 WHERE(
c_e(this%inter_index_x(:,:)))
1596 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1597 (this%inter_index_y(:,:)-1)*nx
1605 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
1611 CALL get_val(in, nx=this%innx, ny=this%inny)
1612 this%outnx = this%innx
1613 this%outny = this%inny
1616 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1617 this%inter_index_y(this%innx,this%inny))
1618 this%inter_index_x(:,:) = imiss
1619 this%inter_index_y(:,:) = 1
1628 inside_x_2:
DO i = 1, this%innx
1629 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1631 DO n = nprev, this%trans%poly%arraysize
1632 IF (
inside(point, this%trans%poly%array(n)))
THEN
1633 this%inter_index_x(i,j) = n
1638 DO n = nprev-1, 1, -1
1639 IF (
inside(point, this%trans%poly%array(n)))
THEN
1640 this%inter_index_x(i,j) = n
1653 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
1656 CALL get_val(in, nx=this%innx, ny=this%inny)
1657 this%outnx = this%innx
1658 this%outny = this%inny
1660 IF (this%trans%sub_type ==
'maskvalid' .OR. this%trans%sub_type ==
'maskinvalid')
THEN
1662 IF (.NOT.
PRESENT(maskgrid))
THEN
1663 CALL l4f_category_log(this%category,l4f_error, &
1664 'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1665 trim(this%trans%sub_type)//
' transformation')
1670 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
1671 CALL l4f_category_log(this%category,l4f_error, &
1672 'grid_transform_init mask non conformal with input field')
1673 CALL l4f_category_log(this%category,l4f_error, &
1674 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
1675 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
1680 ALLOCATE(this%point_mask(this%innx,this%inny))
1682 IF (this%trans%sub_type ==
'maskvalid')
THEN
1685 IF (.NOT.
PRESENT(maskbounds))
THEN
1686 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1687 ELSE IF (
SIZE(maskbounds) < 2)
THEN
1688 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1690 this%point_mask(:,:) =
c_e(maskgrid(:,:)) .AND. &
1691 maskgrid(:,:) > maskbounds(1) .AND. &
1692 maskgrid(:,:) <= maskbounds(
SIZE(maskbounds))
1695 this%point_mask(:,:) = .NOT.
c_e(maskgrid(:,:))
1700 ELSE IF (this%trans%sub_type ==
'lemaskinvalid' .OR. &
1701 this%trans%sub_type ==
'ltmaskinvalid' .OR. &
1702 this%trans%sub_type ==
'gemaskinvalid' .OR. &
1703 this%trans%sub_type ==
'gtmaskinvalid')
THEN
1706 this%val_mask = maskgrid
1709 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
1711 IF (.NOT.
PRESENT(maskbounds))
THEN
1712 CALL l4f_category_log(this%category,l4f_error, &
1713 'grid_transform_init maskbounds missing for metamorphosis:'// &
1714 trim(this%trans%sub_type)//
' transformation')
1716 ELSE IF (
SIZE(maskbounds) < 1)
THEN
1717 CALL l4f_category_log(this%category,l4f_error, &
1718 'grid_transform_init maskbounds empty for metamorphosis:'// &
1719 trim(this%trans%sub_type)//
' transformation')
1722 this%val1 = maskbounds(1)
1724 CALL l4f_category_log(this%category, l4f_debug, &
1725 "grid_transform_init setting invalid data to "//t2c(this%val1))
1731 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
1733 IF (.NOT.
PRESENT(maskbounds))
THEN
1734 CALL l4f_category_log(this%category,l4f_error, &
1735 'grid_transform_init maskbounds missing for metamorphosis:'// &
1736 trim(this%trans%sub_type)//
' transformation')
1739 ELSE IF (
SIZE(maskbounds) < 2)
THEN
1740 CALL l4f_category_log(this%category,l4f_error, &
1741 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1742 trim(this%trans%sub_type)//
' transformation')
1746 this%val1 = maskbounds(1)
1747 this%val2 = maskbounds(
SIZE(maskbounds))
1749 CALL l4f_category_log(this%category, l4f_debug, &
1750 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
1751 t2c(this%val2)//
']')
1765 SUBROUTINE outgrid_setup()
1768 CALL griddim_setsteps(out)
1770 CALL get_val(in,
proj=proj_in, component_flag=cf_i)
1771 CALL get_val(out,
proj=proj_out, component_flag=cf_o)
1772 IF (proj_in == proj_out)
THEN
1775 CALL set_val(out, component_flag=cf_i)
1780 CALL l4f_category_log(this%category,l4f_warn, &
1781 'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1782 CALL l4f_category_log(this%category,l4f_warn, &
1783 'vector fields will probably be wrong')
1785 CALL set_val(out, component_flag=cf_i)
1789 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1791 END SUBROUTINE outgrid_setup
1793 END SUBROUTINE grid_transform_init
1838 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1839 maskgrid, maskbounds, find_index, categoryappend)
1843 TYPE(
vol7d),
INTENT(inout) :: v7d_out
1844 REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1845 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1846 PROCEDURE(basic_find_index),
POINTER,
OPTIONAL :: find_index
1847 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1849 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1851 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1852 DOUBLE PRECISION,
ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1853 REAL,
ALLOCATABLE :: lmaskbounds(:)
1858 IF (
PRESENT(find_index))
THEN
1859 IF (
ASSOCIATED(find_index))
THEN
1860 this%find_index => find_index
1863 CALL grid_transform_init_common(this, trans, categoryappend)
1865 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-v7d")
1869 CALL get_val(trans, time_definition=time_definition)
1870 IF (.NOT.
c_e(time_definition))
THEN
1874 IF (this%trans%trans_type ==
'inter')
THEN
1876 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin' &
1877 .OR. this%trans%sub_type ==
'shapiro_near')
THEN
1880 CALL get_val(lin, nx=this%innx, ny=this%inny)
1881 this%outnx =
SIZE(v7d_out%ana)
1884 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1885 this%inter_index_y(this%outnx,this%outny))
1886 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1888 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1890 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1891 CALL griddim_set_central_lon(lin, lonref)
1892 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1894 IF (this%trans%sub_type ==
'bilin')
THEN
1895 CALL this%find_index(lin, .false., &
1896 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1897 lon, lat, this%trans%extrap, &
1898 this%inter_index_x, this%inter_index_y)
1900 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1901 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1903 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1904 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1907 CALL this%find_index(lin, .true., &
1908 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1909 lon, lat, this%trans%extrap, &
1910 this%inter_index_x, this%inter_index_y)
1921 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
1923 CALL get_val(in, nx=this%innx, ny=this%inny)
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
1934 this%outnx = this%trans%poly%arraysize
1937 CALL init(v7d_out, time_definition=time_definition)
1938 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1941 ALLOCATE(lon(this%outnx,1))
1942 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
1943 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1944 CALL griddim_set_central_lon(lin, lonref)
1949 CALL poly_to_coordinates(this%trans%poly, v7d_out)
1954 DO iy = 1, this%inny
1955 inside_x:
DO ix = 1, this%innx
1956 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
1958 DO n = nprev, this%trans%poly%arraysize
1959 IF (
inside(point, this%trans%poly%array(n)))
THEN
1960 this%inter_index_x(ix,iy) = n
1965 DO n = nprev-1, 1, -1
1966 IF (
inside(point, this%trans%poly%array(n)))
THEN
1967 this%inter_index_x(ix,iy) = n
1979 DO n = 1, this%trans%poly%arraysize
1980 CALL l4f_category_log(this%category, l4f_debug, &
1981 'Polygon: '//t2c(n)//
' grid points: '// &
1982 t2c(count(this%inter_index_x(:,:) == n)))
1989 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
1993 CALL get_val(lin, nx=this%innx, ny=this%inny)
1994 this%outnx =
SIZE(v7d_out%ana)
1997 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1998 this%inter_index_y(this%outnx,this%outny))
1999 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
2001 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2003 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2004 CALL griddim_set_central_lon(lin, lonref)
2006 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2008 CALL this%find_index(lin, .true., &
2009 this%innx, this%inny, xmin, xmax, ymin, ymax, &
2010 lon, lat, this%trans%extrap, &
2011 this%inter_index_x, this%inter_index_y)
2014 nr = int(this%trans%area_info%radius)
2017 r2 = this%trans%area_info%radius**2
2018 ALLOCATE(this%stencil(n,n))
2019 this%stencil(:,:) = .true.
2022 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
2029 xnmax = this%innx - nr
2031 ynmax = this%inny - nr
2032 DO iy = 1, this%outny
2033 DO ix = 1, this%outnx
2034 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2035 this%inter_index_x(ix,iy) > xnmax .OR. &
2036 this%inter_index_y(ix,iy) < ynmin .OR. &
2037 this%inter_index_y(ix,iy) > ynmax)
THEN
2038 this%inter_index_x(ix,iy) = imiss
2039 this%inter_index_y(ix,iy) = imiss
2045 CALL l4f_category_log(this%category, l4f_debug, &
2046 'stencilinter: stencil size '//t2c(n*n))
2047 CALL l4f_category_log(this%category, l4f_debug, &
2048 'stencilinter: stencil points '//t2c(count(this%stencil)))
2056 ELSE IF (this%trans%trans_type ==
'maskinter')
THEN
2058 IF (.NOT.
PRESENT(maskgrid))
THEN
2059 CALL l4f_category_log(this%category,l4f_error, &
2060 'grid_transform_init maskgrid argument missing for maskinter transformation')
2061 CALL raise_fatal_error()
2064 CALL get_val(in, nx=this%innx, ny=this%inny)
2065 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2066 CALL l4f_category_log(this%category,l4f_error, &
2067 'grid_transform_init mask non conformal with input field')
2068 CALL l4f_category_log(this%category,l4f_error, &
2069 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2070 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2071 CALL raise_fatal_error()
2074 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2075 this%inter_index_y(this%innx,this%inny))
2076 this%inter_index_x(:,:) = imiss
2077 this%inter_index_y(:,:) = 1
2080 CALL gen_mask_class()
2088 DO iy = 1, this%inny
2089 DO ix = 1, this%innx
2090 IF (
c_e(maskgrid(ix,iy)))
THEN
2091 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2092 DO n = nmaskarea, 1, -1
2093 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2094 this%inter_index_x(ix,iy) = n
2104 this%outnx = nmaskarea
2107 CALL init(v7d_out, time_definition=time_definition)
2108 CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2113 CALL init(v7d_out%ana(n), &
2115 mask=(this%inter_index_x(:,:) == n))), &
2117 mask=(this%inter_index_x(:,:) == n))))
2123 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
2130 CALL get_val(in, nx=this%innx, ny=this%inny)
2132 ALLOCATE(this%point_index(this%innx,this%inny))
2133 this%point_index(:,:) = imiss
2136 CALL init(v7d_out, time_definition=time_definition)
2138 IF (this%trans%sub_type ==
'all' )
THEN
2140 this%outnx = this%innx*this%inny
2142 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2147 CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2148 lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2150 this%point_index(ix,iy) = n
2156 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN
2161 DO iy = 1, this%inny
2162 DO ix = 1, this%innx
2164 IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2165 lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2166 lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2167 lin%dim%lat(ix,iy) < this%trans%rect_coo%flat)
THEN
2168 this%outnx = this%outnx + 1
2169 this%point_index(ix,iy) = this%outnx
2174 IF (this%outnx <= 0)
THEN
2175 CALL l4f_category_log(this%category,l4f_warn, &
2176 "metamorphosis:coordbb: no points inside bounding box "//&
2177 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2178 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2179 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2180 trim(
to_char(this%trans%rect_coo%flat)))
2183 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2187 DO iy = 1, this%inny
2188 DO ix = 1, this%innx
2189 IF (
c_e(this%point_index(ix,iy)))
THEN
2191 CALL init(v7d_out%ana(n), &
2192 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2199 ELSE IF (this%trans%sub_type ==
'poly' )
THEN
2208 DO iy = 1, this%inny
2209 DO ix = 1, this%innx
2210 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2211 DO n = 1, this%trans%poly%arraysize
2212 IF (
inside(point, this%trans%poly%array(n)))
THEN
2213 this%outnx = this%outnx + 1
2214 this%point_index(ix,iy) = n
2223 IF (this%outnx <= 0)
THEN
2224 CALL l4f_category_log(this%category,l4f_warn, &
2225 "metamorphosis:poly: no points inside polygons")
2228 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2231 DO iy = 1, this%inny
2232 DO ix = 1, this%innx
2233 IF (
c_e(this%point_index(ix,iy)))
THEN
2235 CALL init(v7d_out%ana(n), &
2236 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2243 ELSE IF (this%trans%sub_type ==
'mask' )
THEN
2245 IF (.NOT.
PRESENT(maskgrid))
THEN
2246 CALL l4f_category_log(this%category,l4f_error, &
2247 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2252 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2253 CALL l4f_category_log(this%category,l4f_error, &
2254 'grid_transform_init mask non conformal with input field')
2255 CALL l4f_category_log(this%category,l4f_error, &
2256 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2257 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2263 CALL gen_mask_class()
2271 DO iy = 1, this%inny
2272 DO ix = 1, this%innx
2273 IF (
c_e(maskgrid(ix,iy)))
THEN
2274 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2275 DO n = nmaskarea, 1, -1
2276 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2277 this%outnx = this%outnx + 1
2278 this%point_index(ix,iy) = n
2288 IF (this%outnx <= 0)
THEN
2289 CALL l4f_category_log(this%category,l4f_warn, &
2290 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2294 CALL l4f_category_log(this%category,l4f_info, &
2295 "points in subarea "//t2c(n)//
": "// &
2296 t2c(count(this%point_index(:,:) == n)))
2299 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2302 DO iy = 1, this%inny
2303 DO ix = 1, this%innx
2304 IF (
c_e(this%point_index(ix,iy)))
THEN
2306 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2319 SUBROUTINE gen_mask_class()
2320 REAL :: startmaskclass, mmin, mmax
2323 IF (
PRESENT(maskbounds))
THEN
2324 nmaskarea =
SIZE(maskbounds) - 1
2325 IF (nmaskarea > 0)
THEN
2326 lmaskbounds = maskbounds
2328 ELSE IF (nmaskarea == 0)
THEN
2329 CALL l4f_category_log(this%category,l4f_warn, &
2330 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2331 //
', it will be ignored')
2332 CALL l4f_category_log(this%category,l4f_warn, &
2333 'at least 2 values are required for maskbounds')
2337 mmin = minval(maskgrid, mask=
c_e(maskgrid))
2338 mmax = maxval(maskgrid, mask=
c_e(maskgrid))
2340 nmaskarea = int(mmax - mmin + 1.5)
2341 startmaskclass = mmin - 0.5
2343 ALLOCATE(lmaskbounds(nmaskarea+1))
2344 lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2346 CALL l4f_category_log(this%category,l4f_debug, &
2347 'maskinter, '//t2c(nmaskarea)//
' subareas defined, with boundaries:')
2348 DO i = 1,
SIZE(lmaskbounds)
2349 CALL l4f_category_log(this%category,l4f_debug, &
2350 'maskinter '//t2c(i)//
' '//t2c(lmaskbounds(i)))
2354 END SUBROUTINE gen_mask_class
2356 END SUBROUTINE grid_transform_grid_vol7d_init
2377 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2380 TYPE(
vol7d),
INTENT(in) :: v7d_in
2382 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2385 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2386 DOUBLE PRECISION,
ALLOCATABLE :: lon(:,:),lat(:,:)
2390 CALL grid_transform_init_common(this, trans, categoryappend)
2392 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-vg6d")
2395 IF (this%trans%trans_type ==
'inter')
THEN
2397 IF ( this%trans%sub_type ==
'linear' )
THEN
2399 this%innx=
SIZE(v7d_in%ana)
2401 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2402 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2403 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2405 CALL copy (out, lout)
2407 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2408 CALL griddim_set_central_lon(lout, lonref)
2410 CALL get_val(lout, nx=nx, ny=ny)
2413 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2415 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2416 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2417 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2426 ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
2428 this%innx=
SIZE(v7d_in%ana)
2431 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2432 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2433 this%inter_index_y(this%innx,this%inny))
2435 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2437 CALL copy (out, lout)
2439 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2440 CALL griddim_set_central_lon(lout, lonref)
2442 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2443 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2446 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
2447 CALL get_val(out, dx=this%trans%area_info%boxdx)
2448 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
2449 CALL get_val(out, dx=this%trans%area_info%boxdy)
2451 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2452 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2455 CALL this%find_index(lout, .true., &
2456 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2457 lon, lat, .false., &
2458 this%inter_index_x, this%inter_index_y)
2467 END SUBROUTINE grid_transform_vol7d_grid_init
2504 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2505 maskbounds, categoryappend)
2508 TYPE(
vol7d),
INTENT(in) :: v7d_in
2509 TYPE(
vol7d),
INTENT(inout) :: v7d_out
2510 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2511 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2514 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2516 DOUBLE PRECISION :: lon1, lat1
2520 CALL grid_transform_init_common(this, trans, categoryappend)
2522 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-v7d")
2525 IF (this%trans%trans_type ==
'inter')
THEN
2527 IF ( this%trans%sub_type ==
'linear' )
THEN
2529 CALL vol7d_alloc_vol(v7d_out)
2530 this%outnx=
SIZE(v7d_out%ana)
2533 this%innx=
SIZE(v7d_in%ana)
2536 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2537 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2539 CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2540 CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2546 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
2548 this%innx=
SIZE(v7d_in%ana)
2551 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2552 this%inter_index_y(this%innx,this%inny))
2553 this%inter_index_x(:,:) = imiss
2554 this%inter_index_y(:,:) = 1
2556 DO i = 1,
SIZE(v7d_in%ana)
2558 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2559 point = georef_coord_new(x=lon1, y=lat1)
2561 DO n = 1, this%trans%poly%arraysize
2562 IF (
inside(point, this%trans%poly%array(n)))
THEN
2563 this%inter_index_x(i,1) = n
2569 this%outnx=this%trans%poly%arraysize
2571 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2575 CALL poly_to_coordinates(this%trans%poly, v7d_out)
2579 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
2582 this%innx =
SIZE(v7d_in%ana)
2585 ALLOCATE(this%point_index(this%innx,this%inny))
2586 this%point_index(:,:) = imiss
2588 IF (this%trans%sub_type ==
'all' )
THEN
2590 CALL metamorphosis_all_setup()
2592 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN
2594 ALLOCATE(lon(this%innx),lat(this%innx))
2599 CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2602 IF (lon(i) > this%trans%rect_coo%ilon .AND. &
2603 lon(i) < this%trans%rect_coo%flon .AND. &
2604 lat(i) > this%trans%rect_coo%ilat .AND. &
2605 lat(i) < this%trans%rect_coo%flat)
THEN
2606 this%outnx = this%outnx + 1
2607 this%point_index(i,1) = this%outnx
2611 IF (this%outnx <= 0)
THEN
2612 CALL l4f_category_log(this%category,l4f_warn, &
2613 "metamorphosis:coordbb: no points inside bounding box "//&
2614 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2615 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2616 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2617 trim(
to_char(this%trans%rect_coo%flat)))
2620 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2625 IF (
c_e(this%point_index(i,1)))
THEN
2627 CALL init(v7d_out%ana(n),lon=lon(i),lat=lat(i))
2630 DEALLOCATE(lon, lat)
2634 ELSE IF (this%trans%sub_type ==
'poly' )
THEN
2641 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2642 point = georef_coord_new(x=lon1, y=lat1)
2643 DO n = 1, this%trans%poly%arraysize
2644 IF (
inside(point, this%trans%poly%array(n)))
THEN
2645 this%outnx = this%outnx + 1
2646 this%point_index(i,1) = n
2653 IF (this%outnx <= 0)
THEN
2654 CALL l4f_category_log(this%category,l4f_warn, &
2655 "metamorphosis:poly: no points inside polygons")
2658 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2663 IF (
c_e(this%point_index(i,1)))
THEN
2666 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2667 CALL init(v7d_out%ana(n),lon=lon1,lat=lat1)
2673 ELSE IF (this%trans%sub_type ==
'setinvalidto' )
THEN
2675 IF (.NOT.
PRESENT(maskbounds))
THEN
2676 CALL l4f_category_log(this%category,l4f_error, &
2677 'grid_transform_init maskbounds missing for metamorphosis:'// &
2678 trim(this%trans%sub_type)//
' transformation')
2680 ELSE IF (
SIZE(maskbounds) < 1)
THEN
2681 CALL l4f_category_log(this%category,l4f_error, &
2682 'grid_transform_init maskbounds empty for metamorphosis:'// &
2683 trim(this%trans%sub_type)//
' transformation')
2686 this%val1 = maskbounds(1)
2688 CALL l4f_category_log(this%category, l4f_debug, &
2689 "grid_transform_init setting invalid data to "//t2c(this%val1))
2693 CALL metamorphosis_all_setup()
2695 ELSE IF (this%trans%sub_type ==
'settoinvalid' )
THEN
2697 IF (.NOT.
PRESENT(maskbounds))
THEN
2698 CALL l4f_category_log(this%category,l4f_error, &
2699 'grid_transform_init maskbounds missing for metamorphosis:'// &
2700 trim(this%trans%sub_type)//
' transformation')
2703 ELSE IF (
SIZE(maskbounds) < 2)
THEN
2704 CALL l4f_category_log(this%category,l4f_error, &
2705 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
2706 trim(this%trans%sub_type)//
' transformation')
2710 this%val1 = maskbounds(1)
2711 this%val2 = maskbounds(
SIZE(maskbounds))
2713 CALL l4f_category_log(this%category, l4f_debug, &
2714 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
2715 t2c(this%val2)//
']')
2719 CALL metamorphosis_all_setup()
2728 SUBROUTINE metamorphosis_all_setup()
2730 this%outnx =
SIZE(v7d_in%ana)
2732 this%point_index(:,1) = (/(i,i=1,this%innx)/)
2733 CALL vol7d_alloc(v7d_out, nana=
SIZE(v7d_in%ana))
2734 v7d_out%ana = v7d_in%ana
2738 END SUBROUTINE metamorphosis_all_setup
2740 END SUBROUTINE grid_transform_vol7d_vol7d_init
2744 SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
2747 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2749 CHARACTER(len=512) :: a_name
2751 IF (
PRESENT(categoryappend))
THEN
2752 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
2753 trim(categoryappend))
2755 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2757 this%category=l4f_category_get(a_name)
2760 CALL l4f_category_log(this%category,l4f_debug,
"start init_grid_transform")
2765 END SUBROUTINE grid_transform_init_common
2769 SUBROUTINE poly_to_coordinates(poly, v7d_out)
2771 TYPE(
vol7d),
INTENT(inout) :: v7d_out
2774 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2776 DO n = 1, poly%arraysize
2777 CALL getval(poly%array(n), x=lon, y=lat)
2778 sz = min(
SIZE(lon),
SIZE(lat))
2779 IF (lon(1) == lon(sz) .AND. lat(1) == lat(sz))
THEN
2785 END SUBROUTINE poly_to_coordinates
2790 SUBROUTINE grid_transform_delete(this)
2808 if (
associated(this%inter_index_x))
deallocate (this%inter_index_x)
2809 if (
associated(this%inter_index_y))
deallocate (this%inter_index_y)
2810 if (
associated(this%inter_index_z))
deallocate (this%inter_index_z)
2811 if (
associated(this%point_index))
deallocate (this%point_index)
2813 if (
associated(this%inter_x))
deallocate (this%inter_x)
2814 if (
associated(this%inter_y))
deallocate (this%inter_y)
2816 if (
associated(this%inter_xp))
deallocate (this%inter_xp)
2817 if (
associated(this%inter_yp))
deallocate (this%inter_yp)
2818 if (
associated(this%inter_zp))
deallocate (this%inter_zp)
2819 if (
associated(this%vcoord_in))
deallocate (this%vcoord_in)
2820 if (
associated(this%vcoord_out))
deallocate (this%vcoord_out)
2821 if (
associated(this%point_mask))
deallocate (this%point_mask)
2822 if (
associated(this%stencil))
deallocate (this%stencil)
2823 if (
associated(this%output_level_auto))
deallocate (this%output_level_auto)
2824 IF (
ALLOCATED(this%coord_3d_in))
DEALLOCATE(this%coord_3d_in)
2825 this%valid = .false.
2828 call l4f_category_delete(this%category)
2830 END SUBROUTINE grid_transform_delete
2837 SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
2838 point_index, output_point_index, levshift, levused)
2840 TYPE(vol7d_level),
POINTER,
OPTIONAL :: output_level_auto(:)
2841 LOGICAL,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_mask(:)
2842 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_index(:)
2843 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: output_point_index(:)
2844 INTEGER,
INTENT(out),
OPTIONAL :: levshift
2845 INTEGER,
INTENT(out),
OPTIONAL :: levused
2849 IF (
PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
2850 IF (
PRESENT(point_mask))
THEN
2851 IF (
ASSOCIATED(this%point_index))
THEN
2852 point_mask =
c_e(reshape(this%point_index, (/
SIZE(this%point_index)/)))
2855 IF (
PRESENT(point_index))
THEN
2856 IF (
ASSOCIATED(this%point_index))
THEN
2857 point_index = reshape(this%point_index, (/
SIZE(this%point_index)/))
2860 IF (
PRESENT(output_point_index))
THEN
2861 IF (
ASSOCIATED(this%point_index))
THEN
2863 output_point_index = pack(this%point_index(:,:),
c_e(this%point_index))
2864 ELSE IF (this%trans%trans_type ==
'polyinter' .OR. &
2865 this%trans%trans_type ==
'maskinter')
THEN
2867 output_point_index = (/(i,i=1,this%outnx)/)
2870 IF (
PRESENT(levshift)) levshift = this%levshift
2871 IF (
PRESENT(levused)) levused = this%levused
2873 END SUBROUTINE grid_transform_get_val
2878 FUNCTION grid_transform_c_e(this)
2880 LOGICAL :: grid_transform_c_e
2882 grid_transform_c_e = this%valid
2884 END FUNCTION grid_transform_c_e
2896 RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
2899 REAL,
INTENT(in) :: field_in(:,:,:)
2900 REAL,
INTENT(out) :: field_out(:,:,:)
2901 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
2902 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
2904 INTEGER :: i, j, k, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
2905 kfound, kfoundin, inused, i1, i2, j1, j2, np, ns
2906 INTEGER,
ALLOCATABLE :: nval(:,:)
2907 REAL :: z1,z2,z3,z4,z(4)
2908 DOUBLE PRECISION :: x1,x3,y1,y3,xp,yp
2909 INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype
2910 REAL,
ALLOCATABLE :: coord_in(:)
2911 LOGICAL,
ALLOCATABLE :: mask_in(:)
2912 REAL,
ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
2913 REAL,
POINTER :: coord_3d_in_act(:,:,:)
2915 LOGICAL :: alloc_coord_3d_in_act, nm1
2919 CALL l4f_category_log(this%category,l4f_debug,
"start grid_transform_compute")
2922 field_out(:,:,:) = rmiss
2924 IF (.NOT.this%valid)
THEN
2925 CALL l4f_category_log(this%category,l4f_error, &
2926 "refusing to perform a non valid transformation")
2930 IF (this%recur)
THEN
2932 CALL l4f_category_log(this%category,l4f_debug, &
2933 "recursive transformation in grid_transform_compute")
2936 IF (this%trans%trans_type ==
'polyinter')
THEN
2938 likethis%recur = .false.
2939 likethis%outnx = this%trans%poly%arraysize
2941 ALLOCATE(field_tmp(this%trans%poly%arraysize,1,
SIZE(field_out,3)))
2942 CALL grid_transform_compute(likethis, field_in, field_tmp, var)
2944 DO k = 1,
SIZE(field_out,3)
2947 IF (
c_e(this%inter_index_x(i,j)))
THEN
2948 field_out(i,j,k) = field_tmp(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2953 DEALLOCATE(field_tmp)
2960 IF (
PRESENT(var))
THEN
2961 vartype = vol7d_vartype(var)
2964 innx =
SIZE(field_in,1); inny =
SIZE(field_in,2); innz =
SIZE(field_in,3)
2965 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
2968 IF (this%trans%trans_type ==
'vertint')
THEN
2970 IF (innz /= this%innz)
THEN
2971 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2972 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
2973 t2c(this%innz)//
" /= "//t2c(innz))
2978 IF (outnz /= this%outnz)
THEN
2979 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2980 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
2981 t2c(this%outnz)//
" /= "//t2c(outnz))
2986 IF (innx /= outnx .OR. inny /= outny)
THEN
2987 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2988 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
2989 t2c(innx)//
","//t2c(inny)//
" /= "//&
2990 t2c(outnx)//
","//t2c(outny))
2997 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
2998 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
2999 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
3000 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
3001 t2c(innx)//
","//t2c(inny))
3006 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
3007 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3008 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
3009 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
3010 t2c(outnx)//
","//t2c(outny))
3015 IF (innz /= outnz)
THEN
3016 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3017 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
3018 t2c(innz)//
" /= "//t2c(outnz))
3026 call l4f_category_log(this%category,l4f_debug, &
3027 "start grid_transform_compute "//trim(this%trans%trans_type)//
':'// &
3028 trim(this%trans%sub_type))
3031 IF (this%trans%trans_type ==
'zoom')
THEN
3033 field_out(this%outinx:this%outfnx, &
3034 this%outiny:this%outfny,:) = &
3035 field_in(this%iniox:this%infox, &
3036 this%inioy:this%infoy,:)
3038 ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN
3040 IF (this%trans%sub_type ==
'average')
THEN
3041 IF (vartype == var_dir360)
THEN
3044 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3045 je = j+this%trans%box_info%npy-1
3048 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3049 ie = i+this%trans%box_info%npx-1
3051 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3053 field_out(ii,jj,k) = find_prevailing_direction(field_in(i:ie,j:je,k), &
3063 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3064 je = j+this%trans%box_info%npy-1
3067 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3068 ie = i+this%trans%box_info%npx-1
3070 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3072 field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), &
3073 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3081 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3082 this%trans%sub_type ==
'stddevnm1')
THEN
3084 IF (this%trans%sub_type ==
'stddev')
THEN
3090 navg = this%trans%box_info%npx*this%trans%box_info%npy
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
3101 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3106 ELSE IF (this%trans%sub_type ==
'max')
THEN
3109 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3110 je = j+this%trans%box_info%npy-1
3113 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3114 ie = i+this%trans%box_info%npx-1
3116 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3118 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3119 mask=(field_in(i:ie,j:je,k) /= rmiss))
3125 ELSE IF (this%trans%sub_type ==
'min')
THEN
3128 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3129 je = j+this%trans%box_info%npy-1
3132 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3133 ie = i+this%trans%box_info%npx-1
3135 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3137 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3138 mask=(field_in(i:ie,j:je,k) /= rmiss))
3144 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3146 navg = this%trans%box_info%npx*this%trans%box_info%npy
3149 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3150 je = j+this%trans%box_info%npy-1
3153 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3154 ie = i+this%trans%box_info%npx-1
3157 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3158 (/real(this%trans%stat_info%percentile)/))
3163 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3167 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3168 je = j+this%trans%box_info%npy-1
3171 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3172 ie = i+this%trans%box_info%npx-1
3174 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3176 field_out(ii,jj,k) = sum(interval_info_valid( &
3177 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3178 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3186 ELSE IF (this%trans%trans_type ==
'inter')
THEN
3188 IF (this%trans%sub_type ==
'near')
THEN
3191 DO j = 1, this%outny
3192 DO i = 1, this%outnx
3194 IF (
c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3195 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3201 ELSE IF (this%trans%sub_type ==
'bilin')
THEN
3204 DO j = 1, this%outny
3205 DO i = 1, this%outnx
3207 IF (
c_e(this%inter_index_x(i,j)))
THEN
3209 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3210 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3211 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3212 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3214 IF (
c_e(z1) .AND.
c_e(z2) .AND.
c_e(z3) .AND.
c_e(z4))
THEN
3216 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3217 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3218 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3219 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3221 xp=this%inter_xp(i,j)
3222 yp=this%inter_yp(i,j)
3224 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3232 ELSE IF (this%trans%sub_type ==
'shapiro_near')
THEN
3234 DO j = 1, this%outny
3235 DO i = 1, this%outnx
3237 IF (
c_e(this%inter_index_x(i,j)))
THEN
3239 IF(this%inter_index_x(i,j)-1>0)
THEN
3240 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3244 IF(this%inter_index_x(i,j)+1<=this%outnx)
THEN
3245 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3249 IF(this%inter_index_y(i,j)+1<=this%outny)
THEN
3250 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3254 IF(this%inter_index_y(i,j)-1>0)
THEN
3255 z(4)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)-1,k)
3259 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3268 ELSE IF (this%trans%trans_type ==
'boxinter' &
3269 .OR. this%trans%trans_type ==
'polyinter' &
3270 .OR. this%trans%trans_type ==
'maskinter')
THEN
3272 IF (this%trans%sub_type ==
'average')
THEN
3274 IF (vartype == var_dir360)
THEN
3276 DO j = 1, this%outny
3277 DO i = 1, this%outnx
3278 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3280 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3286 ALLOCATE(nval(this%outnx, this%outny))
3287 field_out(:,:,:) = 0.0
3292 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3293 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3294 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3296 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3297 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3301 WHERE (nval(:,:) /= 0)
3302 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3304 field_out(:,:,k) = rmiss
3310 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3311 this%trans%sub_type ==
'stddevnm1')
THEN
3313 IF (this%trans%sub_type ==
'stddev')
THEN
3319 DO j = 1, this%outny
3320 DO i = 1, this%outnx
3323 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3324 mask=reshape((this%inter_index_x == i .AND. &
3325 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)), nm1=nm1)
3330 ELSE IF (this%trans%sub_type ==
'max')
THEN
3335 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3336 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3337 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3338 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3341 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3350 ELSE IF (this%trans%sub_type ==
'min')
THEN
3355 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3356 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3357 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3358 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3361 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3369 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3372 DO j = 1, this%outny
3373 DO i = 1, this%outnx
3376 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3377 (/real(this%trans%stat_info%percentile)/), &
3378 mask=reshape((this%inter_index_x == i .AND. &
3379 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)))
3384 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3386 ALLOCATE(nval(this%outnx, this%outny))
3387 field_out(:,:,:) = 0.0
3392 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3393 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3394 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3395 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3396 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3397 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3401 WHERE (nval(:,:) /= 0)
3402 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3404 field_out(:,:,k) = rmiss
3411 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
3412 np =
SIZE(this%stencil,1)/2
3413 ns =
SIZE(this%stencil)
3415 IF (this%trans%sub_type ==
'average')
THEN
3417 IF (vartype == var_dir360)
THEN
3419 DO j = 1, this%outny
3420 DO i = 1, this%outnx
3421 IF (
c_e(this%inter_index_x(i,j)))
THEN
3422 i1 = this%inter_index_x(i,j) - np
3423 i2 = this%inter_index_x(i,j) + np
3424 j1 = this%inter_index_y(i,j) - np
3425 j2 = this%inter_index_y(i,j) + np
3426 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3428 mask=this%stencil(:,:))
3439 DO j = 1, this%outny
3440 DO i = 1, this%outnx
3441 IF (
c_e(this%inter_index_x(i,j)))
THEN
3442 i1 = this%inter_index_x(i,j) - np
3443 i2 = this%inter_index_x(i,j) + np
3444 j1 = this%inter_index_y(i,j) - np
3445 j2 = this%inter_index_y(i,j) + np
3446 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3448 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3449 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3458 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3459 this%trans%sub_type ==
'stddevnm1')
THEN
3461 IF (this%trans%sub_type ==
'stddev')
THEN
3470 DO j = 1, this%outny
3471 DO i = 1, this%outnx
3472 IF (
c_e(this%inter_index_x(i,j)))
THEN
3473 i1 = this%inter_index_x(i,j) - np
3474 i2 = this%inter_index_x(i,j) + np
3475 j1 = this%inter_index_y(i,j) - np
3476 j2 = this%inter_index_y(i,j) + np
3479 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3480 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3481 this%stencil(:,:), (/ns/)), nm1=nm1)
3488 ELSE IF (this%trans%sub_type ==
'max')
THEN
3493 DO j = 1, this%outny
3494 DO i = 1, this%outnx
3495 IF (
c_e(this%inter_index_x(i,j)))
THEN
3496 i1 = this%inter_index_x(i,j) - np
3497 i2 = this%inter_index_x(i,j) + np
3498 j1 = this%inter_index_y(i,j) - np
3499 j2 = this%inter_index_y(i,j) + np
3500 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3502 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3503 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3511 ELSE IF (this%trans%sub_type ==
'min')
THEN
3516 DO j = 1, this%outny
3517 DO i = 1, this%outnx
3518 IF (
c_e(this%inter_index_x(i,j)))
THEN
3519 i1 = this%inter_index_x(i,j) - np
3520 i2 = this%inter_index_x(i,j) + np
3521 j1 = this%inter_index_y(i,j) - np
3522 j2 = this%inter_index_y(i,j) + np
3523 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3525 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3526 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3534 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3539 DO j = 1, this%outny
3540 DO i = 1, this%outnx
3541 IF (
c_e(this%inter_index_x(i,j)))
THEN
3542 i1 = this%inter_index_x(i,j) - np
3543 i2 = this%inter_index_x(i,j) + np
3544 j1 = this%inter_index_y(i,j) - np
3545 j2 = this%inter_index_y(i,j) + np
3548 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3549 (/real(this%trans%stat_info%percentile)/), &
3550 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3551 this%stencil(:,:), (/ns/)))
3558 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3563 DO j = 1, this%outny
3564 DO i = 1, this%outnx
3565 IF (
c_e(this%inter_index_x(i,j)))
THEN
3566 i1 = this%inter_index_x(i,j) - np
3567 i2 = this%inter_index_x(i,j) + np
3568 j1 = this%inter_index_y(i,j) - np
3569 j2 = this%inter_index_y(i,j) + np
3570 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3572 field_out(i,j,k) = sum(interval_info_valid( &
3573 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3574 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3584 ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
3587 WHERE(
c_e(this%inter_index_x(:,:)))
3588 field_out(:,:,k) = real(this%inter_index_x(:,:))
3592 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
3594 IF (this%trans%sub_type ==
'all')
THEN
3596 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3598 ELSE IF (this%trans%sub_type ==
'coordbb' .OR. this%trans%sub_type ==
'poly' &
3599 .OR. this%trans%sub_type ==
'mask')
THEN
3603 field_out(:,1,k) = pack(field_in(:,:,k),
c_e(this%point_index(:,:)))
3606 ELSE IF (this%trans%sub_type ==
'maskvalid' .OR. &
3607 this%trans%sub_type ==
'maskinvalid')
THEN
3610 WHERE (this%point_mask(:,:))
3611 field_out(:,:,k) = field_in(:,:,k)
3615 ELSE IF (this%trans%sub_type ==
'lemaskinvalid')
THEN
3618 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3619 field_out(:,:,k) = field_in(:,:,k)
3621 field_out(:,:,k) = rmiss
3625 ELSE IF (this%trans%sub_type ==
'ltmaskinvalid')
THEN
3628 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3629 field_out(:,:,k) = field_in(:,:,k)
3631 field_out(:,:,k) = rmiss
3635 ELSE IF (this%trans%sub_type ==
'gemaskinvalid')
THEN
3638 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3639 field_out(:,:,k) = field_in(:,:,k)
3641 field_out(:,:,k) = rmiss
3645 ELSE IF (this%trans%sub_type ==
'gtmaskinvalid')
THEN
3648 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3649 field_out(:,:,k) = field_in(:,:,k)
3651 field_out(:,:,k) = rmiss
3655 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
3658 WHERE (
c_e(field_in(:,:,k)))
3659 field_out(:,:,k) = field_in(:,:,k)
3661 field_out(:,:,k) = this%val1
3665 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
3666 IF (
c_e(this%val1) .AND.
c_e(this%val2))
THEN
3667 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3668 .AND. field_in(:,:,:) <= this%val2)
3669 field_out(:,:,:) = rmiss
3671 field_out(:,:,:) = field_in(:,:,:)
3673 ELSE IF (
c_e(this%val1))
THEN
3674 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3675 field_out(:,:,:) = rmiss
3677 field_out(:,:,:) = field_in(:,:,:)
3679 ELSE IF (
c_e(this%val2))
THEN
3680 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3681 field_out(:,:,:) = rmiss
3683 field_out(:,:,:) = field_in(:,:,:)
3688 ELSE IF (this%trans%trans_type ==
'vertint')
THEN
3690 IF (this%trans%sub_type ==
'linear')
THEN
3692 alloc_coord_3d_in_act = .false.
3693 IF (
ASSOCIATED(this%inter_index_z))
THEN
3696 IF (
c_e(this%inter_index_z(k)))
THEN
3697 z1 = real(this%inter_zp(k))
3698 z2 = real(1.0d0 - this%inter_zp(k))
3699 SELECT CASE(vartype)
3704 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3705 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3706 field_out(i,j,k) = &
3707 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3708 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3716 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3717 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3718 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3719 field_in(i,j,this%inter_index_z(k)+1) > 0.)
THEN
3720 field_out(i,j,k) = exp( &
3721 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3722 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3730 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3731 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3732 field_out(i,j,k) = &
3733 field_in(i,j,this%inter_index_z(k))*z2 + &
3734 field_in(i,j,this%inter_index_z(k)+1)*z1
3746 IF (
ALLOCATED(this%coord_3d_in))
THEN
3747 coord_3d_in_act => this%coord_3d_in
3749 CALL l4f_category_log(this%category,l4f_debug, &
3750 "Using external vertical coordinate for vertint")
3753 IF (
PRESENT(coord_3d_in))
THEN
3755 CALL l4f_category_log(this%category,l4f_debug, &
3756 "Using internal vertical coordinate for vertint")
3758 IF (this%dolog)
THEN
3759 ALLOCATE(coord_3d_in_act(
SIZE(coord_3d_in,1), &
3760 SIZE(coord_3d_in,2),
SIZE(coord_3d_in,3)))
3761 alloc_coord_3d_in_act = .true.
3762 WHERE (
c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3763 coord_3d_in_act = log(coord_3d_in)
3765 coord_3d_in_act = rmiss
3768 coord_3d_in_act => coord_3d_in
3771 CALL l4f_category_log(this%category,l4f_error, &
3772 "Missing internal and external vertical coordinate for vertint")
3778 inused =
SIZE(coord_3d_in_act,3)
3779 IF (inused < 2)
RETURN
3783 IF (.NOT.
c_e(this%vcoord_out(k))) cycle
3787 DO kk = 1, max(inused-kkcache-1, kkcache)
3789 kkdown = kkcache - kk + 1
3791 IF (kkdown >= 1)
THEN
3792 IF (this%vcoord_out(k) >= &
3793 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3794 this%vcoord_out(k) < &
3795 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)))
THEN
3798 kfound = kkcache + this%levshift
3802 IF (kkup < inused)
THEN
3803 IF (this%vcoord_out(k) >= &
3804 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3805 this%vcoord_out(k) < &
3806 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)))
THEN
3809 kfound = kkcache + this%levshift
3816 IF (
c_e(kfound))
THEN
3817 IF (
c_e(field_in(i,j,kfound)) .AND.
c_e(field_in(i,j,kfound+1)))
THEN
3818 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3819 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3821 SELECT CASE(vartype)
3824 field_out(i,j,k) = &
3825 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3827 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3828 log(field_in(i,j,kfound+1))*z1)
3831 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3840 IF (alloc_coord_3d_in_act)
DEALLOCATE(coord_3d_in_act)
3843 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN
3846 IF (.NOT.
ASSOCIATED(this%vcoord_in) .AND. .NOT.
PRESENT(coord_3d_in))
THEN
3847 CALL l4f_category_log(this%category,l4f_error, &
3848 "linearsparse interpolation, no input vert coord available")
3852 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3856 IF (
ASSOCIATED(this%vcoord_in))
THEN
3857 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3858 .AND.
c_e(this%vcoord_in(:))
3860 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3861 .AND.
c_e(coord_3d_in(i,j,:))
3864 IF (vartype == var_press)
THEN
3865 mask_in(:) = mask_in(:) .AND. &
3866 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3868 inused = count(mask_in)
3869 IF (inused > 1)
THEN
3870 IF (
ASSOCIATED(this%vcoord_in))
THEN
3871 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3873 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3875 IF (vartype == var_press)
THEN
3876 val_in(1:inused) = log(pack( &
3877 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3880 val_in(1:inused) = pack( &
3881 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3888 DO kk = 1, max(inused-kkcache-1, kkcache)
3890 kkdown = kkcache - kk + 1
3892 IF (kkdown >= 1)
THEN
3893 IF (this%vcoord_out(k) >= &
3894 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
3895 this%vcoord_out(k) < &
3896 max(coord_in(kkdown), coord_in(kkdown+1)))
THEN
3903 IF (kkup < inused)
THEN
3904 IF (this%vcoord_out(k) >= &
3905 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
3906 this%vcoord_out(k) < &
3907 max(coord_in(kkup), coord_in(kkup+1)))
THEN
3917 IF (
c_e(kfound))
THEN
3918 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
3919 (coord_in(kfound) - coord_in(kfound-1)))
3921 IF (vartype == var_dir360)
THEN
3922 field_out(i,j,k) = &
3923 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
3924 ELSE IF (vartype == var_press)
THEN
3925 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
3927 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
3937 DEALLOCATE(coord_in,val_in)
3942 ELSE IF (this%trans%trans_type ==
'' .OR. this%trans%trans_type ==
'none')
THEN
3944 field_out(:,:,:) = field_in(:,:,:)
3954 FUNCTION interp_var_360(v1, v2, w1, w2)
3955 REAL,
INTENT(in) :: v1, v2, w1, w2
3956 REAL :: interp_var_360
3960 IF (abs(v1 - v2) > 180.)
THEN
3968 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
3970 interp_var_360 = v1*w2 + v2*w1
3973 END FUNCTION interp_var_360
3976 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask)
RESULT(prev)
3977 REAL,
INTENT(in) :: v1(:,:)
3978 REAL,
INTENT(in) :: l, h, res
3979 LOGICAL,
INTENT(in),
OPTIONAL :: mask(:,:)
3987 IF ((h - l) <= res)
THEN
3992 IF (
PRESENT(mask))
THEN
3993 nl = count(v1 >= l .AND. v1 < m .AND. mask)
3994 nh = count(v1 >= m .AND. v1 < h .AND. mask)
3996 nl = count(v1 >= l .AND. v1 < m)
3997 nh = count(v1 >= m .AND. v1 < h)
4000 prev = find_prevailing_direction(v1, m, h, res)
4001 ELSE IF (nl > nh)
THEN
4002 prev = find_prevailing_direction(v1, l, m, res)
4003 ELSE IF (nl == 0 .AND. nh == 0)
THEN
4009 END FUNCTION find_prevailing_direction
4012 END SUBROUTINE grid_transform_compute
4020 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4023 REAL,
INTENT(in) :: field_in(:,:)
4024 REAL,
INTENT(out):: field_out(:,:,:)
4025 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
4026 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
4028 real,
allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4029 INTEGER :: inn_p, ier, k
4030 INTEGER :: innx, inny, innz, outnx, outny, outnz
4033 call l4f_category_log(this%category,l4f_debug,
"start v7d_grid_transform_compute")
4036 field_out(:,:,:) = rmiss
4038 IF (.NOT.this%valid)
THEN
4039 call l4f_category_log(this%category,l4f_error, &
4040 "refusing to perform a non valid transformation")
4045 innx =
SIZE(field_in,1); inny = 1; innz =
SIZE(field_in,2)
4046 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
4049 IF (this%trans%trans_type ==
'vertint')
THEN
4051 IF (innz /= this%innz)
THEN
4052 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4053 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4054 t2c(this%innz)//
" /= "//t2c(innz))
4059 IF (outnz /= this%outnz)
THEN
4060 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4061 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4062 t2c(this%outnz)//
" /= "//t2c(outnz))
4067 IF (innx /= outnx .OR. inny /= outny)
THEN
4068 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4069 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
4070 t2c(innx)//
","//t2c(inny)//
" /= "//&
4071 t2c(outnx)//
","//t2c(outny))
4078 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
4079 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4080 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4081 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
4082 t2c(innx)//
","//t2c(inny))
4087 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
4088 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4089 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4090 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
4091 t2c(outnx)//
","//t2c(outny))
4096 IF (innz /= outnz)
THEN
4097 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4098 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
4099 t2c(innz)//
" /= "//t2c(outnz))
4107 call l4f_category_log(this%category,l4f_debug, &
4108 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//
':'// &
4109 trim(this%trans%sub_type))
4112 IF (this%trans%trans_type ==
'inter')
THEN
4114 IF (this%trans%sub_type ==
'linear')
THEN
4116 #ifdef HAVE_LIBNGMATH
4118 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4120 inn_p = count(
c_e(field_in(:,k)))
4122 CALL l4f_category_log(this%category,l4f_info, &
4123 "Number of sparse data points: "//t2c(inn_p)//
','//t2c(
SIZE(field_in(:,k))))
4127 field_in_p(1:inn_p) = pack(field_in(:,k),
c_e(field_in(:,k)))
4128 x_in_p(1:inn_p) = pack(this%inter_xp(:,1),
c_e(field_in(:,k)))
4129 y_in_p(1:inn_p) = pack(this%inter_yp(:,1),
c_e(field_in(:,k)))
4131 IF (.NOT.this%trans%extrap)
THEN
4132 CALL nnseti(
'ext', 0)
4133 CALL nnsetr(
'nul', rmiss)
4136 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4137 this%outnx, this%outny, real(this%inter_x(:,1)), &
4138 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4141 CALL l4f_category_log(this%category,l4f_error, &
4142 "Error code from NCAR natgrids: "//t2c(ier))
4148 CALL l4f_category_log(this%category,l4f_info, &
4149 "insufficient data in gridded region to triangulate")
4153 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4156 CALL l4f_category_log(this%category,l4f_error, &
4157 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4164 ELSE IF (this%trans%trans_type ==
'boxinter' .OR. &
4165 this%trans%trans_type ==
'polyinter' .OR. &
4166 this%trans%trans_type ==
'vertint' .OR. &
4167 this%trans%trans_type ==
'metamorphosis')
THEN
4170 reshape(field_in, (/
SIZE(field_in,1), 1,
SIZE(field_in,2)/)), field_out, var, &
4175 END SUBROUTINE grid_transform_v7d_grid_compute
4190 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
RESULT(zp)
4191 REAL,
INTENT(in) :: z1,z2,z3,z4
4192 DOUBLE PRECISION,
INTENT(in):: x1,y1
4193 DOUBLE PRECISION,
INTENT(in):: x3,y3
4194 DOUBLE PRECISION,
INTENT(in):: xp,yp
4201 p2 = real((yp-y1)/(y3-y1))
4202 p1 = real((xp-x1)/(x3-x1))
4207 zp = (z6-z5)*(p1)+z5
4213 FUNCTION shapiro (z,zp)
RESULT(zs)
4216 REAL,
INTENT(in) :: z(:)
4217 REAL,
INTENT(in) :: zp
4223 zs = zp - 0.5* ( n*zp - sum(z,
c_e(z)) )/n
4228 END FUNCTION shapiro
4232 SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
4233 lon, lat, extrap, index_x, index_y)
4235 logical,
INTENT(in) :: near
4236 INTEGER,
INTENT(in) :: nx,ny
4237 DOUBLE PRECISION,
INTENT(in) :: xmin, xmax, ymin, ymax
4238 DOUBLE PRECISION,
INTENT(in) :: lon(:,:),lat(:,:)
4239 LOGICAL,
INTENT(in) :: extrap
4240 INTEGER,
INTENT(out) :: index_x(:,:),index_y(:,:)
4243 DOUBLE PRECISION :: x(SIZE(lon,1),SIZE(lon,2)),y(SIZE(lon,1),SIZE(lon,2))
4246 CALL proj(this,lon,lat,x,y)
4247 index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4248 index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4252 CALL proj(this,lon,lat,x,y)
4253 index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4254 index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4260 index_x = max(index_x, 1)
4261 index_y = max(index_y, 1)
4262 index_x = min(index_x, lnx)
4263 index_y = min(index_y, lny)
4265 WHERE(index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
4271 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...