220CHARACTER(len=255),
PARAMETER:: subcategory=
"grid_transform_class"
224 double precision :: boxdx
225 double precision :: boxdy
226 double precision :: radius
231 DOUBLE PRECISION :: percentile
240END 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
380PUBLIC interval_info, interval_info_new, interval_info_valid, basic_find_index
385FUNCTION interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
RESULT(this)
386REAL,
INTENT(in),
OPTIONAL :: interv_gt
387REAL,
INTENT(in),
OPTIONAL :: interv_ge
388REAL,
INTENT(in),
OPTIONAL :: interv_lt
389REAL,
INTENT(in),
OPTIONAL :: interv_le
391TYPE(interval_info) :: this
393IF (
PRESENT(interv_gt))
THEN
394 IF (
c_e(interv_gt))
THEN
399IF (
PRESENT(interv_ge))
THEN
400 IF (
c_e(interv_ge))
THEN
405IF (
PRESENT(interv_lt))
THEN
406 IF (
c_e(interv_lt))
THEN
411IF (
PRESENT(interv_le))
THEN
412 IF (
c_e(interv_le))
THEN
418END FUNCTION interval_info_new
425ELEMENTAL FUNCTION interval_info_valid(this, val)
426TYPE(interval_info),
INTENT(in) :: this
427REAL,
INTENT(in) :: val
429REAL :: interval_info_valid
431interval_info_valid = 1.0
433IF (
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
439IF (
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
446END FUNCTION interval_info_valid
454SUBROUTINE 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)
460TYPE(transform_def),
INTENT(out) :: this
461CHARACTER(len=*) :: trans_type
462CHARACTER(len=*) :: sub_type
463INTEGER,
INTENT(in),
OPTIONAL :: ix
464INTEGER,
INTENT(in),
OPTIONAL :: iy
465INTEGER,
INTENT(in),
OPTIONAL :: fx
466INTEGER,
INTENT(in),
OPTIONAL :: fy
467DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilon
468DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilat
469DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flon
470DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flat
471INTEGER,
INTENT(IN),
OPTIONAL :: npx
472INTEGER,
INTENT(IN),
OPTIONAL :: npy
473DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdx
474DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdy
475DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: radius
476TYPE(arrayof_georef_coord_array),
OPTIONAL :: poly
477DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: percentile
478REAL,
INTENT(in),
OPTIONAL :: interv_gt
479REAL,
INTENT(in),
OPTIONAL :: interv_ge
480REAL,
INTENT(in),
OPTIONAL :: interv_lt
481REAL,
INTENT(in),
OPTIONAL :: interv_le
482LOGICAL,
INTENT(IN),
OPTIONAL :: extrap
483INTEGER,
INTENT(IN),
OPTIONAL :: time_definition
484TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: input_levtype
485TYPE(vol7d_var),
INTENT(IN),
OPTIONAL :: input_coordvar
486TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: output_levtype
487character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
489character(len=512) :: a_name
491IF (
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))
497this%category=l4f_category_get(a_name)
499this%trans_type = trans_type
500this%sub_type = sub_type
502CALL optio(extrap,this%extrap)
504call optio(ix,this%rect_ind%ix)
505call optio(iy,this%rect_ind%iy)
506call optio(fx,this%rect_ind%fx)
507call optio(fy,this%rect_ind%fy)
509call optio(ilon,this%rect_coo%ilon)
510call optio(ilat,this%rect_coo%ilat)
511call optio(flon,this%rect_coo%flon)
512call optio(flat,this%rect_coo%flat)
514CALL optio(boxdx,this%area_info%boxdx)
515CALL optio(boxdy,this%area_info%boxdy)
516CALL optio(radius,this%area_info%radius)
517IF (
PRESENT(poly)) this%poly = poly
518CALL optio(percentile,this%stat_info%percentile)
520this%interval_info = interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
522CALL optio(npx,this%box_info%npx)
523CALL optio(npy,this%box_info%npy)
525IF (
PRESENT(input_levtype))
THEN
526 this%vertint%input_levtype = input_levtype
528 this%vertint%input_levtype = vol7d_level_miss
530IF (
PRESENT(input_coordvar))
THEN
531 this%vertint%input_coordvar = input_coordvar
533 this%vertint%input_coordvar = vol7d_var_miss
535IF (
PRESENT(output_levtype))
THEN
536 this%vertint%output_levtype = output_levtype
538 this%vertint%output_levtype = vol7d_level_miss
541call optio(time_definition,this%time_definition)
542if (
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()
549IF (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()
622ELSE 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()
632ELSE 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()
694ELSE 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()
711ELSE 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()
732ELSE 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()
772SUBROUTINE sub_type_error()
774CALL l4f_category_log(this%category, l4f_error, trim(this%trans_type) &
775 //
': sub_type '//trim(this%sub_type)//
' is not defined')
776CALL raise_fatal_error()
778END SUBROUTINE sub_type_error
780SUBROUTINE trans_type_error()
782CALL l4f_category_log(this%category, l4f_error,
'trans_type '//this%trans_type &
784CALL raise_fatal_error()
786END SUBROUTINE trans_type_error
789END SUBROUTINE transform_init
795SUBROUTINE transform_delete(this)
801this%rect_ind%ix=imiss
802this%rect_ind%iy=imiss
803this%rect_ind%fx=imiss
804this%rect_ind%fy=imiss
806this%rect_coo%ilon=dmiss
807this%rect_coo%ilat=dmiss
808this%rect_coo%flon=dmiss
809this%rect_coo%flat=dmiss
811this%box_info%npx=imiss
812this%box_info%npy=imiss
817CALL l4f_category_delete(this%category)
819END SUBROUTINE transform_delete
823SUBROUTINE transform_get_val(this, time_definition, trans_type, sub_type, &
824 input_levtype, output_levtype)
826INTEGER,
INTENT(out),
OPTIONAL :: time_definition
827CHARACTER(len=*),
INTENT(out),
OPTIONAL :: trans_type
828CHARACTER(len=*),
INTENT(out),
OPTIONAL :: sub_type
829TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: input_levtype
831TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: output_levtype
834IF (
PRESENT(time_definition)) time_definition=this%time_definition
835IF (
PRESENT(trans_type)) trans_type = this%trans_type
836IF (
PRESENT(sub_type)) sub_type = this%sub_type
837IF (
PRESENT(input_levtype)) input_levtype = this%vertint%input_levtype
838IF (
PRESENT(output_levtype)) output_levtype = this%vertint%output_levtype
841END SUBROUTINE transform_get_val
887SUBROUTINE grid_transform_levtype_levtype_init(this, trans, lev_in, lev_out, &
888 coord_3d_in, categoryappend)
891TYPE(vol7d_level),
INTENT(in) :: lev_in(:)
892TYPE(vol7d_level),
INTENT(in) :: lev_out(:)
893REAL,
INTENT(inout),
OPTIONAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
894CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
896DOUBLE PRECISION :: coord_in(SIZE(lev_in))
897DOUBLE PRECISION,
ALLOCATABLE :: coord_out(:)
898LOGICAL :: mask_in(SIZE(lev_in))
899LOGICAL,
ALLOCATABLE :: mask_out(:)
901INTEGER :: i, j, icache, inused, istart, iend, ostart, oend
904CALL grid_transform_init_common(this, trans, categoryappend)
906CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vertint")
909IF (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)
1145END SUBROUTINE grid_transform_levtype_levtype_init
1150SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1151TYPE(vol7d_level),
INTENT(in) :: lev(:)
1152LOGICAL,
INTENT(inout) :: mask(:)
1153DOUBLE PRECISION,
INTENT(out) :: coord(:)
1154LOGICAL,
INTENT(out) :: dolog
1157DOUBLE PRECISION :: fact
1164IF (any(lev(k)%level1 == height_level))
THEN
1166ELSE IF (any(lev(k)%level1 == thermo_level))
THEN
1168ELSE IF (any(lev(k)%level1 == sigma_level))
THEN
1174IF (
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
1199mask(:) = mask(:) .AND.
c_e(coord(:))
1201END SUBROUTINE make_vert_coord
1221SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1227REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1228REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1229CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1231INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1232 xnmin, xnmax, ynmin, ynmax
1233DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1234 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1235TYPE(geo_proj) :: proj_in, proj_out
1237LOGICAL,
ALLOCATABLE :: point_mask(:,:)
1241CALL grid_transform_init_common(this, trans, categoryappend)
1243CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-vg6d")
1247CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1248CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1250IF (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)
1363ELSE 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)
1392ELSE 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)
1436ELSE 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)
1467ELSE 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)))
1525ELSE 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
1605ELSE 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
1653ELSE 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)//
']')
1765SUBROUTINE outgrid_setup()
1768CALL griddim_setsteps(out)
1770CALL get_val(in,
proj=proj_in, component_flag=cf_i)
1771CALL get_val(out,
proj=proj_out, component_flag=cf_o)
1772IF (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)
1789CALL griddim_set_central_lon(in, griddim_central_lon(out))
1791END SUBROUTINE outgrid_setup
1793END SUBROUTINE grid_transform_init
1838SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1839 maskgrid, maskbounds, find_index, categoryappend)
1843TYPE(
vol7d),
INTENT(inout) :: v7d_out
1844REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1845REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1846PROCEDURE(basic_find_index),
POINTER,
OPTIONAL :: find_index
1847CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1849INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1851DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1852DOUBLE PRECISION,
ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1853REAL,
ALLOCATABLE :: lmaskbounds(:)
1858IF (
PRESENT(find_index))
THEN
1859 IF (
ASSOCIATED(find_index))
THEN
1860 this%find_index => find_index
1863CALL grid_transform_init_common(this, trans, categoryappend)
1865CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-v7d")
1869CALL get_val(trans, time_definition=time_definition)
1870IF (.NOT.
c_e(time_definition))
THEN
1874IF (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)
1921ELSE 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)))
1989ELSE 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)))
2056ELSE 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))))
2123ELSE 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))
2319SUBROUTINE gen_mask_class()
2320REAL :: startmaskclass, mmin, mmax
2323IF (
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')
2337mmin = minval(maskgrid, mask=
c_e(maskgrid))
2338mmax = maxval(maskgrid, mask=
c_e(maskgrid))
2340nmaskarea = int(mmax - mmin + 1.5)
2341startmaskclass = mmin - 0.5
2343ALLOCATE(lmaskbounds(nmaskarea+1))
2344lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2346CALL l4f_category_log(this%category,l4f_debug, &
2347 'maskinter, '//t2c(nmaskarea)//
' subareas defined, with boundaries:')
2348DO i = 1,
SIZE(lmaskbounds)
2349 CALL l4f_category_log(this%category,l4f_debug, &
2350 'maskinter '//t2c(i)//
' '//t2c(lmaskbounds(i)))
2354END SUBROUTINE gen_mask_class
2356END SUBROUTINE grid_transform_grid_vol7d_init
2377SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2380TYPE(
vol7d),
INTENT(in) :: v7d_in
2382character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2385DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2386DOUBLE PRECISION,
ALLOCATABLE :: lon(:,:),lat(:,:)
2390CALL grid_transform_init_common(this, trans, categoryappend)
2392CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-vg6d")
2395IF (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)
2426ELSE 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)
2467END SUBROUTINE grid_transform_vol7d_grid_init
2504SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2505 maskbounds, categoryappend)
2508TYPE(
vol7d),
INTENT(in) :: v7d_in
2509TYPE(
vol7d),
INTENT(inout) :: v7d_out
2510REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2511CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2514DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2516DOUBLE PRECISION :: lon1, lat1
2520CALL grid_transform_init_common(this, trans, categoryappend)
2522CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-v7d")
2525IF (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))
2546ELSE 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)
2579ELSE 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()
2728SUBROUTINE metamorphosis_all_setup()
2730this%outnx =
SIZE(v7d_in%ana)
2732this%point_index(:,1) = (/(i,i=1,this%innx)/)
2733CALL vol7d_alloc(v7d_out, nana=
SIZE(v7d_in%ana))
2734v7d_out%ana = v7d_in%ana
2738END SUBROUTINE metamorphosis_all_setup
2740END SUBROUTINE grid_transform_vol7d_vol7d_init
2744SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
2747CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2749CHARACTER(len=512) :: a_name
2751IF (
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))
2757this%category=l4f_category_get(a_name)
2760CALL l4f_category_log(this%category,l4f_debug,
"start init_grid_transform")
2765END SUBROUTINE grid_transform_init_common
2769SUBROUTINE poly_to_coordinates(poly, v7d_out)
2771TYPE(
vol7d),
INTENT(inout) :: v7d_out
2774DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2776DO 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
2785END SUBROUTINE poly_to_coordinates
2790SUBROUTINE grid_transform_delete(this)
2808if (
associated(this%inter_index_x))
deallocate (this%inter_index_x)
2809if (
associated(this%inter_index_y))
deallocate (this%inter_index_y)
2810if (
associated(this%inter_index_z))
deallocate (this%inter_index_z)
2811if (
associated(this%point_index))
deallocate (this%point_index)
2813if (
associated(this%inter_x))
deallocate (this%inter_x)
2814if (
associated(this%inter_y))
deallocate (this%inter_y)
2816if (
associated(this%inter_xp))
deallocate (this%inter_xp)
2817if (
associated(this%inter_yp))
deallocate (this%inter_yp)
2818if (
associated(this%inter_zp))
deallocate (this%inter_zp)
2819if (
associated(this%vcoord_in))
deallocate (this%vcoord_in)
2820if (
associated(this%vcoord_out))
deallocate (this%vcoord_out)
2821if (
associated(this%point_mask))
deallocate (this%point_mask)
2822if (
associated(this%stencil))
deallocate (this%stencil)
2823if (
associated(this%output_level_auto))
deallocate (this%output_level_auto)
2824IF (
ALLOCATED(this%coord_3d_in))
DEALLOCATE(this%coord_3d_in)
2828call l4f_category_delete(this%category)
2830END SUBROUTINE grid_transform_delete
2837SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
2838 point_index, output_point_index, levshift, levused)
2840TYPE(vol7d_level),
POINTER,
OPTIONAL :: output_level_auto(:)
2841LOGICAL,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_mask(:)
2842INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_index(:)
2843INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: output_point_index(:)
2844INTEGER,
INTENT(out),
OPTIONAL :: levshift
2845INTEGER,
INTENT(out),
OPTIONAL :: levused
2849IF (
PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
2850IF (
PRESENT(point_mask))
THEN
2851 IF (
ASSOCIATED(this%point_index))
THEN
2852 point_mask =
c_e(reshape(this%point_index, (/
SIZE(this%point_index)/)))
2855IF (
PRESENT(point_index))
THEN
2856 IF (
ASSOCIATED(this%point_index))
THEN
2857 point_index = reshape(this%point_index, (/
SIZE(this%point_index)/))
2860IF (
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)/)
2870IF (
PRESENT(levshift)) levshift = this%levshift
2871IF (
PRESENT(levused)) levused = this%levused
2873END SUBROUTINE grid_transform_get_val
2878FUNCTION grid_transform_c_e(this)
2880LOGICAL :: grid_transform_c_e
2882grid_transform_c_e = this%valid
2884END FUNCTION grid_transform_c_e
2896RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
2899REAL,
INTENT(in) :: field_in(:,:,:)
2900REAL,
INTENT(out) :: field_out(:,:,:)
2901TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
2902REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
2904INTEGER :: i, j, k, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
2905 kfound, kfoundin, inused, i1, i2, j1, j2, np, ns
2906INTEGER,
ALLOCATABLE :: nval(:,:)
2907REAL :: z1,z2,z3,z4,z(4)
2908DOUBLE PRECISION :: x1,x3,y1,y3,xp,yp
2909INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype
2910REAL,
ALLOCATABLE :: coord_in(:)
2911LOGICAL,
ALLOCATABLE :: mask_in(:)
2912REAL,
ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
2913REAL,
POINTER :: coord_3d_in_act(:,:,:)
2915LOGICAL :: alloc_coord_3d_in_act, nm1
2919CALL l4f_category_log(this%category,l4f_debug,
"start grid_transform_compute")
2922field_out(:,:,:) = rmiss
2924IF (.NOT.this%valid)
THEN
2925 CALL l4f_category_log(this%category,l4f_error, &
2926 "refusing to perform a non valid transformation")
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)
2960IF (
PRESENT(var))
THEN
2961 vartype = vol7d_vartype(var)
2964innx =
SIZE(field_in,1); inny =
SIZE(field_in,2); innz =
SIZE(field_in,3)
2965outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
2968IF (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))
3026call l4f_category_log(this%category,l4f_debug, &
3027 "start grid_transform_compute "//trim(this%trans%trans_type)//
':'// &
3028 trim(this%trans%sub_type))
3031IF (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,:)
3038ELSE 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
3186ELSE 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))
3268ELSE 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
3411ELSE 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
3584ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
3587 WHERE(
c_e(this%inter_index_x(:,:)))
3588 field_out(:,:,k) = real(this%inter_index_x(:,:))
3592ELSE 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(:,:,:)
3688ELSE 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)
3942ELSE IF (this%trans%trans_type ==
'' .OR. this%trans%trans_type ==
'none')
THEN
3944 field_out(:,:,:) = field_in(:,:,:)
3954FUNCTION interp_var_360(v1, v2, w1, w2)
3955REAL,
INTENT(in) :: v1, v2, w1, w2
3956REAL :: interp_var_360
3960IF (abs(v1 - v2) > 180.)
THEN
3968 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
3970 interp_var_360 = v1*w2 + v2*w1
3973END FUNCTION interp_var_360
3976RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask)
RESULT(prev)
3977REAL,
INTENT(in) :: v1(:,:)
3978REAL,
INTENT(in) :: l, h, res
3979LOGICAL,
INTENT(in),
OPTIONAL :: mask(:,:)
3987IF ((h - l) <= res)
THEN
3992IF (
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)
4001ELSE IF (nl > nh)
THEN
4002 prev = find_prevailing_direction(v1, l, m, res)
4003ELSE IF (nl == 0 .AND. nh == 0)
THEN
4009END FUNCTION find_prevailing_direction
4012END SUBROUTINE grid_transform_compute
4020SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4023REAL,
INTENT(in) :: field_in(:,:)
4024REAL,
INTENT(out):: field_out(:,:,:)
4025TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
4026REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
4028real,
allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4029INTEGER :: inn_p, ier, k
4030INTEGER :: innx, inny, innz, outnx, outny, outnz
4033call l4f_category_log(this%category,l4f_debug,
"start v7d_grid_transform_compute")
4036field_out(:,:,:) = rmiss
4038IF (.NOT.this%valid)
THEN
4039 call l4f_category_log(this%category,l4f_error, &
4040 "refusing to perform a non valid transformation")
4045innx =
SIZE(field_in,1); inny = 1; innz =
SIZE(field_in,2)
4046outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
4049IF (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))
4112IF (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)")
4164ELSE 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, &
4175END SUBROUTINE grid_transform_v7d_grid_compute
4190ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
RESULT(zp)
4191REAL,
INTENT(in) :: z1,z2,z3,z4
4192DOUBLE PRECISION,
INTENT(in):: x1,y1
4193DOUBLE PRECISION,
INTENT(in):: x3,y3
4194DOUBLE PRECISION,
INTENT(in):: xp,yp
4201p2 = real((yp-y1)/(y3-y1))
4202p1 = real((xp-x1)/(x3-x1))
4213FUNCTION shapiro (z,zp)
RESULT(zs)
4216REAL,
INTENT(in) :: z(:)
4217REAL,
INTENT(in) :: zp
4223 zs = zp - 0.5* ( n*zp - sum(z,
c_e(z)) )/n
4232SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
4233 lon, lat, extrap, index_x, index_y)
4235logical,
INTENT(in) :: near
4236INTEGER,
INTENT(in) :: nx,ny
4237DOUBLE PRECISION,
INTENT(in) :: xmin, xmax, ymin, ymax
4238DOUBLE PRECISION,
INTENT(in) :: lon(:,:),lat(:,:)
4239LOGICAL,
INTENT(in) :: extrap
4240INTEGER,
INTENT(out) :: index_x(:,:),index_y(:,:)
4243DOUBLE 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)
4271END SUBROUTINE basic_find_index
Copy an object, creating a fully new instance.
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...