|
◆ grid_transform_compute()
recursive subroutine grid_transform_class::grid_transform_compute |
( |
type(grid_transform), intent(in), target |
this, |
|
|
real, dimension(:,:,:), intent(in) |
field_in, |
|
|
real, dimension(:,:,:), intent(out) |
field_out, |
|
|
type(vol7d_var), intent(in), optional |
var, |
|
|
real, dimension(:,:,:), intent(in), optional, target |
coord_3d_in |
|
) |
| |
|
private |
Compute the output data array from input data array according to the defined transformation.
The grid_transform object this must have been properly initialised, so that it contains all the information needed for computing the transformation. This is the grid-to-grid and grid-to-sparse points version. In the case of grid-to-sparse points transformation, the output argument is still a 2-d array with shape (np,1), which may have to be reshaped and assigned to the target 1-d array after the subroutine call by means of the RESHAPE() intrinsic function. - Parametri
-
[in] | this | grid_transformation object |
[in] | field_in | input array |
[out] | field_out | output array |
[in] | var | physical variable to be interpolated, if provided, some ad-hoc algorithms may be used where possible |
[in] | coord_3d_in | input vertical coordinate for vertical interpolation, if not provided by other means |
Definizione alla linea 3083 del file grid_transform_class.F90.
3090 navg = this%trans%box_info%npx*this%trans%box_info%npy
3093 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3094 je = j+this%trans%box_info%npy-1
3097 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3098 ie = i+this%trans%box_info%npx-1
3100 field_out(ii,jj,k) = stat_stddev( &
3101 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3106 ELSE IF (this%trans%sub_type == 'max') THEN
3109 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3110 je = j+this%trans%box_info%npy-1
3113 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3114 ie = i+this%trans%box_info%npx-1
3116 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3118 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3119 mask=(field_in(i:ie,j:je,k) /= rmiss))
3125 ELSE IF (this%trans%sub_type == 'min') THEN
3128 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3129 je = j+this%trans%box_info%npy-1
3132 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3133 ie = i+this%trans%box_info%npx-1
3135 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3137 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3138 mask=(field_in(i:ie,j:je,k) /= rmiss))
3144 ELSE IF (this%trans%sub_type == 'percentile') THEN
3146 navg = this%trans%box_info%npx*this%trans%box_info%npy
3149 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3150 je = j+this%trans%box_info%npy-1
3153 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3154 ie = i+this%trans%box_info%npx-1
3156 field_out(ii:ii,jj,k) = stat_percentile( &
3157 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3158 (/real(this%trans%stat_info%percentile)/))
3163 ELSE IF (this%trans%sub_type == 'frequency') THEN
3167 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3168 je = j+this%trans%box_info%npy-1
3171 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3172 ie = i+this%trans%box_info%npx-1
3174 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3176 field_out(ii,jj,k) = sum(interval_info_valid( &
3177 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3178 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3186 ELSE IF (this%trans%trans_type == 'inter') THEN
3188 IF (this%trans%sub_type == 'near') THEN
3191 DO j = 1, this%outny
3192 DO i = 1, this%outnx
3194 IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3195 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3201 ELSE IF (this%trans%sub_type == 'bilin') THEN
3204 DO j = 1, this%outny
3205 DO i = 1, this%outnx
3207 IF (c_e(this%inter_index_x(i,j))) THEN
3209 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3210 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j
3211 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j
3212 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+
3214 IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN
3216 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y
3217 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y
3218 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y
3219 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y
3221 xp=this%inter_xp(i,j)
3222 yp=this%inter_yp(i,j)
3224 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3232 ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
3234 DO j = 1, this%outny
3235 DO i = 1, this%outnx
3237 IF (c_e(this%inter_index_x(i,j))) THEN
3239 IF(this%inter_index_x(i,j)-1>0) THEN
3240 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y
3244 IF(this%inter_index_x(i,j)+1<=this%outnx) THEN
3245 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y
3249 IF(this%inter_index_y(i,j)+1<=this%outny) THEN
3250 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y
3254 IF(this%inter_index_y(i,j)-1>0) THEN
3255 z(4)=field_in(this%inter_index_x(i,j) ,this%inter_index_y
3259 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j
3268 ELSE IF (this%trans%trans_type == 'boxinter' &
3269 .OR. this%trans%trans_type == 'polyinter' &
3270 .OR. this%trans%trans_type == 'maskinter') THEN
3272 IF (this%trans%sub_type == 'average') THEN
3274 IF (vartype == var_dir360) THEN
3276 DO j = 1, this%outny
3277 DO i = 1, this%outnx
3278 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k)
3280 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(
3286 ALLOCATE(nval(this%outnx, this%outny))
3287 field_out(:,:,:) = 0.0
3292 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3293 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3294 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j)
3296 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3297 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3301 WHERE (nval(:,:) /= 0)
3302 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3304 field_out(:,:,k) = rmiss
3310 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3311 this%trans%sub_type == 'stddevnm1') THEN
3313 IF (this%trans%sub_type == 'stddev') THEN
3319 DO j = 1, this%outny
3320 DO i = 1, this%outnx
3322 field_out(i:i,j,k) = stat_stddev( &
3323 reshape(field_in(:,:,k), (/ SIZE(field_in(:,:,k))/)), &
3324 mask=reshape((this%inter_index_x == i .AND. &
3325 this%inter_index_y == j), (/ SIZE(field_in(:,:,k))/)), nm1=nm1
3330 ELSE IF (this%trans%sub_type == 'max') THEN
3335 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3336 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y THEN
3337 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3338 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i
3341 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3350 ELSE IF (this%trans%sub_type == 'min') THEN
3355 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3356 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y THEN
3357 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3358 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i
3361 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3369 ELSE IF (this%trans%sub_type == 'percentile') THEN
3372 DO j = 1, this%outny
3373 DO i = 1, this%outnx
3375 field_out(i:i,j,k) = stat_percentile( &
3376 reshape(field_in(:,:,k), (/ SIZE(field_in(:,:,k))/)), &
3377 (/real(this%trans%stat_info%percentile)/), &
3378 mask=reshape((this%inter_index_x == i .AND. &
3379 this%inter_index_y == j), (/ SIZE(field_in(:,:,k))/)))
3384 ELSE IF (this%trans%sub_type == 'frequency') THEN
3386 ALLOCATE(nval(this%outnx, this%outny))
3387 field_out(:,:,:) = 0.0
3392 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3393 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3394 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3395 interval_info_valid(this%trans%interval_info, field_in(i,j,k
3396 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3397 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3401 WHERE (nval(:,:) /= 0)
3402 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3404 field_out(:,:,k) = rmiss
3411 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
3412 np = SIZE(this%stencil,1)/2
3413 ns = SIZE(this%stencil)
3415 IF (this%trans%sub_type == 'average') THEN
3417 IF (vartype == var_dir360) THEN
3419 DO j = 1, this%outny
3420 DO i = 1, this%outnx
3421 IF (c_e(this%inter_index_x(i,j))) THEN
3422 i1 = this%inter_index_x(i,j) - np
3423 i2 = this%inter_index_x(i,j) + np
3424 j1 = this%inter_index_y(i,j) - np
3425 j2 = this%inter_index_y(i,j) + np
3426 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2
3428 mask=this%stencil(:,:))
3439 DO j = 1, this%outny
3440 DO i = 1, this%outnx
3441 IF (c_e(this%inter_index_x(i,j))) THEN
3442 i1 = this%inter_index_x(i,j) - np
3443 i2 = this%inter_index_x(i,j) + np
3444 j1 = this%inter_index_y(i,j) - np
3445 j2 = this%inter_index_y(i,j) + np
3446 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil
3448 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3449 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil
3458 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3459 this%trans%sub_type == 'stddevnm1') THEN
3461 IF (this%trans%sub_type == 'stddev') THEN
3470 DO j = 1, this%outny
3471 DO i = 1, this%outnx
3472 IF (c_e(this%inter_index_x(i,j))) THEN
3473 i1 = this%inter_index_x(i,j) - np
3474 i2 = this%inter_index_x(i,j) + np
3475 j1 = this%inter_index_y(i,j) - np
3476 j2 = this%inter_index_y(i,j) + np
3478 field_out(i:i,j,k) = stat_stddev( &
3479 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3480 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3481 this%stencil(:,:), (/ns/)), nm1=nm1)
3488 ELSE IF (this%trans%sub_type == 'max') THEN
3493 DO j = 1, this%outny
3494 DO i = 1, this%outnx
3495 IF (c_e(this%inter_index_x(i,j))) THEN
3496 i1 = this%inter_index_x(i,j) - np
3497 i2 = this%inter_index_x(i,j) + np
3498 j1 = this%inter_index_y(i,j) - np
3499 j2 = this%inter_index_y(i,j) + np
3500 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil
3502 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3503 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(
3511 ELSE IF (this%trans%sub_type == 'min') THEN
3516 DO j = 1, this%outny
3517 DO i = 1, this%outnx
3518 IF (c_e(this%inter_index_x(i,j))) THEN
3519 i1 = this%inter_index_x(i,j) - np
3520 i2 = this%inter_index_x(i,j) + np
3521 j1 = this%inter_index_y(i,j) - np
3522 j2 = this%inter_index_y(i,j) + np
3523 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil
3525 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3526 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(
3534 ELSE IF (this%trans%sub_type == 'percentile') THEN
3539 DO j = 1, this%outny
3540 DO i = 1, this%outnx
3541 IF (c_e(this%inter_index_x(i,j))) THEN
3542 i1 = this%inter_index_x(i,j) - np
3543 i2 = this%inter_index_x(i,j) + np
3544 j1 = this%inter_index_y(i,j) - np
3545 j2 = this%inter_index_y(i,j) + np
3547 field_out(i:i,j,k) = stat_percentile( &
3548 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3549 (/real(this%trans%stat_info%percentile)/), &
3550 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3551 this%stencil(:,:), (/ns/)))
3558 ELSE IF (this%trans%sub_type == 'frequency') THEN
3563 DO j = 1, this%outny
3564 DO i = 1, this%outnx
3565 IF (c_e(this%inter_index_x(i,j))) THEN
3566 i1 = this%inter_index_x(i,j) - np
3567 i2 = this%inter_index_x(i,j) + np
3568 j1 = this%inter_index_y(i,j) - np
3569 j2 = this%inter_index_y(i,j) + np
3570 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil
3572 field_out(i,j,k) = sum(interval_info_valid( &
3573 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3574 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(
3584 ELSE IF (this%trans%trans_type == 'maskgen') THEN
3587 WHERE(c_e(this%inter_index_x(:,:)))
3588 field_out(:,:,k) = real(this%inter_index_x(:,:))
3592 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
3594 IF (this%trans%sub_type == 'all') THEN
3596 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz
3598 ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly'
3599 .OR. this%trans%sub_type == 'mask') THEN
3603 field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)
3606 ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
3607 this%trans%sub_type == 'maskinvalid') THEN
3610 WHERE (this%point_mask(:,:))
3611 field_out(:,:,k) = field_in(:,:,k)
3615 ELSE IF (this%trans%sub_type == 'lemaskinvalid') THEN
3618 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(
3619 field_out(:,:,k) = field_in(:,:,k)
3621 field_out(:,:,k) = rmiss
3625 ELSE IF (this%trans%sub_type == 'ltmaskinvalid') THEN
3628 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask
3629 field_out(:,:,k) = field_in(:,:,k)
3631 field_out(:,:,k) = rmiss
3635 ELSE IF (this%trans%sub_type == 'gemaskinvalid') THEN
3638 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(
3639 field_out(:,:,k) = field_in(:,:,k)
3641 field_out(:,:,k) = rmiss
3645 ELSE IF (this%trans%sub_type == 'gtmaskinvalid') THEN
3648 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask
3649 field_out(:,:,k) = field_in(:,:,k)
3651 field_out(:,:,k) = rmiss
3655 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
3658 WHERE (c_e(field_in(:,:,k)))
3659 field_out(:,:,k) = field_in(:,:,k)
3661 field_out(:,:,k) = this%val1
3665 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
3666 IF (c_e(this%val1) .AND. c_e(this%val2)) THEN
3667 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3668 .AND. field_in(:,:,:) <= this%val2)
3669 field_out(:,:,:) = rmiss
3671 field_out(:,:,:) = field_in(:,:,:)
3673 ELSE IF (c_e(this%val1)) THEN
3674 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3675 field_out(:,:,:) = rmiss
3677 field_out(:,:,:) = field_in(:,:,:)
3679 ELSE IF (c_e(this%val2)) THEN
3680 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3681 field_out(:,:,:) = rmiss
3683 field_out(:,:,:) = field_in(:,:,:)
3688 ELSE IF (this%trans%trans_type == 'vertint') THEN
3690 IF (this%trans%sub_type == 'linear') THEN
3692 alloc_coord_3d_in_act = .false.
3693 IF ( ASSOCIATED(this%inter_index_z)) THEN
3696 IF (c_e(this%inter_index_z(k))) THEN
3697 z1 = real(this%inter_zp(k))
3698 z2 = real(1.0d0 - this%inter_zp(k))
3699 SELECT CASE(vartype)
3704 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3705 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3706 field_out(i,j,k) = &
3707 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3708 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3716 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3717 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3718 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3719 field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN
3720 field_out(i,j,k) = exp( &
3721 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3722 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3730 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3731 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3732 field_out(i,j,k) = &
3733 field_in(i,j,this%inter_index_z(k))*z2 + &
3734 field_in(i,j,this%inter_index_z(k)+1)*z1
3746 IF ( ALLOCATED(this%coord_3d_in)) THEN
3747 coord_3d_in_act => this%coord_3d_in
3749 CALL l4f_category_log(this%category,l4f_debug, &
3750 "Using external vertical coordinate for vertint")
3753 IF ( PRESENT(coord_3d_in)) THEN
3755 CALL l4f_category_log(this%category,l4f_debug, &
3756 "Using internal vertical coordinate for vertint")
3758 IF (this%dolog) THEN
3759 ALLOCATE(coord_3d_in_act( SIZE(coord_3d_in,1), &
3760 SIZE(coord_3d_in,2), SIZE(coord_3d_in,3)))
3761 alloc_coord_3d_in_act = .true.
3762 WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3763 coord_3d_in_act = log(coord_3d_in)
3765 coord_3d_in_act = rmiss
3768 coord_3d_in_act => coord_3d_in
3771 CALL l4f_category_log(this%category,l4f_error, &
3772 "Missing internal and external vertical coordinate for vertint"
3778 inused = SIZE(coord_3d_in_act,3)
3779 IF (inused < 2) RETURN
3783 IF (.NOT.c_e(this%vcoord_out(k))) cycle
3787 DO kk = 1, max(inused-kkcache-1, kkcache)
3789 kkdown = kkcache - kk + 1
3791 IF (kkdown >= 1) THEN
3792 IF (this%vcoord_out(k) >= &
3793 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown
3794 this%vcoord_out(k) < &
3795 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown THEN
3798 kfound = kkcache + this%levshift
3802 IF (kkup < inused) THEN
3803 IF (this%vcoord_out(k) >= &
3804 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup
3805 this%vcoord_out(k) < &
3806 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup THEN
3809 kfound = kkcache + this%levshift
3816 IF (c_e(kfound)) THEN
3817 IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound THEN
3818 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin
3819 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin
3821 SELECT CASE(vartype)
3824 field_out(i,j,k) = &
3825 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound
3827 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 +
3828 log(field_in(i,j,kfound+1))*z1)
3831 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i
3840 IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
3843 ELSE IF (this%trans%sub_type == 'linearsparse') THEN
3846 IF (.NOT. ASSOCIATED(this%vcoord_in) .AND. .NOT. PRESENT(coord_3d_in)) THEN
3847 CALL l4f_category_log(this%category,l4f_error, &
3848 "linearsparse interpolation, no input vert coord available")
3852 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3856 IF ( ASSOCIATED(this%vcoord_in)) THEN
3857 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused
3858 .AND. c_e(this%vcoord_in(:))
3860 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused
3861 .AND. c_e(coord_3d_in(i,j,:))
3864 IF (vartype == var_press) THEN
3865 mask_in(:) = mask_in(:) .AND. &
3866 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0
3868 inused = count(mask_in)
3869 IF (inused > 1) THEN
3870 IF ( ASSOCIATED(this%vcoord_in)) THEN
3871 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3873 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3875 IF (vartype == var_press) THEN
3876 val_in(1:inused) = log(pack( &
3877 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3880 val_in(1:inused) = pack( &
3881 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3888 DO kk = 1, max(inused-kkcache-1, kkcache)
3890 kkdown = kkcache - kk + 1
3892 IF (kkdown >= 1) THEN
3893 IF (this%vcoord_out(k) >= &
3894 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
3895 this%vcoord_out(k) < &
3896 max(coord_in(kkdown), coord_in(kkdown+1))) THEN
3903 IF (kkup < inused) THEN
3904 IF (this%vcoord_out(k) >= &
3905 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
3906 this%vcoord_out(k) < &
3907 max(coord_in(kkup), coord_in(kkup+1))) THEN
3917 IF (c_e(kfound)) THEN
3918 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
3919 (coord_in(kfound) - coord_in(kfound-1)))
3921 IF (vartype == var_dir360) THEN
3922 field_out(i,j,k) = &
3923 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2
3924 ELSE IF (vartype == var_press) THEN
3925 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound
3927 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
3937 DEALLOCATE(coord_in,val_in)
3942 ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none'THEN
3944 field_out(:,:,:) = field_in(:,:,:)
3954 FUNCTION interp_var_360(v1, v2, w1, w2)
3955 REAL, INTENT(in) :: v1, v2, w1, w2
3956 REAL :: interp_var_360
3960 IF (abs(v1 - v2) > 180.) THEN
3968 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
3970 interp_var_360 = v1*w2 + v2*w1
3973 END FUNCTION interp_var_360
3976 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
3977 REAL, INTENT(in) :: v1(:,:)
3978 REAL, INTENT(in) :: l, h, res
3979 LOGICAL, INTENT(in), OPTIONAL :: mask(:,:)
3987 IF ((h - l) <= res) THEN
3992 IF ( PRESENT(mask)) THEN
3993 nl = count(v1 >= l .AND. v1 < m .AND. mask)
3994 nh = count(v1 >= m .AND. v1 < h .AND. mask)
3996 nl = count(v1 >= l .AND. v1 < m)
3997 nh = count(v1 >= m .AND. v1 < h)
4000 prev = find_prevailing_direction(v1, m, h, res)
4001 ELSE IF (nl > nh) THEN
4002 prev = find_prevailing_direction(v1, l, m, res)
4003 ELSE IF (nl == 0 .AND. nh == 0) THEN
4009 END FUNCTION find_prevailing_direction
4012 END SUBROUTINE grid_transform_compute
4020 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4022 TYPE(grid_transform), INTENT(in) :: this
4023 REAL, INTENT(in) :: field_in(:,:)
4024 REAL, INTENT(out):: field_out(:,:,:)
4025 TYPE(vol7d_var), INTENT(in), OPTIONAL :: var
4026 REAL, INTENT(in), OPTIONAL, TARGET :: coord_3d_in(:,:,:)
4028 real, allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4029 INTEGER :: inn_p, ier, k
4030 INTEGER :: innx, inny, innz, outnx, outny, outnz
4036 field_out(:,:,:) = rmiss
4038 IF (.NOT.this%valid) THEN
4039 call l4f_category_log(this%category,l4f_error, &
4040 "refusing to perform a non valid transformation")
4045 innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
4046 outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out
4049 IF (this%trans%trans_type == 'vertint') THEN
4051 IF (innz /= this%innz) THEN
4052 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation"
4053 CALL l4f_category_log(this%category,l4f_error, "inconsistent input shape: "
4054 t2c(this%innz)// " /= "//t2c(innz))
4059 IF (outnz /= this%outnz) THEN
4060 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation"
4061 CALL l4f_category_log(this%category,l4f_error, "inconsistent output shape: "
4062 t2c(this%outnz)// " /= "//t2c(outnz))
4067 IF (innx /= outnx .OR. inny /= outny) THEN
4068 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation"
4069 CALL l4f_category_log(this%category,l4f_error, "inconsistent hor. sizes: "
4070 t2c(innx)// ","//t2c(inny)// " /= "//&
4071 t2c(outnx)// ","//t2c(outny))
4076 ELSE ! horizontal interpolation
4078 IF (innx /= this%innx .OR. inny /= this%inny) THEN
4079 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation"
4080 CALL l4f_category_log(this%category,l4f_error, "inconsistent input shape: "
4081 t2c(this%innx)// ","//t2c(this%inny)// " /= "//&
4082 t2c(innx)// ","//t2c(inny))
4087 IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
4088 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation"
4089 CALL l4f_category_log(this%category,l4f_error, "inconsistent output shape: "
4090 t2c(this%outnx)// ","//t2c(this%outny)// " /= "//&
4091 t2c(outnx)// ","//t2c(outny))
4096 IF (innz /= outnz) THEN
4097 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation"
4098 CALL l4f_category_log(this%category,l4f_error, "inconsistent vert. sizes: "
4099 t2c(innz)// " /= "//t2c(outnz))
4107 call l4f_category_log(this%category,l4f_debug, &
4108 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type) ':'
4109 trim(this%trans%sub_type))
4112 IF (this%trans%trans_type == 'inter') THEN
4114 IF (this%trans%sub_type == 'linear') THEN
4116 #ifdef HAVE_LIBNGMATH
4118 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny)
4120 inn_p = count(c_e(field_in(:,k)))
4122 CALL l4f_category_log(this%category,l4f_info, &
4123 "Number of sparse data points: "//t2c(inn_p)// ','//t2c( SIZE(field_in
4127 field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k)))
4128 x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k)))
4129 y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k)))
4131 IF (.NOT.this%trans%extrap) THEN
4132 CALL nnseti( 'ext', 0)
4133 CALL nnsetr( 'nul', rmiss)
4136 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4137 this%outnx, this%outny, real(this%inter_x(:,1)), &
4138 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4141 CALL l4f_category_log(this%category,l4f_error, &
4142 "Error code from NCAR natgrids: "//t2c(ier))
4148 CALL l4f_category_log(this%category,l4f_info, &
4149 "insufficient data in gridded region to triangulate")
4153 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4156 CALL l4f_category_log(this%category,l4f_error, &
4157 "libsim compiled without natgridd(ngmath ncarg library) ")
4164 .OR. ELSE IF (this%trans%trans_type == 'boxinter' &
4165 .OR. this%trans%trans_type == 'polyinter' &
4166 .OR. this%trans%trans_type == 'vertint' &
4167 this%trans%trans_type == 'metamorphosis') THEN ! use the grid-to-grid method
4169 CALL compute(this, &
4170 RESHAPE(field_in, (/SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
4175 END SUBROUTINE grid_transform_v7d_grid_compute
4178 ! Bilinear interpolation
4179 ! Effettua interpolazione bilineare dati i valori nei punti 1,2,3,4 e
4180 ! le coordinate dei punti 1 e 3 oltre a quelle del punto p dove viene
4181 ! valutato il campo.
4182 !_____________________________________________________________
4183 ! disposizione punti
4189 ! _____________________________________________________________
4190 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
4191 REAL,INTENT(in) :: z1,z2,z3,z4 ! Z values on the four points
4192 DOUBLE PRECISION,INTENT(in):: x1,y1 ! coordinate of the lower left point
4193 DOUBLE PRECISION,INTENT(in):: x3,y3 ! coordinate of the upper right point
4194 DOUBLE PRECISION,INTENT(in):: xp,yp ! coordinate of point where interpolate
|