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