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))//
'/'// &
567 trim(
to_char(this%rect_coo%flat)))
568 call raise_fatal_error()
573 call l4f_category_log(this%category,l4f_error,
"zoom: coord parameters missing")
574 call raise_fatal_error()
578 else if (this%sub_type ==
'coordbb')
then
580 if (
c_e(this%rect_coo%ilon) .and.
c_e(this%rect_coo%ilat) .and. &
581 c_e(this%rect_coo%flon) .and.
c_e(this%rect_coo%flat))
then
584 call l4f_category_log(this%category,l4f_error,
"zoom: coordbb parameters missing")
585 call raise_fatal_error()
589 else if (this%sub_type ==
'index')
then
591 IF (
c_e(this%rect_ind%ix) .AND.
c_e(this%rect_ind%iy) .AND. &
592 c_e(this%rect_ind%fx) .AND.
c_e(this%rect_ind%fy))
THEN
595 IF (this%rect_ind%ix > this%rect_ind%fx .OR. &
596 this%rect_ind%iy > this%rect_ind%fy)
THEN
598 CALL l4f_category_log(this%category,l4f_error,
'invalid zoom indices: ')
599 CALL l4f_category_log(this%category,l4f_error, &
600 trim(
to_char(this%rect_ind%ix))//
'/'// &
601 trim(
to_char(this%rect_ind%fx)))
602 CALL l4f_category_log(this%category,l4f_error, &
603 trim(
to_char(this%rect_ind%iy))//
'/'// &
604 trim(
to_char(this%rect_ind%fy)))
606 CALL raise_fatal_error()
611 CALL l4f_category_log(this%category,l4f_error,&
612 'zoom: index parameters ix, iy, fx, fy not provided')
613 CALL raise_fatal_error()
618 CALL sub_type_error()
622ELSE IF (this%trans_type ==
'inter' .OR. this%trans_type ==
'intersearch')
THEN
624 IF (this%sub_type ==
'near' .OR. this%sub_type ==
'bilin' .OR. &
625 this%sub_type ==
'linear' .OR. this%sub_type ==
'shapiro_near')
THEN
628 CALL sub_type_error()
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' .OR. this%trans%trans_type ==
'intersearch')
THEN
1394 CALL outgrid_setup()
1396 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin'&
1397 .OR. this%trans%sub_type ==
'shapiro_near')
THEN
1399 CALL get_val(in, nx=this%innx, ny=this%inny, &
1400 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1401 CALL get_val(out, nx=this%outnx, ny=this%outny)
1403 ALLOCATE(this%inter_index_x(this%outnx,this%outny), &
1404 this%inter_index_y(this%outnx,this%outny))
1405 CALL copy(out, lout)
1408 IF (this%trans%sub_type ==
'bilin')
THEN
1409 CALL this%find_index(in, .false., &
1410 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1411 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1412 this%inter_index_x, this%inter_index_y)
1414 ALLOCATE(this%inter_x(this%innx,this%inny), &
1415 this%inter_y(this%innx,this%inny))
1416 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1417 this%inter_yp(this%outnx,this%outny))
1420 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1422 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1425 CALL this%find_index(in, .true., &
1426 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1427 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1428 this%inter_index_x, this%inter_index_y)
1430 IF (this%trans%trans_type ==
'intersearch')
THEN
1431 ALLOCATE(this%inter_x(this%innx,this%inny), &
1432 this%inter_y(this%innx,this%inny))
1433 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1434 this%inter_yp(this%outnx,this%outny))
1437 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1439 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1447ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
1449 CALL outgrid_setup()
1450 CALL get_val(in, nx=this%innx, ny=this%inny)
1451 CALL get_val(out, nx=this%outnx, ny=this%outny, &
1452 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1455 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
1456 CALL get_val(out, dx=this%trans%area_info%boxdx)
1457 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
1458 CALL get_val(out, dx=this%trans%area_info%boxdy)
1460 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1461 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1463 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1464 this%inter_index_y(this%innx,this%inny))
1470 CALL this%find_index(out, .true., &
1471 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1472 lin%dim%lon, lin%dim%lat, .false., &
1473 this%inter_index_x, this%inter_index_y)
1478ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
1480 CALL outgrid_setup()
1482 CALL get_val(in, nx=this%innx, ny=this%inny, &
1483 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1484 CALL get_val(out, nx=this%outnx, ny=this%outny)
1486 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1487 this%inter_index_y(this%outnx,this%outny))
1488 CALL copy(out, lout)
1490 CALL this%find_index(in, .true., &
1491 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1492 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1493 this%inter_index_x, this%inter_index_y)
1496 nr = int(this%trans%area_info%radius)
1499 r2 = this%trans%area_info%radius**2
1500 ALLOCATE(this%stencil(n,n))
1501 this%stencil(:,:) = .true.
1504 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1511 xnmax = this%innx - nr
1513 ynmax = this%inny - nr
1514 DO iy = 1, this%outny
1515 DO ix = 1, this%outnx
1516 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1517 this%inter_index_x(ix,iy) > xnmax .OR. &
1518 this%inter_index_y(ix,iy) < ynmin .OR. &
1519 this%inter_index_y(ix,iy) > ynmax)
THEN
1520 this%inter_index_x(ix,iy) = imiss
1521 this%inter_index_y(ix,iy) = imiss
1527 CALL l4f_category_log(this%category, l4f_debug, &
1528 'stencilinter: stencil size '//t2c(n*n))
1529 CALL l4f_category_log(this%category, l4f_debug, &
1530 'stencilinter: stencil points '//t2c(count(this%stencil)))
1536ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
1538 IF (this%trans%sub_type ==
'poly')
THEN
1541 CALL get_val(in, nx=this%innx, ny=this%inny)
1542 this%outnx = this%innx
1543 this%outny = this%inny
1546 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1547 this%inter_index_y(this%innx,this%inny))
1548 this%inter_index_x(:,:) = imiss
1549 this%inter_index_y(:,:) = 1
1558 inside_x:
DO i = 1, this%innx
1559 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1561 DO n = nprev, this%trans%poly%arraysize
1562 IF (
inside(point, this%trans%poly%array(n)))
THEN
1563 this%inter_index_x(i,j) = n
1568 DO n = nprev-1, 1, -1
1569 IF (
inside(point, this%trans%poly%array(n)))
THEN
1570 this%inter_index_x(i,j) = n
1581 ELSE IF (this%trans%sub_type ==
'grid')
THEN
1584 CALL copy(out, lout)
1587 CALL get_val(in, nx=this%innx, ny=this%inny)
1588 this%outnx = this%innx
1589 this%outny = this%inny
1590 CALL get_val(lout, nx=nx, ny=ny, &
1591 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1594 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1595 this%inter_index_y(this%innx,this%inny))
1601 CALL this%find_index(lout, .true., &
1602 nx, ny, xmin, xmax, ymin, ymax, &
1603 out%dim%lon, out%dim%lat, .false., &
1604 this%inter_index_x, this%inter_index_y)
1606 WHERE(
c_e(this%inter_index_x(:,:)))
1607 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1608 (this%inter_index_y(:,:)-1)*nx
1616ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
1622 CALL get_val(in, nx=this%innx, ny=this%inny)
1623 this%outnx = this%innx
1624 this%outny = this%inny
1627 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1628 this%inter_index_y(this%innx,this%inny))
1629 this%inter_index_x(:,:) = imiss
1630 this%inter_index_y(:,:) = 1
1639 inside_x_2:
DO i = 1, this%innx
1640 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1642 DO n = nprev, this%trans%poly%arraysize
1643 IF (
inside(point, this%trans%poly%array(n)))
THEN
1644 this%inter_index_x(i,j) = n
1649 DO n = nprev-1, 1, -1
1650 IF (
inside(point, this%trans%poly%array(n)))
THEN
1651 this%inter_index_x(i,j) = n
1664ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
1667 CALL get_val(in, nx=this%innx, ny=this%inny)
1668 this%outnx = this%innx
1669 this%outny = this%inny
1671 IF (this%trans%sub_type ==
'maskvalid' .OR. this%trans%sub_type ==
'maskinvalid')
THEN
1673 IF (.NOT.
PRESENT(maskgrid))
THEN
1674 CALL l4f_category_log(this%category,l4f_error, &
1675 'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1676 trim(this%trans%sub_type)//
' transformation')
1681 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
1682 CALL l4f_category_log(this%category,l4f_error, &
1683 'grid_transform_init mask non conformal with input field')
1684 CALL l4f_category_log(this%category,l4f_error, &
1685 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
1686 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
1691 ALLOCATE(this%point_mask(this%innx,this%inny))
1693 IF (this%trans%sub_type ==
'maskvalid')
THEN
1696 IF (.NOT.
PRESENT(maskbounds))
THEN
1697 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1698 ELSE IF (
SIZE(maskbounds) < 2)
THEN
1699 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1701 this%point_mask(:,:) =
c_e(maskgrid(:,:)) .AND. &
1702 maskgrid(:,:) > maskbounds(1) .AND. &
1703 maskgrid(:,:) <= maskbounds(
SIZE(maskbounds))
1706 this%point_mask(:,:) = .NOT.
c_e(maskgrid(:,:))
1711 ELSE IF (this%trans%sub_type ==
'lemaskinvalid' .OR. &
1712 this%trans%sub_type ==
'ltmaskinvalid' .OR. &
1713 this%trans%sub_type ==
'gemaskinvalid' .OR. &
1714 this%trans%sub_type ==
'gtmaskinvalid')
THEN
1717 this%val_mask = maskgrid
1720 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
1722 IF (.NOT.
PRESENT(maskbounds))
THEN
1723 CALL l4f_category_log(this%category,l4f_error, &
1724 'grid_transform_init maskbounds missing for metamorphosis:'// &
1725 trim(this%trans%sub_type)//
' transformation')
1727 ELSE IF (
SIZE(maskbounds) < 1)
THEN
1728 CALL l4f_category_log(this%category,l4f_error, &
1729 'grid_transform_init maskbounds empty for metamorphosis:'// &
1730 trim(this%trans%sub_type)//
' transformation')
1733 this%val1 = maskbounds(1)
1735 CALL l4f_category_log(this%category, l4f_debug, &
1736 "grid_transform_init setting invalid data to "//t2c(this%val1))
1742 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
1744 IF (.NOT.
PRESENT(maskbounds))
THEN
1745 CALL l4f_category_log(this%category,l4f_error, &
1746 'grid_transform_init maskbounds missing for metamorphosis:'// &
1747 trim(this%trans%sub_type)//
' transformation')
1750 ELSE IF (
SIZE(maskbounds) < 2)
THEN
1751 CALL l4f_category_log(this%category,l4f_error, &
1752 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1753 trim(this%trans%sub_type)//
' transformation')
1757 this%val1 = maskbounds(1)
1758 this%val2 = maskbounds(
SIZE(maskbounds))
1760 CALL l4f_category_log(this%category, l4f_debug, &
1761 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
1762 t2c(this%val2)//
']')
1776SUBROUTINE outgrid_setup()
1779CALL griddim_setsteps(out)
1781CALL get_val(in,
proj=proj_in, component_flag=cf_i)
1782CALL get_val(out,
proj=proj_out, component_flag=cf_o)
1783IF (proj_in == proj_out)
THEN
1786 CALL set_val(out, component_flag=cf_i)
1791 CALL l4f_category_log(this%category,l4f_warn, &
1792 'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1793 CALL l4f_category_log(this%category,l4f_warn, &
1794 'vector fields will probably be wrong')
1796 CALL set_val(out, component_flag=cf_i)
1800CALL griddim_set_central_lon(in, griddim_central_lon(out))
1802END SUBROUTINE outgrid_setup
1804END SUBROUTINE grid_transform_init
1849SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1850 maskgrid, maskbounds, find_index, categoryappend)
1854TYPE(
vol7d),
INTENT(inout) :: v7d_out
1855REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1856REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1857PROCEDURE(basic_find_index),
POINTER,
OPTIONAL :: find_index
1858CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1860INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1862DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1863DOUBLE PRECISION,
ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1864REAL,
ALLOCATABLE :: lmaskbounds(:)
1869IF (
PRESENT(find_index))
THEN
1870 IF (
ASSOCIATED(find_index))
THEN
1871 this%find_index => find_index
1874CALL grid_transform_init_common(this, trans, categoryappend)
1876CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-v7d")
1880CALL get_val(trans, time_definition=time_definition)
1881IF (.NOT.
c_e(time_definition))
THEN
1885IF (this%trans%trans_type ==
'inter')
THEN
1887 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin' &
1888 .OR. this%trans%sub_type ==
'shapiro_near')
THEN
1891 CALL get_val(lin, nx=this%innx, ny=this%inny)
1892 this%outnx =
SIZE(v7d_out%ana)
1895 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1896 this%inter_index_y(this%outnx,this%outny))
1897 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1899 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1901 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1902 CALL griddim_set_central_lon(lin, lonref)
1903 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1905 IF (this%trans%sub_type ==
'bilin')
THEN
1906 CALL this%find_index(lin, .false., &
1907 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1908 lon, lat, this%trans%extrap, &
1909 this%inter_index_x, this%inter_index_y)
1911 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1912 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1914 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1915 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1918 CALL this%find_index(lin, .true., &
1919 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1920 lon, lat, this%trans%extrap, &
1921 this%inter_index_x, this%inter_index_y)
1923 IF (this%trans%trans_type ==
'intersearch')
THEN
1924 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1925 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1927 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1928 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1939ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
1941 CALL get_val(in, nx=this%innx, ny=this%inny)
1943 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1944 this%inter_index_y(this%innx,this%inny))
1945 this%inter_index_x(:,:) = imiss
1946 this%inter_index_y(:,:) = 1
1952 this%outnx = this%trans%poly%arraysize
1955 CALL init(v7d_out, time_definition=time_definition)
1956 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1959 ALLOCATE(lon(this%outnx,1))
1960 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
1961 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1962 CALL griddim_set_central_lon(lin, lonref)
1967 CALL poly_to_coordinates(this%trans%poly, v7d_out)
1972 DO iy = 1, this%inny
1973 inside_x:
DO ix = 1, this%innx
1974 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
1976 DO n = nprev, this%trans%poly%arraysize
1977 IF (
inside(point, this%trans%poly%array(n)))
THEN
1978 this%inter_index_x(ix,iy) = n
1983 DO n = nprev-1, 1, -1
1984 IF (
inside(point, this%trans%poly%array(n)))
THEN
1985 this%inter_index_x(ix,iy) = n
1997 DO n = 1, this%trans%poly%arraysize
1998 CALL l4f_category_log(this%category, l4f_debug, &
1999 'Polygon: '//t2c(n)//
' grid points: '// &
2000 t2c(count(this%inter_index_x(:,:) == n)))
2007ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
2011 CALL get_val(lin, nx=this%innx, ny=this%inny)
2012 this%outnx =
SIZE(v7d_out%ana)
2015 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
2016 this%inter_index_y(this%outnx,this%outny))
2017 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
2019 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2021 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2022 CALL griddim_set_central_lon(lin, lonref)
2024 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2026 CALL this%find_index(lin, .true., &
2027 this%innx, this%inny, xmin, xmax, ymin, ymax, &
2028 lon, lat, this%trans%extrap, &
2029 this%inter_index_x, this%inter_index_y)
2032 nr = int(this%trans%area_info%radius)
2035 r2 = this%trans%area_info%radius**2
2036 ALLOCATE(this%stencil(n,n))
2037 this%stencil(:,:) = .true.
2040 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
2047 xnmax = this%innx - nr
2049 ynmax = this%inny - nr
2050 DO iy = 1, this%outny
2051 DO ix = 1, this%outnx
2052 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2053 this%inter_index_x(ix,iy) > xnmax .OR. &
2054 this%inter_index_y(ix,iy) < ynmin .OR. &
2055 this%inter_index_y(ix,iy) > ynmax)
THEN
2056 this%inter_index_x(ix,iy) = imiss
2057 this%inter_index_y(ix,iy) = imiss
2063 CALL l4f_category_log(this%category, l4f_debug, &
2064 'stencilinter: stencil size '//t2c(n*n))
2065 CALL l4f_category_log(this%category, l4f_debug, &
2066 'stencilinter: stencil points '//t2c(count(this%stencil)))
2074ELSE IF (this%trans%trans_type ==
'maskinter')
THEN
2076 IF (.NOT.
PRESENT(maskgrid))
THEN
2077 CALL l4f_category_log(this%category,l4f_error, &
2078 'grid_transform_init maskgrid argument missing for maskinter transformation')
2079 CALL raise_fatal_error()
2082 CALL get_val(in, nx=this%innx, ny=this%inny)
2083 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2084 CALL l4f_category_log(this%category,l4f_error, &
2085 'grid_transform_init mask non conformal with input field')
2086 CALL l4f_category_log(this%category,l4f_error, &
2087 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2088 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2089 CALL raise_fatal_error()
2092 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2093 this%inter_index_y(this%innx,this%inny))
2094 this%inter_index_x(:,:) = imiss
2095 this%inter_index_y(:,:) = 1
2098 CALL gen_mask_class()
2106 DO iy = 1, this%inny
2107 DO ix = 1, this%innx
2108 IF (
c_e(maskgrid(ix,iy)))
THEN
2109 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2110 DO n = nmaskarea, 1, -1
2111 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2112 this%inter_index_x(ix,iy) = n
2122 this%outnx = nmaskarea
2125 CALL init(v7d_out, time_definition=time_definition)
2126 CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2131 CALL init(v7d_out%ana(n), &
2133 mask=(this%inter_index_x(:,:) == n))), &
2135 mask=(this%inter_index_x(:,:) == n))))
2141ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
2148 CALL get_val(in, nx=this%innx, ny=this%inny)
2150 ALLOCATE(this%point_index(this%innx,this%inny))
2151 this%point_index(:,:) = imiss
2154 CALL init(v7d_out, time_definition=time_definition)
2156 IF (this%trans%sub_type ==
'all' )
THEN
2158 this%outnx = this%innx*this%inny
2160 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2165 CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2166 lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2168 this%point_index(ix,iy) = n
2174 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN
2179 DO iy = 1, this%inny
2180 DO ix = 1, this%innx
2182 IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2183 lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2184 lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2185 lin%dim%lat(ix,iy) < this%trans%rect_coo%flat)
THEN
2186 this%outnx = this%outnx + 1
2187 this%point_index(ix,iy) = this%outnx
2192 IF (this%outnx <= 0)
THEN
2193 CALL l4f_category_log(this%category,l4f_warn, &
2194 "metamorphosis:coordbb: no points inside bounding box "//&
2195 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2196 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2197 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2198 trim(
to_char(this%trans%rect_coo%flat)))
2201 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2205 DO iy = 1, this%inny
2206 DO ix = 1, this%innx
2207 IF (
c_e(this%point_index(ix,iy)))
THEN
2209 CALL init(v7d_out%ana(n), &
2210 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2217 ELSE IF (this%trans%sub_type ==
'poly' )
THEN
2226 DO iy = 1, this%inny
2227 DO ix = 1, this%innx
2228 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2229 DO n = 1, this%trans%poly%arraysize
2230 IF (
inside(point, this%trans%poly%array(n)))
THEN
2231 this%outnx = this%outnx + 1
2232 this%point_index(ix,iy) = n
2241 IF (this%outnx <= 0)
THEN
2242 CALL l4f_category_log(this%category,l4f_warn, &
2243 "metamorphosis:poly: no points inside polygons")
2246 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2249 DO iy = 1, this%inny
2250 DO ix = 1, this%innx
2251 IF (
c_e(this%point_index(ix,iy)))
THEN
2253 CALL init(v7d_out%ana(n), &
2254 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2261 ELSE IF (this%trans%sub_type ==
'mask' )
THEN
2263 IF (.NOT.
PRESENT(maskgrid))
THEN
2264 CALL l4f_category_log(this%category,l4f_error, &
2265 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2270 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2271 CALL l4f_category_log(this%category,l4f_error, &
2272 'grid_transform_init mask non conformal with input field')
2273 CALL l4f_category_log(this%category,l4f_error, &
2274 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2275 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2281 CALL gen_mask_class()
2289 DO iy = 1, this%inny
2290 DO ix = 1, this%innx
2291 IF (
c_e(maskgrid(ix,iy)))
THEN
2292 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2293 DO n = nmaskarea, 1, -1
2294 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2295 this%outnx = this%outnx + 1
2296 this%point_index(ix,iy) = n
2306 IF (this%outnx <= 0)
THEN
2307 CALL l4f_category_log(this%category,l4f_warn, &
2308 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2312 CALL l4f_category_log(this%category,l4f_info, &
2313 "points in subarea "//t2c(n)//
": "// &
2314 t2c(count(this%point_index(:,:) == n)))
2317 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2320 DO iy = 1, this%inny
2321 DO ix = 1, this%innx
2322 IF (
c_e(this%point_index(ix,iy)))
THEN
2324 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2337SUBROUTINE gen_mask_class()
2338REAL :: startmaskclass, mmin, mmax
2341IF (
PRESENT(maskbounds))
THEN
2342 nmaskarea =
SIZE(maskbounds) - 1
2343 IF (nmaskarea > 0)
THEN
2344 lmaskbounds = maskbounds
2346 ELSE IF (nmaskarea == 0)
THEN
2347 CALL l4f_category_log(this%category,l4f_warn, &
2348 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2349 //
', it will be ignored')
2350 CALL l4f_category_log(this%category,l4f_warn, &
2351 'at least 2 values are required for maskbounds')
2355mmin = minval(maskgrid, mask=
c_e(maskgrid))
2356mmax = maxval(maskgrid, mask=
c_e(maskgrid))
2358nmaskarea = int(mmax - mmin + 1.5)
2359startmaskclass = mmin - 0.5
2361ALLOCATE(lmaskbounds(nmaskarea+1))
2362lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2364CALL l4f_category_log(this%category,l4f_debug, &
2365 'maskinter, '//t2c(nmaskarea)//
' subareas defined, with boundaries:')
2366DO i = 1,
SIZE(lmaskbounds)
2367 CALL l4f_category_log(this%category,l4f_debug, &
2368 'maskinter '//t2c(i)//
' '//t2c(lmaskbounds(i)))
2372END SUBROUTINE gen_mask_class
2374END SUBROUTINE grid_transform_grid_vol7d_init
2395SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2398TYPE(
vol7d),
INTENT(in) :: v7d_in
2400character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2403DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2404DOUBLE PRECISION,
ALLOCATABLE :: lon(:,:),lat(:,:)
2408CALL grid_transform_init_common(this, trans, categoryappend)
2410CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-vg6d")
2413IF (this%trans%trans_type ==
'inter')
THEN
2415 IF ( this%trans%sub_type ==
'linear' )
THEN
2417 this%innx=
SIZE(v7d_in%ana)
2419 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2420 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2421 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2423 CALL copy (out, lout)
2425 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2426 CALL griddim_set_central_lon(lout, lonref)
2428 CALL get_val(lout, nx=nx, ny=ny)
2431 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2433 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2434 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2435 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2444ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
2446 this%innx=
SIZE(v7d_in%ana)
2449 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2450 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2451 this%inter_index_y(this%innx,this%inny))
2453 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2455 CALL copy (out, lout)
2457 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2458 CALL griddim_set_central_lon(lout, lonref)
2460 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2461 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2464 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
2465 CALL get_val(out, dx=this%trans%area_info%boxdx)
2466 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
2467 CALL get_val(out, dx=this%trans%area_info%boxdy)
2469 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2470 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2473 CALL this%find_index(lout, .true., &
2474 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2475 lon, lat, .false., &
2476 this%inter_index_x, this%inter_index_y)
2485END SUBROUTINE grid_transform_vol7d_grid_init
2522SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2523 maskbounds, categoryappend)
2526TYPE(
vol7d),
INTENT(in) :: v7d_in
2527TYPE(
vol7d),
INTENT(inout) :: v7d_out
2528REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2529CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2532DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2534DOUBLE PRECISION :: lon1, lat1
2538CALL grid_transform_init_common(this, trans, categoryappend)
2540CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-v7d")
2543IF (this%trans%trans_type ==
'inter')
THEN
2545 IF ( this%trans%sub_type ==
'linear' )
THEN
2547 CALL vol7d_alloc_vol(v7d_out)
2548 this%outnx=
SIZE(v7d_out%ana)
2551 this%innx=
SIZE(v7d_in%ana)
2554 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2555 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2557 CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2558 CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2564ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
2566 this%innx=
SIZE(v7d_in%ana)
2569 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2570 this%inter_index_y(this%innx,this%inny))
2571 this%inter_index_x(:,:) = imiss
2572 this%inter_index_y(:,:) = 1
2574 DO i = 1,
SIZE(v7d_in%ana)
2576 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2577 point = georef_coord_new(x=lon1, y=lat1)
2579 DO n = 1, this%trans%poly%arraysize
2580 IF (
inside(point, this%trans%poly%array(n)))
THEN
2581 this%inter_index_x(i,1) = n
2587 this%outnx=this%trans%poly%arraysize
2589 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2593 CALL poly_to_coordinates(this%trans%poly, v7d_out)
2597ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
2600 this%innx =
SIZE(v7d_in%ana)
2603 ALLOCATE(this%point_index(this%innx,this%inny))
2604 this%point_index(:,:) = imiss
2606 IF (this%trans%sub_type ==
'all' )
THEN
2608 CALL metamorphosis_all_setup()
2610 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN
2612 ALLOCATE(lon(this%innx),lat(this%innx))
2617 CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2620 IF (lon(i) > this%trans%rect_coo%ilon .AND. &
2621 lon(i) < this%trans%rect_coo%flon .AND. &
2622 lat(i) > this%trans%rect_coo%ilat .AND. &
2623 lat(i) < this%trans%rect_coo%flat)
THEN
2624 this%outnx = this%outnx + 1
2625 this%point_index(i,1) = this%outnx
2629 IF (this%outnx <= 0)
THEN
2630 CALL l4f_category_log(this%category,l4f_warn, &
2631 "metamorphosis:coordbb: no points inside bounding box "//&
2632 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2633 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2634 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2635 trim(
to_char(this%trans%rect_coo%flat)))
2638 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2643 IF (
c_e(this%point_index(i,1)))
THEN
2645 CALL init(v7d_out%ana(n),lon=lon(i),lat=lat(i))
2648 DEALLOCATE(lon, lat)
2652 ELSE IF (this%trans%sub_type ==
'poly' )
THEN
2659 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2660 point = georef_coord_new(x=lon1, y=lat1)
2661 DO n = 1, this%trans%poly%arraysize
2662 IF (
inside(point, this%trans%poly%array(n)))
THEN
2663 this%outnx = this%outnx + 1
2664 this%point_index(i,1) = n
2671 IF (this%outnx <= 0)
THEN
2672 CALL l4f_category_log(this%category,l4f_warn, &
2673 "metamorphosis:poly: no points inside polygons")
2676 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2681 IF (
c_e(this%point_index(i,1)))
THEN
2684 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2685 CALL init(v7d_out%ana(n),lon=lon1,lat=lat1)
2691 ELSE IF (this%trans%sub_type ==
'setinvalidto' )
THEN
2693 IF (.NOT.
PRESENT(maskbounds))
THEN
2694 CALL l4f_category_log(this%category,l4f_error, &
2695 'grid_transform_init maskbounds missing for metamorphosis:'// &
2696 trim(this%trans%sub_type)//
' transformation')
2698 ELSE IF (
SIZE(maskbounds) < 1)
THEN
2699 CALL l4f_category_log(this%category,l4f_error, &
2700 'grid_transform_init maskbounds empty for metamorphosis:'// &
2701 trim(this%trans%sub_type)//
' transformation')
2704 this%val1 = maskbounds(1)
2706 CALL l4f_category_log(this%category, l4f_debug, &
2707 "grid_transform_init setting invalid data to "//t2c(this%val1))
2711 CALL metamorphosis_all_setup()
2713 ELSE IF (this%trans%sub_type ==
'settoinvalid' )
THEN
2715 IF (.NOT.
PRESENT(maskbounds))
THEN
2716 CALL l4f_category_log(this%category,l4f_error, &
2717 'grid_transform_init maskbounds missing for metamorphosis:'// &
2718 trim(this%trans%sub_type)//
' transformation')
2721 ELSE IF (
SIZE(maskbounds) < 2)
THEN
2722 CALL l4f_category_log(this%category,l4f_error, &
2723 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
2724 trim(this%trans%sub_type)//
' transformation')
2728 this%val1 = maskbounds(1)
2729 this%val2 = maskbounds(
SIZE(maskbounds))
2731 CALL l4f_category_log(this%category, l4f_debug, &
2732 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
2733 t2c(this%val2)//
']')
2737 CALL metamorphosis_all_setup()
2746SUBROUTINE metamorphosis_all_setup()
2748this%outnx =
SIZE(v7d_in%ana)
2750this%point_index(:,1) = (/(i,i=1,this%innx)/)
2751CALL vol7d_alloc(v7d_out, nana=
SIZE(v7d_in%ana))
2752v7d_out%ana = v7d_in%ana
2756END SUBROUTINE metamorphosis_all_setup
2758END SUBROUTINE grid_transform_vol7d_vol7d_init
2762SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
2765CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2767CHARACTER(len=512) :: a_name
2769IF (
PRESENT(categoryappend))
THEN
2770 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
2771 trim(categoryappend))
2773 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2775this%category=l4f_category_get(a_name)
2778CALL l4f_category_log(this%category,l4f_debug,
"start init_grid_transform")
2783END SUBROUTINE grid_transform_init_common
2787SUBROUTINE poly_to_coordinates(poly, v7d_out)
2789TYPE(
vol7d),
INTENT(inout) :: v7d_out
2792DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2794DO n = 1, poly%arraysize
2795 CALL getval(poly%array(n), x=lon, y=lat)
2796 sz = min(
SIZE(lon),
SIZE(lat))
2797 IF (lon(1) == lon(sz) .AND. lat(1) == lat(sz))
THEN
2803END SUBROUTINE poly_to_coordinates
2808SUBROUTINE grid_transform_delete(this)
2826if (
associated(this%inter_index_x))
deallocate (this%inter_index_x)
2827if (
associated(this%inter_index_y))
deallocate (this%inter_index_y)
2828if (
associated(this%inter_index_z))
deallocate (this%inter_index_z)
2829if (
associated(this%point_index))
deallocate (this%point_index)
2831if (
associated(this%inter_x))
deallocate (this%inter_x)
2832if (
associated(this%inter_y))
deallocate (this%inter_y)
2834if (
associated(this%inter_xp))
deallocate (this%inter_xp)
2835if (
associated(this%inter_yp))
deallocate (this%inter_yp)
2836if (
associated(this%inter_zp))
deallocate (this%inter_zp)
2837if (
associated(this%vcoord_in))
deallocate (this%vcoord_in)
2838if (
associated(this%vcoord_out))
deallocate (this%vcoord_out)
2839if (
associated(this%point_mask))
deallocate (this%point_mask)
2840if (
associated(this%stencil))
deallocate (this%stencil)
2841if (
associated(this%output_level_auto))
deallocate (this%output_level_auto)
2842IF (
ALLOCATED(this%coord_3d_in))
DEALLOCATE(this%coord_3d_in)
2846call l4f_category_delete(this%category)
2848END SUBROUTINE grid_transform_delete
2855SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
2856 point_index, output_point_index, levshift, levused)
2858TYPE(vol7d_level),
POINTER,
OPTIONAL :: output_level_auto(:)
2859LOGICAL,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_mask(:)
2860INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_index(:)
2861INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: output_point_index(:)
2862INTEGER,
INTENT(out),
OPTIONAL :: levshift
2863INTEGER,
INTENT(out),
OPTIONAL :: levused
2867IF (
PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
2868IF (
PRESENT(point_mask))
THEN
2869 IF (
ASSOCIATED(this%point_index))
THEN
2870 point_mask =
c_e(reshape(this%point_index, (/
SIZE(this%point_index)/)))
2873IF (
PRESENT(point_index))
THEN
2874 IF (
ASSOCIATED(this%point_index))
THEN
2875 point_index = reshape(this%point_index, (/
SIZE(this%point_index)/))
2878IF (
PRESENT(output_point_index))
THEN
2879 IF (
ASSOCIATED(this%point_index))
THEN
2881 output_point_index = pack(this%point_index(:,:),
c_e(this%point_index))
2882 ELSE IF (this%trans%trans_type ==
'polyinter' .OR. &
2883 this%trans%trans_type ==
'maskinter')
THEN
2885 output_point_index = (/(i,i=1,this%outnx)/)
2888IF (
PRESENT(levshift)) levshift = this%levshift
2889IF (
PRESENT(levused)) levused = this%levused
2891END SUBROUTINE grid_transform_get_val
2896FUNCTION grid_transform_c_e(this)
2898LOGICAL :: grid_transform_c_e
2900grid_transform_c_e = this%valid
2902END FUNCTION grid_transform_c_e
2914RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
2917REAL,
INTENT(in) :: field_in(:,:,:)
2918REAL,
INTENT(out) :: field_out(:,:,:)
2919TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
2920REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
2922INTEGER :: i, j, k, l, m, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
2923 kfound, kfoundin, inused, i1, i2, j1, j2, np, ns
2924INTEGER,
ALLOCATABLE :: nval(:,:)
2925REAL :: z1,z2,z3,z4,z(4)
2926DOUBLE PRECISION :: x1,x3,y1,y3,xp,yp, disttmp, dist
2927INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype
2928REAL,
ALLOCATABLE :: coord_in(:)
2929LOGICAL,
ALLOCATABLE :: mask_in(:)
2930REAL,
ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
2931REAL,
POINTER :: coord_3d_in_act(:,:,:)
2933LOGICAL :: alloc_coord_3d_in_act, nm1
2937CALL l4f_category_log(this%category,l4f_debug,
"start grid_transform_compute")
2940field_out(:,:,:) = rmiss
2942IF (.NOT.this%valid)
THEN
2943 CALL l4f_category_log(this%category,l4f_error, &
2944 "refusing to perform a non valid transformation")
2950 CALL l4f_category_log(this%category,l4f_debug, &
2951 "recursive transformation in grid_transform_compute")
2954 IF (this%trans%trans_type ==
'polyinter')
THEN
2956 likethis%recur = .false.
2957 likethis%outnx = this%trans%poly%arraysize
2959 ALLOCATE(field_tmp(this%trans%poly%arraysize,1,
SIZE(field_out,3)))
2960 CALL grid_transform_compute(likethis, field_in, field_tmp, var, coord_3d_in)
2962 DO k = 1,
SIZE(field_out,3)
2965 IF (
c_e(this%inter_index_x(i,j)))
THEN
2966 field_out(i,j,k) = field_tmp(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2971 DEALLOCATE(field_tmp)
2978IF (
PRESENT(var))
THEN
2979 vartype = vol7d_vartype(var)
2982innx =
SIZE(field_in,1); inny =
SIZE(field_in,2); innz =
SIZE(field_in,3)
2983outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
2986IF (this%trans%trans_type ==
'vertint')
THEN
2988 IF (innz /= this%innz)
THEN
2989 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2990 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
2991 t2c(this%innz)//
" /= "//t2c(innz))
2996 IF (outnz /= this%outnz)
THEN
2997 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2998 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
2999 t2c(this%outnz)//
" /= "//t2c(outnz))
3004 IF (innx /= outnx .OR. inny /= outny)
THEN
3005 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
3006 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
3007 t2c(innx)//
","//t2c(inny)//
" /= "//&
3008 t2c(outnx)//
","//t2c(outny))
3015 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
3016 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3017 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
3018 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
3019 t2c(innx)//
","//t2c(inny))
3024 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
3025 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3026 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
3027 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
3028 t2c(outnx)//
","//t2c(outny))
3033 IF (innz /= outnz)
THEN
3034 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3035 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
3036 t2c(innz)//
" /= "//t2c(outnz))
3044call l4f_category_log(this%category,l4f_debug, &
3045 "start grid_transform_compute "//trim(this%trans%trans_type)//
':'// &
3046 trim(this%trans%sub_type))
3049IF (this%trans%trans_type ==
'zoom')
THEN
3051 field_out(this%outinx:this%outfnx, &
3052 this%outiny:this%outfny,:) = &
3053 field_in(this%iniox:this%infox, &
3054 this%inioy:this%infoy,:)
3056ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN
3058 IF (this%trans%sub_type ==
'average')
THEN
3059 IF (vartype == var_dir360)
THEN
3062 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3063 je = j+this%trans%box_info%npy-1
3066 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3067 ie = i+this%trans%box_info%npx-1
3069 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3071 field_out(ii,jj,k) = find_prevailing_direction(field_in(i:ie,j:je,k), &
3081 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3082 je = j+this%trans%box_info%npy-1
3085 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3086 ie = i+this%trans%box_info%npx-1
3088 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3090 field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), &
3091 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3099 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3100 this%trans%sub_type ==
'stddevnm1')
THEN
3102 IF (this%trans%sub_type ==
'stddev')
THEN
3108 navg = this%trans%box_info%npx*this%trans%box_info%npy
3111 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3112 je = j+this%trans%box_info%npy-1
3115 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3116 ie = i+this%trans%box_info%npx-1
3119 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3124 ELSE IF (this%trans%sub_type ==
'max')
THEN
3127 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3128 je = j+this%trans%box_info%npy-1
3131 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3132 ie = i+this%trans%box_info%npx-1
3134 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3136 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3137 mask=(field_in(i:ie,j:je,k) /= rmiss))
3143 ELSE IF (this%trans%sub_type ==
'min')
THEN
3146 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3147 je = j+this%trans%box_info%npy-1
3150 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3151 ie = i+this%trans%box_info%npx-1
3153 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3155 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3156 mask=(field_in(i:ie,j:je,k) /= rmiss))
3162 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3164 navg = this%trans%box_info%npx*this%trans%box_info%npy
3167 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3168 je = j+this%trans%box_info%npy-1
3171 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3172 ie = i+this%trans%box_info%npx-1
3175 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3176 (/real(this%trans%stat_info%percentile)/))
3181 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3185 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3186 je = j+this%trans%box_info%npy-1
3189 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3190 ie = i+this%trans%box_info%npx-1
3192 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3194 field_out(ii,jj,k) = sum(interval_info_valid( &
3195 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3196 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3204ELSE IF (this%trans%trans_type ==
'inter')
THEN
3206 IF (this%trans%sub_type ==
'near')
THEN
3209 DO j = 1, this%outny
3210 DO i = 1, this%outnx
3212 IF (
c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3213 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3219 ELSE IF (this%trans%sub_type ==
'bilin')
THEN
3222 DO j = 1, this%outny
3223 DO i = 1, this%outnx
3225 IF (
c_e(this%inter_index_x(i,j)))
THEN
3227 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3228 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3229 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3230 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3232 IF (
c_e(z1) .AND.
c_e(z2) .AND.
c_e(z3) .AND.
c_e(z4))
THEN
3234 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3235 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3236 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3237 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3239 xp=this%inter_xp(i,j)
3240 yp=this%inter_yp(i,j)
3242 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3250 ELSE IF (this%trans%sub_type ==
'shapiro_near')
THEN
3252 DO j = 1, this%outny
3253 DO i = 1, this%outnx
3255 IF (
c_e(this%inter_index_x(i,j)))
THEN
3257 IF(this%inter_index_x(i,j)-1>0)
THEN
3258 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3262 IF(this%inter_index_x(i,j)+1<=this%outnx)
THEN
3263 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3267 IF(this%inter_index_y(i,j)+1<=this%outny)
THEN
3268 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3272 IF(this%inter_index_y(i,j)-1>0)
THEN
3273 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3277 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3286ELSE IF (this%trans%trans_type ==
'intersearch')
THEN
3289 likethis%trans%trans_type =
'inter'
3290 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3293 IF ((.NOT.all(
c_e(field_out(:,:,k)))) .AND. (any(
c_e(field_in(:,:,k)))))
THEN
3294 DO j = 1, this%outny
3295 DO i = 1, this%outnx
3296 IF (.NOT.
c_e(field_out(i,j,k)))
THEN
3300 IF (
c_e(field_in(l,m,k)))
THEN
3301 disttmp = (this%inter_xp(l,m) - this%inter_x(i,j))**2 + (this%inter_yp(l,m) - this%inter_y(i,j))**2
3302 IF (disttmp <
dist)
THEN
3304 field_out(i,j,k) = field_in(l,m,k)
3315ELSE IF (this%trans%trans_type ==
'boxinter' &
3316 .OR. this%trans%trans_type ==
'polyinter' &
3317 .OR. this%trans%trans_type ==
'maskinter')
THEN
3319 IF (this%trans%sub_type ==
'average')
THEN
3321 IF (vartype == var_dir360)
THEN
3323 DO j = 1, this%outny
3324 DO i = 1, this%outnx
3325 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3327 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3333 ALLOCATE(nval(this%outnx, this%outny))
3334 field_out(:,:,:) = 0.0
3339 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3340 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3341 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3343 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3344 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3348 WHERE (nval(:,:) /= 0)
3349 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3351 field_out(:,:,k) = rmiss
3357 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3358 this%trans%sub_type ==
'stddevnm1')
THEN
3360 IF (this%trans%sub_type ==
'stddev')
THEN
3366 DO j = 1, this%outny
3367 DO i = 1, this%outnx
3370 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3371 mask=reshape((this%inter_index_x == i .AND. &
3372 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)), nm1=nm1)
3377 ELSE IF (this%trans%sub_type ==
'max')
THEN
3382 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3383 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3384 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3385 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3388 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3397 ELSE IF (this%trans%sub_type ==
'min')
THEN
3402 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3403 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3404 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3405 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3408 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3416 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3419 DO j = 1, this%outny
3420 DO i = 1, this%outnx
3423 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3424 (/real(this%trans%stat_info%percentile)/), &
3425 mask=reshape((this%inter_index_x == i .AND. &
3426 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)))
3431 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3433 ALLOCATE(nval(this%outnx, this%outny))
3434 field_out(:,:,:) = 0.0
3439 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3440 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3441 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3442 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3443 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3444 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3448 WHERE (nval(:,:) /= 0)
3449 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3451 field_out(:,:,k) = rmiss
3458ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
3459 np =
SIZE(this%stencil,1)/2
3460 ns =
SIZE(this%stencil)
3462 IF (this%trans%sub_type ==
'average')
THEN
3464 IF (vartype == var_dir360)
THEN
3466 DO j = 1, this%outny
3467 DO i = 1, this%outnx
3468 IF (
c_e(this%inter_index_x(i,j)))
THEN
3469 i1 = this%inter_index_x(i,j) - np
3470 i2 = this%inter_index_x(i,j) + np
3471 j1 = this%inter_index_y(i,j) - np
3472 j2 = this%inter_index_y(i,j) + np
3473 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3475 mask=this%stencil(:,:))
3486 DO j = 1, this%outny
3487 DO i = 1, this%outnx
3488 IF (
c_e(this%inter_index_x(i,j)))
THEN
3489 i1 = this%inter_index_x(i,j) - np
3490 i2 = this%inter_index_x(i,j) + np
3491 j1 = this%inter_index_y(i,j) - np
3492 j2 = this%inter_index_y(i,j) + np
3493 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3495 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3496 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3505 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3506 this%trans%sub_type ==
'stddevnm1')
THEN
3508 IF (this%trans%sub_type ==
'stddev')
THEN
3517 DO j = 1, this%outny
3518 DO i = 1, this%outnx
3519 IF (
c_e(this%inter_index_x(i,j)))
THEN
3520 i1 = this%inter_index_x(i,j) - np
3521 i2 = this%inter_index_x(i,j) + np
3522 j1 = this%inter_index_y(i,j) - np
3523 j2 = this%inter_index_y(i,j) + np
3526 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3527 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3528 this%stencil(:,:), (/ns/)), nm1=nm1)
3535 ELSE IF (this%trans%sub_type ==
'max')
THEN
3540 DO j = 1, this%outny
3541 DO i = 1, this%outnx
3542 IF (
c_e(this%inter_index_x(i,j)))
THEN
3543 i1 = this%inter_index_x(i,j) - np
3544 i2 = this%inter_index_x(i,j) + np
3545 j1 = this%inter_index_y(i,j) - np
3546 j2 = this%inter_index_y(i,j) + np
3547 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3549 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3550 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3558 ELSE IF (this%trans%sub_type ==
'min')
THEN
3563 DO j = 1, this%outny
3564 DO i = 1, this%outnx
3565 IF (
c_e(this%inter_index_x(i,j)))
THEN
3566 i1 = this%inter_index_x(i,j) - np
3567 i2 = this%inter_index_x(i,j) + np
3568 j1 = this%inter_index_y(i,j) - np
3569 j2 = this%inter_index_y(i,j) + np
3570 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3572 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3573 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3581 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3586 DO j = 1, this%outny
3587 DO i = 1, this%outnx
3588 IF (
c_e(this%inter_index_x(i,j)))
THEN
3589 i1 = this%inter_index_x(i,j) - np
3590 i2 = this%inter_index_x(i,j) + np
3591 j1 = this%inter_index_y(i,j) - np
3592 j2 = this%inter_index_y(i,j) + np
3595 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3596 (/real(this%trans%stat_info%percentile)/), &
3597 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3598 this%stencil(:,:), (/ns/)))
3605 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3610 DO j = 1, this%outny
3611 DO i = 1, this%outnx
3612 IF (
c_e(this%inter_index_x(i,j)))
THEN
3613 i1 = this%inter_index_x(i,j) - np
3614 i2 = this%inter_index_x(i,j) + np
3615 j1 = this%inter_index_y(i,j) - np
3616 j2 = this%inter_index_y(i,j) + np
3617 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3619 field_out(i,j,k) = sum(interval_info_valid( &
3620 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3621 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3631ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
3634 WHERE(
c_e(this%inter_index_x(:,:)))
3635 field_out(:,:,k) = real(this%inter_index_x(:,:))
3639ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
3641 IF (this%trans%sub_type ==
'all')
THEN
3643 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3645 ELSE IF (this%trans%sub_type ==
'coordbb' .OR. this%trans%sub_type ==
'poly' &
3646 .OR. this%trans%sub_type ==
'mask')
THEN
3650 field_out(:,1,k) = pack(field_in(:,:,k),
c_e(this%point_index(:,:)))
3653 ELSE IF (this%trans%sub_type ==
'maskvalid' .OR. &
3654 this%trans%sub_type ==
'maskinvalid')
THEN
3657 WHERE (this%point_mask(:,:))
3658 field_out(:,:,k) = field_in(:,:,k)
3662 ELSE IF (this%trans%sub_type ==
'lemaskinvalid')
THEN
3665 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3666 field_out(:,:,k) = field_in(:,:,k)
3668 field_out(:,:,k) = rmiss
3672 ELSE IF (this%trans%sub_type ==
'ltmaskinvalid')
THEN
3675 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3676 field_out(:,:,k) = field_in(:,:,k)
3678 field_out(:,:,k) = rmiss
3682 ELSE IF (this%trans%sub_type ==
'gemaskinvalid')
THEN
3685 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3686 field_out(:,:,k) = field_in(:,:,k)
3688 field_out(:,:,k) = rmiss
3692 ELSE IF (this%trans%sub_type ==
'gtmaskinvalid')
THEN
3695 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3696 field_out(:,:,k) = field_in(:,:,k)
3698 field_out(:,:,k) = rmiss
3702 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
3705 WHERE (
c_e(field_in(:,:,k)))
3706 field_out(:,:,k) = field_in(:,:,k)
3708 field_out(:,:,k) = this%val1
3712 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
3713 IF (
c_e(this%val1) .AND.
c_e(this%val2))
THEN
3714 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3715 .AND. field_in(:,:,:) <= this%val2)
3716 field_out(:,:,:) = rmiss
3718 field_out(:,:,:) = field_in(:,:,:)
3720 ELSE IF (
c_e(this%val1))
THEN
3721 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3722 field_out(:,:,:) = rmiss
3724 field_out(:,:,:) = field_in(:,:,:)
3726 ELSE IF (
c_e(this%val2))
THEN
3727 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3728 field_out(:,:,:) = rmiss
3730 field_out(:,:,:) = field_in(:,:,:)
3735ELSE IF (this%trans%trans_type ==
'vertint')
THEN
3737 IF (this%trans%sub_type ==
'linear')
THEN
3739 alloc_coord_3d_in_act = .false.
3740 IF (
ASSOCIATED(this%inter_index_z))
THEN
3743 IF (
c_e(this%inter_index_z(k)))
THEN
3744 z1 = real(this%inter_zp(k))
3745 z2 = real(1.0d0 - this%inter_zp(k))
3746 SELECT CASE(vartype)
3751 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3752 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3753 field_out(i,j,k) = &
3754 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3755 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3763 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3764 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3765 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3766 field_in(i,j,this%inter_index_z(k)+1) > 0.)
THEN
3767 field_out(i,j,k) = exp( &
3768 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3769 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3777 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3778 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3779 field_out(i,j,k) = &
3780 field_in(i,j,this%inter_index_z(k))*z2 + &
3781 field_in(i,j,this%inter_index_z(k)+1)*z1
3793 IF (
ALLOCATED(this%coord_3d_in))
THEN
3794 coord_3d_in_act => this%coord_3d_in
3796 CALL l4f_category_log(this%category,l4f_debug, &
3797 "Using external vertical coordinate for vertint")
3800 IF (
PRESENT(coord_3d_in))
THEN
3802 CALL l4f_category_log(this%category,l4f_debug, &
3803 "Using internal vertical coordinate for vertint")
3805 IF (this%dolog)
THEN
3806 ALLOCATE(coord_3d_in_act(
SIZE(coord_3d_in,1), &
3807 SIZE(coord_3d_in,2),
SIZE(coord_3d_in,3)))
3808 alloc_coord_3d_in_act = .true.
3809 WHERE (
c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3810 coord_3d_in_act = log(coord_3d_in)
3812 coord_3d_in_act = rmiss
3815 coord_3d_in_act => coord_3d_in
3818 CALL l4f_category_log(this%category,l4f_error, &
3819 "Missing internal and external vertical coordinate for vertint")
3825 inused =
SIZE(coord_3d_in_act,3)
3826 IF (inused < 2)
RETURN
3830 IF (.NOT.
c_e(this%vcoord_out(k))) cycle
3834 DO kk = 1, max(inused-kkcache-1, kkcache)
3836 kkdown = kkcache - kk + 1
3838 IF (kkdown >= 1)
THEN
3839 IF (this%vcoord_out(k) >= &
3840 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3841 this%vcoord_out(k) < &
3842 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)))
THEN
3845 kfound = kkcache + this%levshift
3849 IF (kkup < inused)
THEN
3850 IF (this%vcoord_out(k) >= &
3851 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3852 this%vcoord_out(k) < &
3853 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)))
THEN
3856 kfound = kkcache + this%levshift
3863 IF (
c_e(kfound))
THEN
3864 IF (
c_e(field_in(i,j,kfound)) .AND.
c_e(field_in(i,j,kfound+1)))
THEN
3865 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3866 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3868 SELECT CASE(vartype)
3871 field_out(i,j,k) = &
3872 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3874 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3875 log(field_in(i,j,kfound+1))*z1)
3878 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3887 IF (alloc_coord_3d_in_act)
DEALLOCATE(coord_3d_in_act)
3890 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN
3893 IF (.NOT.
ASSOCIATED(this%vcoord_in) .AND. .NOT.
PRESENT(coord_3d_in))
THEN
3894 CALL l4f_category_log(this%category,l4f_error, &
3895 "linearsparse interpolation, no input vert coord available")
3899 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3903 IF (
ASSOCIATED(this%vcoord_in))
THEN
3904 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3905 .AND.
c_e(this%vcoord_in(:))
3907 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3908 .AND.
c_e(coord_3d_in(i,j,:))
3911 IF (vartype == var_press)
THEN
3912 mask_in(:) = mask_in(:) .AND. &
3913 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3915 inused = count(mask_in)
3916 IF (inused > 1)
THEN
3917 IF (
ASSOCIATED(this%vcoord_in))
THEN
3918 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3920 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3922 IF (vartype == var_press)
THEN
3923 val_in(1:inused) = log(pack( &
3924 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3927 val_in(1:inused) = pack( &
3928 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3935 DO kk = 1, max(inused-kkcache-1, kkcache)
3937 kkdown = kkcache - kk + 1
3939 IF (kkdown >= 1)
THEN
3940 IF (this%vcoord_out(k) >= &
3941 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
3942 this%vcoord_out(k) < &
3943 max(coord_in(kkdown), coord_in(kkdown+1)))
THEN
3950 IF (kkup < inused)
THEN
3951 IF (this%vcoord_out(k) >= &
3952 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
3953 this%vcoord_out(k) < &
3954 max(coord_in(kkup), coord_in(kkup+1)))
THEN
3964 IF (
c_e(kfound))
THEN
3965 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
3966 (coord_in(kfound) - coord_in(kfound-1)))
3968 IF (vartype == var_dir360)
THEN
3969 field_out(i,j,k) = &
3970 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
3971 ELSE IF (vartype == var_press)
THEN
3972 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
3974 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
3984 DEALLOCATE(coord_in,val_in)
3989ELSE IF (this%trans%trans_type ==
'' .OR. this%trans%trans_type ==
'none')
THEN
3991 field_out(:,:,:) = field_in(:,:,:)
4001FUNCTION interp_var_360(v1, v2, w1, w2)
4002REAL,
INTENT(in) :: v1, v2, w1, w2
4003REAL :: interp_var_360
4007IF (abs(v1 - v2) > 180.)
THEN
4015 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4017 interp_var_360 = v1*w2 + v2*w1
4020END FUNCTION interp_var_360
4023RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask)
RESULT(prev)
4024REAL,
INTENT(in) :: v1(:,:)
4025REAL,
INTENT(in) :: l, h, res
4026LOGICAL,
INTENT(in),
OPTIONAL :: mask(:,:)
4034IF ((h - l) <= res)
THEN
4039IF (
PRESENT(mask))
THEN
4040 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4041 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4043 nl = count(v1 >= l .AND. v1 < m)
4044 nh = count(v1 >= m .AND. v1 < h)
4047 prev = find_prevailing_direction(v1, m, h, res)
4048ELSE IF (nl > nh)
THEN
4049 prev = find_prevailing_direction(v1, l, m, res)
4050ELSE IF (nl == 0 .AND. nh == 0)
THEN
4056END FUNCTION find_prevailing_direction
4059END SUBROUTINE grid_transform_compute
4067SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4070REAL,
INTENT(in) :: field_in(:,:)
4071REAL,
INTENT(out):: field_out(:,:,:)
4072TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
4073REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
4075real,
allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4076INTEGER :: inn_p, ier, k
4077INTEGER :: innx, inny, innz, outnx, outny, outnz
4080call l4f_category_log(this%category,l4f_debug,
"start v7d_grid_transform_compute")
4083field_out(:,:,:) = rmiss
4085IF (.NOT.this%valid)
THEN
4086 call l4f_category_log(this%category,l4f_error, &
4087 "refusing to perform a non valid transformation")
4092innx =
SIZE(field_in,1); inny = 1; innz =
SIZE(field_in,2)
4093outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
4096IF (this%trans%trans_type ==
'vertint')
THEN
4098 IF (innz /= this%innz)
THEN
4099 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4100 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4101 t2c(this%innz)//
" /= "//t2c(innz))
4106 IF (outnz /= this%outnz)
THEN
4107 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4108 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4109 t2c(this%outnz)//
" /= "//t2c(outnz))
4114 IF (innx /= outnx .OR. inny /= outny)
THEN
4115 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4116 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
4117 t2c(innx)//
","//t2c(inny)//
" /= "//&
4118 t2c(outnx)//
","//t2c(outny))
4125 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
4126 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4127 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4128 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
4129 t2c(innx)//
","//t2c(inny))
4134 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
4135 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4136 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4137 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
4138 t2c(outnx)//
","//t2c(outny))
4143 IF (innz /= outnz)
THEN
4144 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4145 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
4146 t2c(innz)//
" /= "//t2c(outnz))
4154 call l4f_category_log(this%category,l4f_debug, &
4155 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//
':'// &
4156 trim(this%trans%sub_type))
4159IF (this%trans%trans_type ==
'inter')
THEN
4161 IF (this%trans%sub_type ==
'linear')
THEN
4163#ifdef HAVE_LIBNGMATH
4165 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4167 inn_p = count(
c_e(field_in(:,k)))
4169 CALL l4f_category_log(this%category,l4f_info, &
4170 "Number of sparse data points: "//t2c(inn_p)//
','//t2c(
SIZE(field_in(:,k))))
4174 field_in_p(1:inn_p) = pack(field_in(:,k),
c_e(field_in(:,k)))
4175 x_in_p(1:inn_p) = pack(this%inter_xp(:,1),
c_e(field_in(:,k)))
4176 y_in_p(1:inn_p) = pack(this%inter_yp(:,1),
c_e(field_in(:,k)))
4178 IF (.NOT.this%trans%extrap)
THEN
4179 CALL nnseti(
'ext', 0)
4180 CALL nnsetr(
'nul', rmiss)
4183 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4184 this%outnx, this%outny, real(this%inter_x(:,1)), &
4185 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4188 CALL l4f_category_log(this%category,l4f_error, &
4189 "Error code from NCAR natgrids: "//t2c(ier))
4195 CALL l4f_category_log(this%category,l4f_info, &
4196 "insufficient data in gridded region to triangulate")
4200 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4203 CALL l4f_category_log(this%category,l4f_error, &
4204 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4211ELSE IF (this%trans%trans_type ==
'boxinter' .OR. &
4212 this%trans%trans_type ==
'polyinter' .OR. &
4213 this%trans%trans_type ==
'vertint' .OR. &
4214 this%trans%trans_type ==
'metamorphosis')
THEN
4217 reshape(field_in, (/
SIZE(field_in,1), 1,
SIZE(field_in,2)/)), field_out, var, &
4222END SUBROUTINE grid_transform_v7d_grid_compute
4237ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
RESULT(zp)
4238REAL,
INTENT(in) :: z1,z2,z3,z4
4239DOUBLE PRECISION,
INTENT(in):: x1,y1
4240DOUBLE PRECISION,
INTENT(in):: x3,y3
4241DOUBLE PRECISION,
INTENT(in):: xp,yp
4248p2 = real((yp-y1)/(y3-y1))
4249p1 = real((xp-x1)/(x3-x1))
4260FUNCTION shapiro (z,zp)
RESULT(zs)
4263REAL,
INTENT(in) :: z(:)
4264REAL,
INTENT(in) :: zp
4270 zs = zp - 0.5* ( n*zp - sum(z,
c_e(z)) )/n
4279SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
4280 lon, lat, extrap, index_x, index_y)
4282logical,
INTENT(in) :: near
4283INTEGER,
INTENT(in) :: nx,ny
4284DOUBLE PRECISION,
INTENT(in) :: xmin, xmax, ymin, ymax
4285DOUBLE PRECISION,
INTENT(in) :: lon(:,:),lat(:,:)
4286LOGICAL,
INTENT(in) :: extrap
4287INTEGER,
INTENT(out) :: index_x(:,:),index_y(:,:)
4290DOUBLE PRECISION :: x(SIZE(lon,1),SIZE(lon,2)),y(SIZE(lon,1),SIZE(lon,2))
4293 CALL proj(this,lon,lat,x,y)
4294 index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4295 index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4299 CALL proj(this,lon,lat,x,y)
4300 index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4301 index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4307 index_x = max(index_x, 1)
4308 index_y = max(index_y, 1)
4309 index_x = min(index_x, lnx)
4310 index_y = min(index_y, lny)
4312 WHERE(index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
4318END 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.
Compute the distance in m between two points.
Methods for returning the value of object members.
Determine whether a point lies inside a polygon or a rectangle.
Generic subroutine for checking OPTIONAL parameters.
Compute the average of the random variable provided, taking into account missing data.
Compute a set of percentiles for a random variable.
Compute the standard deviation of the random variable provided, taking into account missing data.
Represent data in a pretty string.
This module defines usefull general purpose function and subroutine.
This module defines objects describing georeferenced sparse points possibly with topology and project...
Module for describing geographically referenced regular grids.
Module for defining the extension and coordinates of a rectangular georeferenced grid.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Module for basic statistical computations taking into account missing data.
Classe per la gestione di un volume completo di dati osservati.
Derived type defining a dynamically extensible array of TYPE(georef_coord_array) elements.
Derive type defining a single georeferenced point, either in geodetic or in projected coordinates.
This object completely describes a grid on a geographic projection.
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...