|
◆ grid_transform_compute()
recursive subroutine 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 |
|
) |
| |
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 3130 del file grid_transform_class.F90.
3134 navg = this%trans%box_info%npx*this%trans%box_info%npy
3137 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3138 je = j+this%trans%box_info%npy-1
3141 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3142 ie = i+this%trans%box_info%npx-1
3144 field_out(ii,jj,k) = stat_stddev( &
3145 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3150 ELSE IF (this%trans%sub_type == 'max') THEN
3153 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3154 je = j+this%trans%box_info%npy-1
3157 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3158 ie = i+this%trans%box_info%npx-1
3160 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3162 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3163 mask=(field_in(i:ie,j:je,k) /= rmiss))
3169 ELSE IF (this%trans%sub_type == 'min') THEN
3172 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3173 je = j+this%trans%box_info%npy-1
3176 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3177 ie = i+this%trans%box_info%npx-1
3179 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3181 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3182 mask=(field_in(i:ie,j:je,k) /= rmiss))
3188 ELSE IF (this%trans%sub_type == 'percentile') THEN
3190 navg = this%trans%box_info%npx*this%trans%box_info%npy
3193 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3194 je = j+this%trans%box_info%npy-1
3197 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3198 ie = i+this%trans%box_info%npx-1
3200 field_out(ii:ii,jj,k) = stat_percentile( &
3201 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3202 (/real(this%trans%stat_info%percentile)/))
3207 ELSE IF (this%trans%sub_type == 'frequency') THEN
3211 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3212 je = j+this%trans%box_info%npy-1
3215 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3216 ie = i+this%trans%box_info%npx-1
3218 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3220 field_out(ii,jj,k) = sum(interval_info_valid( &
3221 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3222 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3230 ELSE IF (this%trans%trans_type == 'inter') THEN
3232 IF (this%trans%sub_type == 'near') THEN
3235 DO j = 1, this%outny
3236 DO i = 1, this%outnx
3238 IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3239 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3245 ELSE IF (this%trans%sub_type == 'bilin') THEN
3248 DO j = 1, this%outny
3249 DO i = 1, this%outnx
3251 IF (c_e(this%inter_index_x(i,j))) THEN
3253 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3254 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3255 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3256 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3258 IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN
3260 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3261 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3262 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3263 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3265 xp=this%inter_xp(i,j)
3266 yp=this%inter_yp(i,j)
3268 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3276 ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
3278 DO j = 1, this%outny
3279 DO i = 1, this%outnx
3281 IF (c_e(this%inter_index_x(i,j))) THEN
3283 IF(this%inter_index_x(i,j)-1>0) THEN
3284 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3288 IF(this%inter_index_x(i,j)+1<=this%outnx) THEN
3289 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3293 IF(this%inter_index_y(i,j)+1<=this%outny) THEN
3294 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3298 IF(this%inter_index_y(i,j)-1>0) THEN
3299 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3303 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3312 ELSE IF (this%trans%trans_type == 'intersearch') THEN
3315 likethis%trans%trans_type = 'inter'
3316 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3317 CALL getenv( 'LIBSIM_DISABLEOPTSEARCH', env_var)
3318 optsearch = len_trim(env_var) == 0
3321 IF ((.NOT.all(c_e(field_out(:,:,k)))) .AND. (any(c_e(field_in(:,:,k))))) THEN
3322 DO j = 1, this%outny
3323 DO i = 1, this%outnx
3324 IF (.NOT.c_e(field_out(i,j,k))) THEN
3328 ix = this%inter_index_x(i,j)
3329 iy = this%inter_index_y(i,j)
3330 DO s = 0, max(this%innx, this%inny)
3332 DO m = iy-s, iy+s, max(2*s, 1)
3333 IF (m < 1 .OR. m > this%inny) cycle
3334 DO l = max(1, ix-s), min(this%innx, ix+s)
3335 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3336 IF (c_e(field_in(l,m,k))) THEN
3337 IF (disttmp < dist) THEN
3339 field_out(i,j,k) = field_in(l,m,k)
3341 ELSE IF (disttmp == dist) THEN
3342 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3343 nearcount = nearcount + 1
3346 IF (disttmp < dist) farenough = .false.
3349 DO m = max(1, iy-s+1), min(this%inny, iy+s-1)
3350 DO l = ix-s, ix+s, 2*s
3351 IF (l < 1 .OR. l > this%innx) cycle
3352 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3353 IF (c_e(field_in(l,m,k))) THEN
3354 IF (disttmp < dist) THEN
3356 field_out(i,j,k) = field_in(l,m,k)
3358 ELSE IF (disttmp == dist) THEN
3359 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3360 nearcount = nearcount + 1
3363 IF (disttmp < dist) farenough = .false.
3366 IF (s > 0 .AND. farenough) EXIT
3371 IF (c_e(field_in(l,m,k))) THEN
3372 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3373 IF (disttmp < dist) THEN
3375 field_out(i,j,k) = field_in(l,m,k)
3377 ELSE IF (disttmp == dist) THEN
3378 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3379 nearcount = nearcount + 1
3386 IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount
3393 ELSE IF (this%trans%trans_type == 'boxinter' &
3394 .OR. this%trans%trans_type == 'polyinter' &
3395 .OR. this%trans%trans_type == 'maskinter') THEN
3397 IF (this%trans%sub_type == 'average') THEN
3399 IF (vartype == var_dir360) THEN
3401 DO j = 1, this%outny
3402 DO i = 1, this%outnx
3403 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3405 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3411 ALLOCATE(nval(this%outnx, this%outny))
3412 field_out(:,:,:) = 0.0
3417 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3418 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3419 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3421 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3422 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3426 WHERE (nval(:,:) /= 0)
3427 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3429 field_out(:,:,k) = rmiss
3435 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3436 this%trans%sub_type == 'stddevnm1') THEN
3438 IF (this%trans%sub_type == 'stddev') THEN
3444 DO j = 1, this%outny
3445 DO i = 1, this%outnx
3447 field_out(i:i,j,k) = stat_stddev( &
3448 reshape(field_in(:,:,k), (/ SIZE(field_in(:,:,k))/)), &
3449 mask=reshape((this%inter_index_x == i .AND. &
3450 this%inter_index_y == j), (/ SIZE(field_in(:,:,k))/)), nm1=nm1)
3455 ELSE IF (this%trans%sub_type == 'max') THEN
3460 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3461 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3462 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3463 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3466 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3475 ELSE IF (this%trans%sub_type == 'min') THEN
3480 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3481 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3482 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3483 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3486 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3494 ELSE IF (this%trans%sub_type == 'percentile') THEN
3497 DO j = 1, this%outny
3498 DO i = 1, this%outnx
3500 field_out(i:i,j,k) = stat_percentile( &
3501 reshape(field_in(:,:,k), (/ SIZE(field_in(:,:,k))/)), &
3502 (/real(this%trans%stat_info%percentile)/), &
3503 mask=reshape((this%inter_index_x == i .AND. &
3504 this%inter_index_y == j), (/ SIZE(field_in(:,:,k))/)))
3509 ELSE IF (this%trans%sub_type == 'frequency') THEN
3511 ALLOCATE(nval(this%outnx, this%outny))
3512 field_out(:,:,:) = 0.0
3517 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3518 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3519 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3520 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3521 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3522 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3526 WHERE (nval(:,:) /= 0)
3527 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3529 field_out(:,:,k) = rmiss
3536 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
3537 np = SIZE(this%stencil,1)/2
3538 ns = SIZE(this%stencil)
3540 IF (this%trans%sub_type == 'average') THEN
3542 IF (vartype == var_dir360) THEN
3544 DO j = 1, this%outny
3545 DO i = 1, this%outnx
3546 IF (c_e(this%inter_index_x(i,j))) THEN
3547 i1 = this%inter_index_x(i,j) - np
3548 i2 = this%inter_index_x(i,j) + np
3549 j1 = this%inter_index_y(i,j) - np
3550 j2 = this%inter_index_y(i,j) + np
3551 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3553 mask=this%stencil(:,:))
3564 DO j = 1, this%outny
3565 DO i = 1, this%outnx
3566 IF (c_e(this%inter_index_x(i,j))) THEN
3567 i1 = this%inter_index_x(i,j) - np
3568 i2 = this%inter_index_x(i,j) + np
3569 j1 = this%inter_index_y(i,j) - np
3570 j2 = this%inter_index_y(i,j) + np
3571 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3573 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3574 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3583 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3584 this%trans%sub_type == 'stddevnm1') THEN
3586 IF (this%trans%sub_type == 'stddev') THEN
3595 DO j = 1, this%outny
3596 DO i = 1, this%outnx
3597 IF (c_e(this%inter_index_x(i,j))) THEN
3598 i1 = this%inter_index_x(i,j) - np
3599 i2 = this%inter_index_x(i,j) + np
3600 j1 = this%inter_index_y(i,j) - np
3601 j2 = this%inter_index_y(i,j) + np
3603 field_out(i:i,j,k) = stat_stddev( &
3604 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3605 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3606 this%stencil(:,:), (/ns/)), nm1=nm1)
3613 ELSE IF (this%trans%sub_type == 'max') THEN
3618 DO j = 1, this%outny
3619 DO i = 1, this%outnx
3620 IF (c_e(this%inter_index_x(i,j))) THEN
3621 i1 = this%inter_index_x(i,j) - np
3622 i2 = this%inter_index_x(i,j) + np
3623 j1 = this%inter_index_y(i,j) - np
3624 j2 = this%inter_index_y(i,j) + np
3625 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3627 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3628 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3636 ELSE IF (this%trans%sub_type == 'min') THEN
3641 DO j = 1, this%outny
3642 DO i = 1, this%outnx
3643 IF (c_e(this%inter_index_x(i,j))) THEN
3644 i1 = this%inter_index_x(i,j) - np
3645 i2 = this%inter_index_x(i,j) + np
3646 j1 = this%inter_index_y(i,j) - np
3647 j2 = this%inter_index_y(i,j) + np
3648 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3650 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3651 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3659 ELSE IF (this%trans%sub_type == 'percentile') THEN
3664 DO j = 1, this%outny
3665 DO i = 1, this%outnx
3666 IF (c_e(this%inter_index_x(i,j))) THEN
3667 i1 = this%inter_index_x(i,j) - np
3668 i2 = this%inter_index_x(i,j) + np
3669 j1 = this%inter_index_y(i,j) - np
3670 j2 = this%inter_index_y(i,j) + np
3672 field_out(i:i,j,k) = stat_percentile( &
3673 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3674 (/real(this%trans%stat_info%percentile)/), &
3675 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3676 this%stencil(:,:), (/ns/)))
3683 ELSE IF (this%trans%sub_type == 'frequency') THEN
3688 DO j = 1, this%outny
3689 DO i = 1, this%outnx
3690 IF (c_e(this%inter_index_x(i,j))) THEN
3691 i1 = this%inter_index_x(i,j) - np
3692 i2 = this%inter_index_x(i,j) + np
3693 j1 = this%inter_index_y(i,j) - np
3694 j2 = this%inter_index_y(i,j) + np
3695 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3697 field_out(i,j,k) = sum(interval_info_valid( &
3698 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3699 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3709 ELSE IF (this%trans%trans_type == 'maskgen') THEN
3712 WHERE(c_e(this%inter_index_x(:,:)))
3713 field_out(:,:,k) = real(this%inter_index_x(:,:))
3717 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
3719 IF (this%trans%sub_type == 'all') THEN
3721 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3723 ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' &
3724 .OR. this%trans%sub_type == 'mask') THEN
3728 field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)))
3731 ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
3732 this%trans%sub_type == 'maskinvalid') THEN
3735 WHERE (this%point_mask(:,:))
3736 field_out(:,:,k) = field_in(:,:,k)
3740 ELSE IF (this%trans%sub_type == 'lemaskinvalid') THEN
3743 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3744 field_out(:,:,k) = field_in(:,:,k)
3746 field_out(:,:,k) = rmiss
3750 ELSE IF (this%trans%sub_type == 'ltmaskinvalid') THEN
3753 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3754 field_out(:,:,k) = field_in(:,:,k)
3756 field_out(:,:,k) = rmiss
3760 ELSE IF (this%trans%sub_type == 'gemaskinvalid') THEN
3763 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3764 field_out(:,:,k) = field_in(:,:,k)
3766 field_out(:,:,k) = rmiss
3770 ELSE IF (this%trans%sub_type == 'gtmaskinvalid') THEN
3773 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3774 field_out(:,:,k) = field_in(:,:,k)
3776 field_out(:,:,k) = rmiss
3780 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
3783 WHERE (c_e(field_in(:,:,k)))
3784 field_out(:,:,k) = field_in(:,:,k)
3786 field_out(:,:,k) = this%val1
3790 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
3791 IF (c_e(this%val1) .AND. c_e(this%val2)) THEN
3792 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3793 .AND. field_in(:,:,:) <= this%val2)
3794 field_out(:,:,:) = rmiss
3796 field_out(:,:,:) = field_in(:,:,:)
3798 ELSE IF (c_e(this%val1)) THEN
3799 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3800 field_out(:,:,:) = rmiss
3802 field_out(:,:,:) = field_in(:,:,:)
3804 ELSE IF (c_e(this%val2)) THEN
3805 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3806 field_out(:,:,:) = rmiss
3808 field_out(:,:,:) = field_in(:,:,:)
3813 ELSE IF (this%trans%trans_type == 'vertint') THEN
3815 IF (this%trans%sub_type == 'linear') THEN
3817 alloc_coord_3d_in_act = .false.
3818 IF ( ASSOCIATED(this%inter_index_z)) THEN
3821 IF (c_e(this%inter_index_z(k))) THEN
3822 z1 = real(this%inter_zp(k))
3823 z2 = real(1.0d0 - this%inter_zp(k))
3824 SELECT CASE(vartype)
3829 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3830 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3831 field_out(i,j,k) = &
3832 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3833 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3841 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3842 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3843 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3844 field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN
3845 field_out(i,j,k) = exp( &
3846 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3847 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3855 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3856 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3857 field_out(i,j,k) = &
3858 field_in(i,j,this%inter_index_z(k))*z2 + &
3859 field_in(i,j,this%inter_index_z(k)+1)*z1
3871 IF ( ALLOCATED(this%coord_3d_in)) THEN
3872 coord_3d_in_act => this%coord_3d_in
3874 CALL l4f_category_log(this%category,l4f_debug, &
3875 "Using external vertical coordinate for vertint")
3878 IF ( PRESENT(coord_3d_in)) THEN
3880 CALL l4f_category_log(this%category,l4f_debug, &
3881 "Using internal vertical coordinate for vertint")
3883 IF (this%dolog) THEN
3884 ALLOCATE(coord_3d_in_act( SIZE(coord_3d_in,1), &
3885 SIZE(coord_3d_in,2), SIZE(coord_3d_in,3)))
3886 alloc_coord_3d_in_act = .true.
3887 WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3888 coord_3d_in_act = log(coord_3d_in)
3890 coord_3d_in_act = rmiss
3893 coord_3d_in_act => coord_3d_in
3896 CALL l4f_category_log(this%category,l4f_error, &
3897 "Missing internal and external vertical coordinate for vertint")
3903 inused = SIZE(coord_3d_in_act,3)
3904 IF (inused < 2) RETURN
3908 IF (.NOT.c_e(this%vcoord_out(k))) cycle
3912 DO kk = 1, max(inused-kkcache-1, kkcache)
3914 kkdown = kkcache - kk + 1
3916 IF (kkdown >= 1) THEN
3917 IF (this%vcoord_out(k) >= &
3918 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3919 this%vcoord_out(k) < &
3920 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1))) THEN
3923 kfound = kkcache + this%levshift
3927 IF (kkup < inused) THEN
3928 IF (this%vcoord_out(k) >= &
3929 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3930 this%vcoord_out(k) < &
3931 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1))) THEN
3934 kfound = kkcache + this%levshift
3941 IF (c_e(kfound)) THEN
3942 IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound+1))) THEN
3943 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3944 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3946 SELECT CASE(vartype)
3949 field_out(i,j,k) = &
3950 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3952 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3953 log(field_in(i,j,kfound+1))*z1)
3956 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3965 IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
3968 ELSE IF (this%trans%sub_type == 'linearsparse') THEN
3971 IF (.NOT. ASSOCIATED(this%vcoord_in) .AND. .NOT. PRESENT(coord_3d_in)) THEN
3972 CALL l4f_category_log(this%category,l4f_error, &
3973 "linearsparse interpolation, no input vert coord available")
3977 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3981 IF ( ASSOCIATED(this%vcoord_in)) THEN
3982 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3983 .AND. c_e(this%vcoord_in(:))
3985 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3986 .AND. c_e(coord_3d_in(i,j,:))
3989 IF (vartype == var_press) THEN
3990 mask_in(:) = mask_in(:) .AND. &
3991 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3993 inused = count(mask_in)
3994 IF (inused > 1) THEN
3995 IF ( ASSOCIATED(this%vcoord_in)) THEN
3996 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3998 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
4000 IF (vartype == var_press) THEN
4001 val_in(1:inused) = log(pack( &
4002 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4005 val_in(1:inused) = pack( &
4006 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4013 DO kk = 1, max(inused-kkcache-1, kkcache)
4015 kkdown = kkcache - kk + 1
4017 IF (kkdown >= 1) THEN
4018 IF (this%vcoord_out(k) >= &
4019 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
4020 this%vcoord_out(k) < &
4021 max(coord_in(kkdown), coord_in(kkdown+1))) THEN
4028 IF (kkup < inused) THEN
4029 IF (this%vcoord_out(k) >= &
4030 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
4031 this%vcoord_out(k) < &
4032 max(coord_in(kkup), coord_in(kkup+1))) THEN
4042 IF (c_e(kfound)) THEN
4043 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
4044 (coord_in(kfound) - coord_in(kfound-1)))
4046 IF (vartype == var_dir360) THEN
4047 field_out(i,j,k) = &
4048 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
4049 ELSE IF (vartype == var_press) THEN
4050 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
4052 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
4062 DEALLOCATE(coord_in,val_in)
4067 ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN
4069 field_out(:,:,:) = field_in(:,:,:)
4079 FUNCTION interp_var_360(v1, v2, w1, w2)
4080 REAL, INTENT(in) :: v1, v2, w1, w2
4081 REAL :: interp_var_360
4085 IF (abs(v1 - v2) > 180.) THEN
4093 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4095 interp_var_360 = v1*w2 + v2*w1
4098 END FUNCTION interp_var_360
4101 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
4102 REAL, INTENT(in) :: v1(:,:)
4103 REAL, INTENT(in) :: l, h, res
4104 LOGICAL, INTENT(in), OPTIONAL :: mask(:,:)
4112 IF ((h - l) <= res) THEN
4117 IF ( PRESENT(mask)) THEN
4118 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4119 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4121 nl = count(v1 >= l .AND. v1 < m)
4122 nh = count(v1 >= m .AND. v1 < h)
4125 prev = find_prevailing_direction(v1, m, h, res)
4126 ELSE IF (nl > nh) THEN
4127 prev = find_prevailing_direction(v1, l, m, res)
4128 ELSE IF (nl == 0 .AND. nh == 0) THEN
4134 END FUNCTION find_prevailing_direction
4137 END SUBROUTINE grid_transform_compute
4145 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4147 TYPE(grid_transform), INTENT(in) :: this
4148 REAL, INTENT(in) :: field_in(:,:)
4149 REAL, INTENT(out):: field_out(:,:,:)
4150 TYPE(vol7d_var), INTENT(in), OPTIONAL :: var
4151 REAL, INTENT(in), OPTIONAL, TARGET :: coord_3d_in(:,:,:)
4153 real, allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4154 INTEGER :: inn_p, ier, k
4155 INTEGER :: innx, inny, innz, outnx, outny, outnz
4158 call l4f_category_log(this%category,l4f_debug, "start v7d_grid_transform_compute")
4161 field_out(:,:,:) = rmiss
4163 IF (.NOT.this%valid) THEN
4164 call l4f_category_log(this%category,l4f_error, &
4165 "refusing to perform a non valid transformation")
4170 innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
4171 outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
4174 IF (this%trans%trans_type == 'vertint') THEN
4176 IF (innz /= this%innz) THEN
4177 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation")
4178 CALL l4f_category_log(this%category,l4f_error, "inconsistent input shape: "//&
4179 t2c(this%innz)// " /= "//t2c(innz))
4184 IF (outnz /= this%outnz) THEN
4185 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation")
4186 CALL l4f_category_log(this%category,l4f_error, "inconsistent output shape: "//&
4187 t2c(this%outnz)// " /= "//t2c(outnz))
4192 IF (innx /= outnx .OR. inny /= outny) THEN
4193 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation")
4194 CALL l4f_category_log(this%category,l4f_error, "inconsistent hor. sizes: "//&
4195 t2c(innx)// ","//t2c(inny)// " /= "//&
4196 t2c(outnx)// ","//t2c(outny))
4203 IF (innx /= this%innx .OR. inny /= this%inny) THEN
4204 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation")
4205 CALL l4f_category_log(this%category,l4f_error, "inconsistent input shape: "//&
4206 t2c(this%innx)// ","//t2c(this%inny)// " /= "//&
4207 t2c(innx)// ","//t2c(inny))
4212 IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
4213 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation")
4214 CALL l4f_category_log(this%category,l4f_error, "inconsistent output shape: "//&
4215 t2c(this%outnx)// ","//t2c(this%outny)// " /= "//&
4216 t2c(outnx)// ","//t2c(outny))
4221 IF (innz /= outnz) THEN
4222 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation")
4223 CALL l4f_category_log(this%category,l4f_error, "inconsistent vert. sizes: "//&
4224 t2c(innz)// " /= "//t2c(outnz))
4232 call l4f_category_log(this%category,l4f_debug, &
4233 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)// ':'// &
4234 trim(this%trans%sub_type))
4237 IF (this%trans%trans_type == 'inter') THEN
4239 IF (this%trans%sub_type == 'linear') THEN
4241 #ifdef HAVE_LIBNGMATH
4243 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4245 inn_p = count(c_e(field_in(:,k)))
4247 CALL l4f_category_log(this%category,l4f_info, &
4248 "Number of sparse data points: "//t2c(inn_p)// ','//t2c( SIZE(field_in(:,k))))
4252 field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k)))
4253 x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k)))
4254 y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k)))
4256 IF (.NOT.this%trans%extrap) THEN
4257 CALL nnseti( 'ext', 0)
4258 CALL nnsetr( 'nul', rmiss)
4261 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4262 this%outnx, this%outny, real(this%inter_x(:,1)), &
4263 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4266 CALL l4f_category_log(this%category,l4f_error, &
4267 "Error code from NCAR natgrids: "//t2c(ier))
4273 CALL l4f_category_log(this%category,l4f_info, &
4274 "insufficient data in gridded region to triangulate")
4278 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4281 CALL l4f_category_log(this%category,l4f_error, &
4282 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4289 ELSE IF (this%trans%trans_type == 'boxinter' .OR. &
4290 this%trans%trans_type == 'polyinter' .OR. &
4291 this%trans%trans_type == 'vertint' .OR. &
4292 this%trans%trans_type == 'metamorphosis') THEN
4294 CALL compute(this, &
4295 reshape(field_in, (/ SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
4300 END SUBROUTINE grid_transform_v7d_grid_compute
4315 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
4316 REAL, INTENT(in) :: z1,z2,z3,z4
4317 DOUBLE PRECISION, INTENT(in):: x1,y1
4318 DOUBLE PRECISION, INTENT(in):: x3,y3
4319 DOUBLE PRECISION, INTENT(in):: xp,yp
4326 p2 = real((yp-y1)/(y3-y1))
4327 p1 = real((xp-x1)/(x3-x1))
|