|
◆ 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 3132 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))
|