|
◆ 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 3105 del file grid_transform_class.F90.
3108 navg = this%trans%box_info%npx*this%trans%box_info%npy
3111 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3112 je = j+this%trans%box_info%npy-1
3115 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3116 ie = i+this%trans%box_info%npx-1
3118 field_out(ii,jj,k) = stat_stddev( &
3119 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3124 ELSE IF (this%trans%sub_type == 'max') THEN
3127 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3128 je = j+this%trans%box_info%npy-1
3131 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3132 ie = i+this%trans%box_info%npx-1
3134 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3136 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3137 mask=(field_in(i:ie,j:je,k) /= rmiss))
3143 ELSE IF (this%trans%sub_type == 'min') THEN
3146 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3147 je = j+this%trans%box_info%npy-1
3150 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3151 ie = i+this%trans%box_info%npx-1
3153 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3155 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3156 mask=(field_in(i:ie,j:je,k) /= rmiss))
3162 ELSE IF (this%trans%sub_type == 'percentile') THEN
3164 navg = this%trans%box_info%npx*this%trans%box_info%npy
3167 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3168 je = j+this%trans%box_info%npy-1
3171 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3172 ie = i+this%trans%box_info%npx-1
3174 field_out(ii:ii,jj,k) = stat_percentile( &
3175 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3176 (/real(this%trans%stat_info%percentile)/))
3181 ELSE IF (this%trans%sub_type == 'frequency') THEN
3185 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3186 je = j+this%trans%box_info%npy-1
3189 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3190 ie = i+this%trans%box_info%npx-1
3192 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3194 field_out(ii,jj,k) = sum(interval_info_valid( &
3195 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3196 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3204 ELSE IF (this%trans%trans_type == 'inter') THEN
3206 IF (this%trans%sub_type == 'near') THEN
3209 DO j = 1, this%outny
3210 DO i = 1, this%outnx
3212 IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3213 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3219 ELSE IF (this%trans%sub_type == 'bilin') THEN
3222 DO j = 1, this%outny
3223 DO i = 1, this%outnx
3225 IF (c_e(this%inter_index_x(i,j))) THEN
3227 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3228 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3229 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3230 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3232 IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN
3234 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3235 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3236 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3237 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3239 xp=this%inter_xp(i,j)
3240 yp=this%inter_yp(i,j)
3242 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3250 ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
3252 DO j = 1, this%outny
3253 DO i = 1, this%outnx
3255 IF (c_e(this%inter_index_x(i,j))) THEN
3257 IF(this%inter_index_x(i,j)-1>0) THEN
3258 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3262 IF(this%inter_index_x(i,j)+1<=this%outnx) THEN
3263 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3267 IF(this%inter_index_y(i,j)+1<=this%outny) THEN
3268 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3272 IF(this%inter_index_y(i,j)-1>0) THEN
3273 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3277 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3286 ELSE IF (this%trans%trans_type == 'intersearch') THEN
3289 likethis%trans%trans_type = 'inter'
3290 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3293 IF ((.NOT.all(c_e(field_out(:,:,k)))) .AND. (any(c_e(field_in(:,:,k))))) THEN
3294 DO j = 1, this%outny
3295 DO i = 1, this%outnx
3296 IF (.NOT.c_e(field_out(i,j,k))) THEN
3300 IF (c_e(field_in(l,m,k))) THEN
3301 disttmp = (this%inter_xp(l,m) - this%inter_x(i,j))**2 + (this%inter_yp(l,m) - this%inter_y(i,j))**2
3302 IF (disttmp < dist) THEN
3304 field_out(i,j,k) = field_in(l,m,k)
3315 ELSE IF (this%trans%trans_type == 'boxinter' &
3316 .OR. this%trans%trans_type == 'polyinter' &
3317 .OR. this%trans%trans_type == 'maskinter') THEN
3319 IF (this%trans%sub_type == 'average') THEN
3321 IF (vartype == var_dir360) THEN
3323 DO j = 1, this%outny
3324 DO i = 1, this%outnx
3325 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3327 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3333 ALLOCATE(nval(this%outnx, this%outny))
3334 field_out(:,:,:) = 0.0
3339 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3340 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3341 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3343 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3344 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3348 WHERE (nval(:,:) /= 0)
3349 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3351 field_out(:,:,k) = rmiss
3357 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3358 this%trans%sub_type == 'stddevnm1') THEN
3360 IF (this%trans%sub_type == 'stddev') THEN
3366 DO j = 1, this%outny
3367 DO i = 1, this%outnx
3369 field_out(i:i,j,k) = stat_stddev( &
3370 reshape(field_in(:,:,k), (/ SIZE(field_in(:,:,k))/)), &
3371 mask=reshape((this%inter_index_x == i .AND. &
3372 this%inter_index_y == j), (/ SIZE(field_in(:,:,k))/)), nm1=nm1)
3377 ELSE IF (this%trans%sub_type == 'max') THEN
3382 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3383 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3384 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3385 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3388 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3397 ELSE IF (this%trans%sub_type == 'min') THEN
3402 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3403 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3404 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3405 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3408 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3416 ELSE IF (this%trans%sub_type == 'percentile') THEN
3419 DO j = 1, this%outny
3420 DO i = 1, this%outnx
3422 field_out(i:i,j,k) = stat_percentile( &
3423 reshape(field_in(:,:,k), (/ SIZE(field_in(:,:,k))/)), &
3424 (/real(this%trans%stat_info%percentile)/), &
3425 mask=reshape((this%inter_index_x == i .AND. &
3426 this%inter_index_y == j), (/ SIZE(field_in(:,:,k))/)))
3431 ELSE IF (this%trans%sub_type == 'frequency') THEN
3433 ALLOCATE(nval(this%outnx, this%outny))
3434 field_out(:,:,:) = 0.0
3439 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3440 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3441 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3442 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3443 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3444 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3448 WHERE (nval(:,:) /= 0)
3449 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3451 field_out(:,:,k) = rmiss
3458 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
3459 np = SIZE(this%stencil,1)/2
3460 ns = SIZE(this%stencil)
3462 IF (this%trans%sub_type == 'average') THEN
3464 IF (vartype == var_dir360) THEN
3466 DO j = 1, this%outny
3467 DO i = 1, this%outnx
3468 IF (c_e(this%inter_index_x(i,j))) THEN
3469 i1 = this%inter_index_x(i,j) - np
3470 i2 = this%inter_index_x(i,j) + np
3471 j1 = this%inter_index_y(i,j) - np
3472 j2 = this%inter_index_y(i,j) + np
3473 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3475 mask=this%stencil(:,:))
3486 DO j = 1, this%outny
3487 DO i = 1, this%outnx
3488 IF (c_e(this%inter_index_x(i,j))) THEN
3489 i1 = this%inter_index_x(i,j) - np
3490 i2 = this%inter_index_x(i,j) + np
3491 j1 = this%inter_index_y(i,j) - np
3492 j2 = this%inter_index_y(i,j) + np
3493 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3495 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3496 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3505 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3506 this%trans%sub_type == 'stddevnm1') THEN
3508 IF (this%trans%sub_type == 'stddev') THEN
3517 DO j = 1, this%outny
3518 DO i = 1, this%outnx
3519 IF (c_e(this%inter_index_x(i,j))) THEN
3520 i1 = this%inter_index_x(i,j) - np
3521 i2 = this%inter_index_x(i,j) + np
3522 j1 = this%inter_index_y(i,j) - np
3523 j2 = this%inter_index_y(i,j) + np
3525 field_out(i:i,j,k) = stat_stddev( &
3526 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3527 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3528 this%stencil(:,:), (/ns/)), nm1=nm1)
3535 ELSE IF (this%trans%sub_type == 'max') THEN
3540 DO j = 1, this%outny
3541 DO i = 1, this%outnx
3542 IF (c_e(this%inter_index_x(i,j))) THEN
3543 i1 = this%inter_index_x(i,j) - np
3544 i2 = this%inter_index_x(i,j) + np
3545 j1 = this%inter_index_y(i,j) - np
3546 j2 = this%inter_index_y(i,j) + np
3547 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3549 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3550 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3558 ELSE IF (this%trans%sub_type == 'min') THEN
3563 DO j = 1, this%outny
3564 DO i = 1, this%outnx
3565 IF (c_e(this%inter_index_x(i,j))) THEN
3566 i1 = this%inter_index_x(i,j) - np
3567 i2 = this%inter_index_x(i,j) + np
3568 j1 = this%inter_index_y(i,j) - np
3569 j2 = this%inter_index_y(i,j) + np
3570 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3572 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3573 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3581 ELSE IF (this%trans%sub_type == 'percentile') THEN
3586 DO j = 1, this%outny
3587 DO i = 1, this%outnx
3588 IF (c_e(this%inter_index_x(i,j))) THEN
3589 i1 = this%inter_index_x(i,j) - np
3590 i2 = this%inter_index_x(i,j) + np
3591 j1 = this%inter_index_y(i,j) - np
3592 j2 = this%inter_index_y(i,j) + np
3594 field_out(i:i,j,k) = stat_percentile( &
3595 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3596 (/real(this%trans%stat_info%percentile)/), &
3597 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3598 this%stencil(:,:), (/ns/)))
3605 ELSE IF (this%trans%sub_type == 'frequency') THEN
3610 DO j = 1, this%outny
3611 DO i = 1, this%outnx
3612 IF (c_e(this%inter_index_x(i,j))) THEN
3613 i1 = this%inter_index_x(i,j) - np
3614 i2 = this%inter_index_x(i,j) + np
3615 j1 = this%inter_index_y(i,j) - np
3616 j2 = this%inter_index_y(i,j) + np
3617 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3619 field_out(i,j,k) = sum(interval_info_valid( &
3620 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3621 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3631 ELSE IF (this%trans%trans_type == 'maskgen') THEN
3634 WHERE(c_e(this%inter_index_x(:,:)))
3635 field_out(:,:,k) = real(this%inter_index_x(:,:))
3639 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
3641 IF (this%trans%sub_type == 'all') THEN
3643 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3645 ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' &
3646 .OR. this%trans%sub_type == 'mask') THEN
3650 field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)))
3653 ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
3654 this%trans%sub_type == 'maskinvalid') THEN
3657 WHERE (this%point_mask(:,:))
3658 field_out(:,:,k) = field_in(:,:,k)
3662 ELSE IF (this%trans%sub_type == 'lemaskinvalid') THEN
3665 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3666 field_out(:,:,k) = field_in(:,:,k)
3668 field_out(:,:,k) = rmiss
3672 ELSE IF (this%trans%sub_type == 'ltmaskinvalid') THEN
3675 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3676 field_out(:,:,k) = field_in(:,:,k)
3678 field_out(:,:,k) = rmiss
3682 ELSE IF (this%trans%sub_type == 'gemaskinvalid') THEN
3685 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3686 field_out(:,:,k) = field_in(:,:,k)
3688 field_out(:,:,k) = rmiss
3692 ELSE IF (this%trans%sub_type == 'gtmaskinvalid') THEN
3695 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3696 field_out(:,:,k) = field_in(:,:,k)
3698 field_out(:,:,k) = rmiss
3702 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
3705 WHERE (c_e(field_in(:,:,k)))
3706 field_out(:,:,k) = field_in(:,:,k)
3708 field_out(:,:,k) = this%val1
3712 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
3713 IF (c_e(this%val1) .AND. c_e(this%val2)) THEN
3714 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3715 .AND. field_in(:,:,:) <= this%val2)
3716 field_out(:,:,:) = rmiss
3718 field_out(:,:,:) = field_in(:,:,:)
3720 ELSE IF (c_e(this%val1)) THEN
3721 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3722 field_out(:,:,:) = rmiss
3724 field_out(:,:,:) = field_in(:,:,:)
3726 ELSE IF (c_e(this%val2)) THEN
3727 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3728 field_out(:,:,:) = rmiss
3730 field_out(:,:,:) = field_in(:,:,:)
3735 ELSE IF (this%trans%trans_type == 'vertint') THEN
3737 IF (this%trans%sub_type == 'linear') THEN
3739 alloc_coord_3d_in_act = .false.
3740 IF ( ASSOCIATED(this%inter_index_z)) THEN
3743 IF (c_e(this%inter_index_z(k))) THEN
3744 z1 = real(this%inter_zp(k))
3745 z2 = real(1.0d0 - this%inter_zp(k))
3746 SELECT CASE(vartype)
3751 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3752 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3753 field_out(i,j,k) = &
3754 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3755 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3763 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3764 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3765 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3766 field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN
3767 field_out(i,j,k) = exp( &
3768 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3769 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3777 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3778 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3779 field_out(i,j,k) = &
3780 field_in(i,j,this%inter_index_z(k))*z2 + &
3781 field_in(i,j,this%inter_index_z(k)+1)*z1
3793 IF ( ALLOCATED(this%coord_3d_in)) THEN
3794 coord_3d_in_act => this%coord_3d_in
3796 CALL l4f_category_log(this%category,l4f_debug, &
3797 "Using external vertical coordinate for vertint")
3800 IF ( PRESENT(coord_3d_in)) THEN
3802 CALL l4f_category_log(this%category,l4f_debug, &
3803 "Using internal vertical coordinate for vertint")
3805 IF (this%dolog) THEN
3806 ALLOCATE(coord_3d_in_act( SIZE(coord_3d_in,1), &
3807 SIZE(coord_3d_in,2), SIZE(coord_3d_in,3)))
3808 alloc_coord_3d_in_act = .true.
3809 WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3810 coord_3d_in_act = log(coord_3d_in)
3812 coord_3d_in_act = rmiss
3815 coord_3d_in_act => coord_3d_in
3818 CALL l4f_category_log(this%category,l4f_error, &
3819 "Missing internal and external vertical coordinate for vertint")
3825 inused = SIZE(coord_3d_in_act,3)
3826 IF (inused < 2) RETURN
3830 IF (.NOT.c_e(this%vcoord_out(k))) cycle
3834 DO kk = 1, max(inused-kkcache-1, kkcache)
3836 kkdown = kkcache - kk + 1
3838 IF (kkdown >= 1) THEN
3839 IF (this%vcoord_out(k) >= &
3840 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3841 this%vcoord_out(k) < &
3842 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1))) THEN
3845 kfound = kkcache + this%levshift
3849 IF (kkup < inused) THEN
3850 IF (this%vcoord_out(k) >= &
3851 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3852 this%vcoord_out(k) < &
3853 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1))) THEN
3856 kfound = kkcache + this%levshift
3863 IF (c_e(kfound)) THEN
3864 IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound+1))) THEN
3865 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3866 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3868 SELECT CASE(vartype)
3871 field_out(i,j,k) = &
3872 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3874 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3875 log(field_in(i,j,kfound+1))*z1)
3878 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3887 IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
3890 ELSE IF (this%trans%sub_type == 'linearsparse') THEN
3893 IF (.NOT. ASSOCIATED(this%vcoord_in) .AND. .NOT. PRESENT(coord_3d_in)) THEN
3894 CALL l4f_category_log(this%category,l4f_error, &
3895 "linearsparse interpolation, no input vert coord available")
3899 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3903 IF ( ASSOCIATED(this%vcoord_in)) THEN
3904 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3905 .AND. c_e(this%vcoord_in(:))
3907 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3908 .AND. c_e(coord_3d_in(i,j,:))
3911 IF (vartype == var_press) THEN
3912 mask_in(:) = mask_in(:) .AND. &
3913 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3915 inused = count(mask_in)
3916 IF (inused > 1) THEN
3917 IF ( ASSOCIATED(this%vcoord_in)) THEN
3918 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3920 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3922 IF (vartype == var_press) THEN
3923 val_in(1:inused) = log(pack( &
3924 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3927 val_in(1:inused) = pack( &
3928 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3935 DO kk = 1, max(inused-kkcache-1, kkcache)
3937 kkdown = kkcache - kk + 1
3939 IF (kkdown >= 1) THEN
3940 IF (this%vcoord_out(k) >= &
3941 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
3942 this%vcoord_out(k) < &
3943 max(coord_in(kkdown), coord_in(kkdown+1))) THEN
3950 IF (kkup < inused) THEN
3951 IF (this%vcoord_out(k) >= &
3952 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
3953 this%vcoord_out(k) < &
3954 max(coord_in(kkup), coord_in(kkup+1))) THEN
3964 IF (c_e(kfound)) THEN
3965 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
3966 (coord_in(kfound) - coord_in(kfound-1)))
3968 IF (vartype == var_dir360) THEN
3969 field_out(i,j,k) = &
3970 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
3971 ELSE IF (vartype == var_press) THEN
3972 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
3974 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
3984 DEALLOCATE(coord_in,val_in)
3989 ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN
3991 field_out(:,:,:) = field_in(:,:,:)
4001 FUNCTION interp_var_360(v1, v2, w1, w2)
4002 REAL, INTENT(in) :: v1, v2, w1, w2
4003 REAL :: interp_var_360
4007 IF (abs(v1 - v2) > 180.) THEN
4015 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4017 interp_var_360 = v1*w2 + v2*w1
4020 END FUNCTION interp_var_360
4023 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
4024 REAL, INTENT(in) :: v1(:,:)
4025 REAL, INTENT(in) :: l, h, res
4026 LOGICAL, INTENT(in), OPTIONAL :: mask(:,:)
4034 IF ((h - l) <= res) THEN
4039 IF ( PRESENT(mask)) THEN
4040 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4041 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4043 nl = count(v1 >= l .AND. v1 < m)
4044 nh = count(v1 >= m .AND. v1 < h)
4047 prev = find_prevailing_direction(v1, m, h, res)
4048 ELSE IF (nl > nh) THEN
4049 prev = find_prevailing_direction(v1, l, m, res)
4050 ELSE IF (nl == 0 .AND. nh == 0) THEN
4056 END FUNCTION find_prevailing_direction
4059 END SUBROUTINE grid_transform_compute
4067 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4069 TYPE(grid_transform), INTENT(in) :: this
4070 REAL, INTENT(in) :: field_in(:,:)
4071 REAL, INTENT(out):: field_out(:,:,:)
4072 TYPE(vol7d_var), INTENT(in), OPTIONAL :: var
4073 REAL, INTENT(in), OPTIONAL, TARGET :: coord_3d_in(:,:,:)
4075 real, allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4076 INTEGER :: inn_p, ier, k
4077 INTEGER :: innx, inny, innz, outnx, outny, outnz
4080 call l4f_category_log(this%category,l4f_debug, "start v7d_grid_transform_compute")
4083 field_out(:,:,:) = rmiss
4085 IF (.NOT.this%valid) THEN
4086 call l4f_category_log(this%category,l4f_error, &
4087 "refusing to perform a non valid transformation")
4092 innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
4093 outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
4096 IF (this%trans%trans_type == 'vertint') THEN
4098 IF (innz /= this%innz) THEN
4099 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation")
4100 CALL l4f_category_log(this%category,l4f_error, "inconsistent input shape: "//&
4101 t2c(this%innz)// " /= "//t2c(innz))
4106 IF (outnz /= this%outnz) THEN
4107 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation")
4108 CALL l4f_category_log(this%category,l4f_error, "inconsistent output shape: "//&
4109 t2c(this%outnz)// " /= "//t2c(outnz))
4114 IF (innx /= outnx .OR. inny /= outny) THEN
4115 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation")
4116 CALL l4f_category_log(this%category,l4f_error, "inconsistent hor. sizes: "//&
4117 t2c(innx)// ","//t2c(inny)// " /= "//&
4118 t2c(outnx)// ","//t2c(outny))
4125 IF (innx /= this%innx .OR. inny /= this%inny) THEN
4126 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation")
4127 CALL l4f_category_log(this%category,l4f_error, "inconsistent input shape: "//&
4128 t2c(this%innx)// ","//t2c(this%inny)// " /= "//&
4129 t2c(innx)// ","//t2c(inny))
4134 IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
4135 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation")
4136 CALL l4f_category_log(this%category,l4f_error, "inconsistent output shape: "//&
4137 t2c(this%outnx)// ","//t2c(this%outny)// " /= "//&
4138 t2c(outnx)// ","//t2c(outny))
4143 IF (innz /= outnz) THEN
4144 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation")
4145 CALL l4f_category_log(this%category,l4f_error, "inconsistent vert. sizes: "//&
4146 t2c(innz)// " /= "//t2c(outnz))
4154 call l4f_category_log(this%category,l4f_debug, &
4155 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)// ':'// &
4156 trim(this%trans%sub_type))
4159 IF (this%trans%trans_type == 'inter') THEN
4161 IF (this%trans%sub_type == 'linear') THEN
4163 #ifdef HAVE_LIBNGMATH
4165 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4167 inn_p = count(c_e(field_in(:,k)))
4169 CALL l4f_category_log(this%category,l4f_info, &
4170 "Number of sparse data points: "//t2c(inn_p)// ','//t2c( SIZE(field_in(:,k))))
4174 field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k)))
4175 x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k)))
4176 y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k)))
4178 IF (.NOT.this%trans%extrap) THEN
4179 CALL nnseti( 'ext', 0)
4180 CALL nnsetr( 'nul', rmiss)
4183 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4184 this%outnx, this%outny, real(this%inter_x(:,1)), &
4185 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4188 CALL l4f_category_log(this%category,l4f_error, &
4189 "Error code from NCAR natgrids: "//t2c(ier))
4195 CALL l4f_category_log(this%category,l4f_info, &
4196 "insufficient data in gridded region to triangulate")
4200 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4203 CALL l4f_category_log(this%category,l4f_error, &
4204 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4211 ELSE IF (this%trans%trans_type == 'boxinter' .OR. &
4212 this%trans%trans_type == 'polyinter' .OR. &
4213 this%trans%trans_type == 'vertint' .OR. &
4214 this%trans%trans_type == 'metamorphosis') THEN
4216 CALL compute(this, &
4217 reshape(field_in, (/ SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
4222 END SUBROUTINE grid_transform_v7d_grid_compute
4237 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
4238 REAL, INTENT(in) :: z1,z2,z3,z4
4239 DOUBLE PRECISION, INTENT(in):: x1,y1
4240 DOUBLE PRECISION, INTENT(in):: x3,y3
4241 DOUBLE PRECISION, INTENT(in):: xp,yp
4248 p2 = real((yp-y1)/(y3-y1))
4249 p1 = real((xp-x1)/(x3-x1))
|