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