|
◆ 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 3113 del file grid_transform_class.F90.
3117 navg = this%trans%box_info%npx*this%trans%box_info%npy
3120 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3121 je = j+this%trans%box_info%npy-1
3124 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3125 ie = i+this%trans%box_info%npx-1
3127 field_out(ii,jj,k) = stat_stddev( &
3128 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3133 ELSE IF (this%trans%sub_type == 'max') THEN
3136 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3137 je = j+this%trans%box_info%npy-1
3140 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3141 ie = i+this%trans%box_info%npx-1
3143 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3145 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3146 mask=(field_in(i:ie,j:je,k) /= rmiss))
3152 ELSE IF (this%trans%sub_type == 'min') THEN
3155 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3156 je = j+this%trans%box_info%npy-1
3159 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3160 ie = i+this%trans%box_info%npx-1
3162 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3164 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3165 mask=(field_in(i:ie,j:je,k) /= rmiss))
3171 ELSE IF (this%trans%sub_type == 'percentile') THEN
3173 navg = this%trans%box_info%npx*this%trans%box_info%npy
3176 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3177 je = j+this%trans%box_info%npy-1
3180 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3181 ie = i+this%trans%box_info%npx-1
3183 field_out(ii:ii,jj,k) = stat_percentile( &
3184 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3185 (/real(this%trans%stat_info%percentile)/))
3190 ELSE IF (this%trans%sub_type == 'frequency') THEN
3194 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3195 je = j+this%trans%box_info%npy-1
3198 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3199 ie = i+this%trans%box_info%npx-1
3201 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3203 field_out(ii,jj,k) = sum(interval_info_valid( &
3204 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3205 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3213 ELSE IF (this%trans%trans_type == 'inter') THEN
3215 IF (this%trans%sub_type == 'near') THEN
3218 DO j = 1, this%outny
3219 DO i = 1, this%outnx
3221 IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3222 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3228 ELSE IF (this%trans%sub_type == 'bilin') THEN
3231 DO j = 1, this%outny
3232 DO i = 1, this%outnx
3234 IF (c_e(this%inter_index_x(i,j))) THEN
3236 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3237 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3238 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3239 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3241 IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN
3243 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3244 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3245 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3246 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3248 xp=this%inter_xp(i,j)
3249 yp=this%inter_yp(i,j)
3251 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3259 ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
3261 DO j = 1, this%outny
3262 DO i = 1, this%outnx
3264 IF (c_e(this%inter_index_x(i,j))) THEN
3266 IF(this%inter_index_x(i,j)-1>0) THEN
3267 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3271 IF(this%inter_index_x(i,j)+1<=this%outnx) THEN
3272 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3276 IF(this%inter_index_y(i,j)+1<=this%outny) THEN
3277 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3281 IF(this%inter_index_y(i,j)-1>0) THEN
3282 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3286 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3295 ELSE IF (this%trans%trans_type == 'intersearch') THEN
3298 likethis%trans%trans_type = 'inter'
3299 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3300 CALL getenv( 'LIBSIM_DISABLEOPTSEARCH', env_var)
3301 optsearch = len_trim(env_var) == 0
3304 IF ((.NOT.all(c_e(field_out(:,:,k)))) .AND. (any(c_e(field_in(:,:,k))))) THEN
3305 DO j = 1, this%outny
3306 DO i = 1, this%outnx
3307 IF (.NOT.c_e(field_out(i,j,k))) THEN
3311 ix = this%inter_index_x(i,j)
3312 iy = this%inter_index_y(i,j)
3313 DO s = 0, max(this%innx, this%inny)
3315 DO m = iy-s, iy+s, max(2*s, 1)
3316 IF (m < 1 .OR. m > this%inny) cycle
3317 DO l = max(1, ix-s), min(this%innx, ix+s)
3318 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3319 IF (c_e(field_in(l,m,k))) THEN
3320 IF (disttmp < dist) THEN
3322 field_out(i,j,k) = field_in(l,m,k)
3324 ELSE IF (disttmp == dist) THEN
3325 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3326 nearcount = nearcount + 1
3329 IF (disttmp < dist) farenough = .false.
3332 DO m = max(1, iy-s+1), min(this%inny, iy+s-1)
3333 DO l = ix-s, ix+s, 2*s
3334 IF (l < 1 .OR. l > this%innx) cycle
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 IF (s > 0 .AND. farenough) EXIT
3354 IF (c_e(field_in(l,m,k))) THEN
3355 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3356 IF (disttmp < dist) THEN
3358 field_out(i,j,k) = field_in(l,m,k)
3360 ELSE IF (disttmp == dist) THEN
3361 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3362 nearcount = nearcount + 1
3369 IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount
3376 ELSE IF (this%trans%trans_type == 'boxinter' &
3377 .OR. this%trans%trans_type == 'polyinter' &
3378 .OR. this%trans%trans_type == 'maskinter') THEN
3380 IF (this%trans%sub_type == 'average') THEN
3382 IF (vartype == var_dir360) THEN
3384 DO j = 1, this%outny
3385 DO i = 1, this%outnx
3386 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3388 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3394 ALLOCATE(nval(this%outnx, this%outny))
3395 field_out(:,:,:) = 0.0
3400 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3401 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3402 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3404 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3405 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3409 WHERE (nval(:,:) /= 0)
3410 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3412 field_out(:,:,k) = rmiss
3418 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3419 this%trans%sub_type == 'stddevnm1') THEN
3421 IF (this%trans%sub_type == 'stddev') THEN
3427 DO j = 1, this%outny
3428 DO i = 1, this%outnx
3430 field_out(i:i,j,k) = stat_stddev( &
3431 reshape(field_in(:,:,k), (/ SIZE(field_in(:,:,k))/)), &
3432 mask=reshape((this%inter_index_x == i .AND. &
3433 this%inter_index_y == j), (/ SIZE(field_in(:,:,k))/)), nm1=nm1)
3438 ELSE IF (this%trans%sub_type == 'max') THEN
3443 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3444 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3445 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3446 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3449 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3458 ELSE IF (this%trans%sub_type == 'min') THEN
3463 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3464 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3465 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3466 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3469 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3477 ELSE IF (this%trans%sub_type == 'percentile') THEN
3480 DO j = 1, this%outny
3481 DO i = 1, this%outnx
3483 field_out(i:i,j,k) = stat_percentile( &
3484 reshape(field_in(:,:,k), (/ SIZE(field_in(:,:,k))/)), &
3485 (/real(this%trans%stat_info%percentile)/), &
3486 mask=reshape((this%inter_index_x == i .AND. &
3487 this%inter_index_y == j), (/ SIZE(field_in(:,:,k))/)))
3492 ELSE IF (this%trans%sub_type == 'frequency') THEN
3494 ALLOCATE(nval(this%outnx, this%outny))
3495 field_out(:,:,:) = 0.0
3500 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3501 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3502 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3503 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3504 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3505 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3509 WHERE (nval(:,:) /= 0)
3510 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3512 field_out(:,:,k) = rmiss
3519 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
3520 np = SIZE(this%stencil,1)/2
3521 ns = SIZE(this%stencil)
3523 IF (this%trans%sub_type == 'average') THEN
3525 IF (vartype == var_dir360) THEN
3527 DO j = 1, this%outny
3528 DO i = 1, this%outnx
3529 IF (c_e(this%inter_index_x(i,j))) THEN
3530 i1 = this%inter_index_x(i,j) - np
3531 i2 = this%inter_index_x(i,j) + np
3532 j1 = this%inter_index_y(i,j) - np
3533 j2 = this%inter_index_y(i,j) + np
3534 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3536 mask=this%stencil(:,:))
3547 DO j = 1, this%outny
3548 DO i = 1, this%outnx
3549 IF (c_e(this%inter_index_x(i,j))) THEN
3550 i1 = this%inter_index_x(i,j) - np
3551 i2 = this%inter_index_x(i,j) + np
3552 j1 = this%inter_index_y(i,j) - np
3553 j2 = this%inter_index_y(i,j) + np
3554 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3556 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3557 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3566 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3567 this%trans%sub_type == 'stddevnm1') THEN
3569 IF (this%trans%sub_type == 'stddev') THEN
3578 DO j = 1, this%outny
3579 DO i = 1, this%outnx
3580 IF (c_e(this%inter_index_x(i,j))) THEN
3581 i1 = this%inter_index_x(i,j) - np
3582 i2 = this%inter_index_x(i,j) + np
3583 j1 = this%inter_index_y(i,j) - np
3584 j2 = this%inter_index_y(i,j) + np
3586 field_out(i:i,j,k) = stat_stddev( &
3587 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3588 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3589 this%stencil(:,:), (/ns/)), nm1=nm1)
3596 ELSE IF (this%trans%sub_type == 'max') THEN
3601 DO j = 1, this%outny
3602 DO i = 1, this%outnx
3603 IF (c_e(this%inter_index_x(i,j))) THEN
3604 i1 = this%inter_index_x(i,j) - np
3605 i2 = this%inter_index_x(i,j) + np
3606 j1 = this%inter_index_y(i,j) - np
3607 j2 = this%inter_index_y(i,j) + np
3608 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3610 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3611 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3619 ELSE IF (this%trans%sub_type == 'min') THEN
3624 DO j = 1, this%outny
3625 DO i = 1, this%outnx
3626 IF (c_e(this%inter_index_x(i,j))) THEN
3627 i1 = this%inter_index_x(i,j) - np
3628 i2 = this%inter_index_x(i,j) + np
3629 j1 = this%inter_index_y(i,j) - np
3630 j2 = this%inter_index_y(i,j) + np
3631 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3633 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3634 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3642 ELSE IF (this%trans%sub_type == 'percentile') THEN
3647 DO j = 1, this%outny
3648 DO i = 1, this%outnx
3649 IF (c_e(this%inter_index_x(i,j))) THEN
3650 i1 = this%inter_index_x(i,j) - np
3651 i2 = this%inter_index_x(i,j) + np
3652 j1 = this%inter_index_y(i,j) - np
3653 j2 = this%inter_index_y(i,j) + np
3655 field_out(i:i,j,k) = stat_percentile( &
3656 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3657 (/real(this%trans%stat_info%percentile)/), &
3658 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3659 this%stencil(:,:), (/ns/)))
3666 ELSE IF (this%trans%sub_type == 'frequency') THEN
3671 DO j = 1, this%outny
3672 DO i = 1, this%outnx
3673 IF (c_e(this%inter_index_x(i,j))) THEN
3674 i1 = this%inter_index_x(i,j) - np
3675 i2 = this%inter_index_x(i,j) + np
3676 j1 = this%inter_index_y(i,j) - np
3677 j2 = this%inter_index_y(i,j) + np
3678 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3680 field_out(i,j,k) = sum(interval_info_valid( &
3681 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3682 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3692 ELSE IF (this%trans%trans_type == 'maskgen') THEN
3695 WHERE(c_e(this%inter_index_x(:,:)))
3696 field_out(:,:,k) = real(this%inter_index_x(:,:))
3700 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
3702 IF (this%trans%sub_type == 'all') THEN
3704 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3706 ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' &
3707 .OR. this%trans%sub_type == 'mask') THEN
3711 field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)))
3714 ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
3715 this%trans%sub_type == 'maskinvalid') THEN
3718 WHERE (this%point_mask(:,:))
3719 field_out(:,:,k) = field_in(:,:,k)
3723 ELSE IF (this%trans%sub_type == 'lemaskinvalid') THEN
3726 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3727 field_out(:,:,k) = field_in(:,:,k)
3729 field_out(:,:,k) = rmiss
3733 ELSE IF (this%trans%sub_type == 'ltmaskinvalid') THEN
3736 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3737 field_out(:,:,k) = field_in(:,:,k)
3739 field_out(:,:,k) = rmiss
3743 ELSE IF (this%trans%sub_type == 'gemaskinvalid') THEN
3746 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3747 field_out(:,:,k) = field_in(:,:,k)
3749 field_out(:,:,k) = rmiss
3753 ELSE IF (this%trans%sub_type == 'gtmaskinvalid') THEN
3756 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3757 field_out(:,:,k) = field_in(:,:,k)
3759 field_out(:,:,k) = rmiss
3763 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
3766 WHERE (c_e(field_in(:,:,k)))
3767 field_out(:,:,k) = field_in(:,:,k)
3769 field_out(:,:,k) = this%val1
3773 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
3774 IF (c_e(this%val1) .AND. c_e(this%val2)) THEN
3775 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3776 .AND. field_in(:,:,:) <= this%val2)
3777 field_out(:,:,:) = rmiss
3779 field_out(:,:,:) = field_in(:,:,:)
3781 ELSE IF (c_e(this%val1)) THEN
3782 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3783 field_out(:,:,:) = rmiss
3785 field_out(:,:,:) = field_in(:,:,:)
3787 ELSE IF (c_e(this%val2)) THEN
3788 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3789 field_out(:,:,:) = rmiss
3791 field_out(:,:,:) = field_in(:,:,:)
3796 ELSE IF (this%trans%trans_type == 'vertint') THEN
3798 IF (this%trans%sub_type == 'linear') THEN
3800 alloc_coord_3d_in_act = .false.
3801 IF ( ASSOCIATED(this%inter_index_z)) THEN
3804 IF (c_e(this%inter_index_z(k))) THEN
3805 z1 = real(this%inter_zp(k))
3806 z2 = real(1.0d0 - this%inter_zp(k))
3807 SELECT CASE(vartype)
3812 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3813 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3814 field_out(i,j,k) = &
3815 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3816 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3824 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3825 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3826 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3827 field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN
3828 field_out(i,j,k) = exp( &
3829 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3830 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3838 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3839 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3840 field_out(i,j,k) = &
3841 field_in(i,j,this%inter_index_z(k))*z2 + &
3842 field_in(i,j,this%inter_index_z(k)+1)*z1
3854 IF ( ALLOCATED(this%coord_3d_in)) THEN
3855 coord_3d_in_act => this%coord_3d_in
3857 CALL l4f_category_log(this%category,l4f_debug, &
3858 "Using external vertical coordinate for vertint")
3861 IF ( PRESENT(coord_3d_in)) THEN
3863 CALL l4f_category_log(this%category,l4f_debug, &
3864 "Using internal vertical coordinate for vertint")
3866 IF (this%dolog) THEN
3867 ALLOCATE(coord_3d_in_act( SIZE(coord_3d_in,1), &
3868 SIZE(coord_3d_in,2), SIZE(coord_3d_in,3)))
3869 alloc_coord_3d_in_act = .true.
3870 WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3871 coord_3d_in_act = log(coord_3d_in)
3873 coord_3d_in_act = rmiss
3876 coord_3d_in_act => coord_3d_in
3879 CALL l4f_category_log(this%category,l4f_error, &
3880 "Missing internal and external vertical coordinate for vertint")
3886 inused = SIZE(coord_3d_in_act,3)
3887 IF (inused < 2) RETURN
3891 IF (.NOT.c_e(this%vcoord_out(k))) cycle
3895 DO kk = 1, max(inused-kkcache-1, kkcache)
3897 kkdown = kkcache - kk + 1
3899 IF (kkdown >= 1) THEN
3900 IF (this%vcoord_out(k) >= &
3901 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3902 this%vcoord_out(k) < &
3903 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1))) THEN
3906 kfound = kkcache + this%levshift
3910 IF (kkup < inused) THEN
3911 IF (this%vcoord_out(k) >= &
3912 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3913 this%vcoord_out(k) < &
3914 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1))) THEN
3917 kfound = kkcache + this%levshift
3924 IF (c_e(kfound)) THEN
3925 IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound+1))) THEN
3926 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3927 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3929 SELECT CASE(vartype)
3932 field_out(i,j,k) = &
3933 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3935 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3936 log(field_in(i,j,kfound+1))*z1)
3939 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3948 IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
3951 ELSE IF (this%trans%sub_type == 'linearsparse') THEN
3954 IF (.NOT. ASSOCIATED(this%vcoord_in) .AND. .NOT. PRESENT(coord_3d_in)) THEN
3955 CALL l4f_category_log(this%category,l4f_error, &
3956 "linearsparse interpolation, no input vert coord available")
3960 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3964 IF ( ASSOCIATED(this%vcoord_in)) THEN
3965 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3966 .AND. c_e(this%vcoord_in(:))
3968 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3969 .AND. c_e(coord_3d_in(i,j,:))
3972 IF (vartype == var_press) THEN
3973 mask_in(:) = mask_in(:) .AND. &
3974 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3976 inused = count(mask_in)
3977 IF (inused > 1) THEN
3978 IF ( ASSOCIATED(this%vcoord_in)) THEN
3979 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3981 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3983 IF (vartype == var_press) THEN
3984 val_in(1:inused) = log(pack( &
3985 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3988 val_in(1:inused) = pack( &
3989 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3996 DO kk = 1, max(inused-kkcache-1, kkcache)
3998 kkdown = kkcache - kk + 1
4000 IF (kkdown >= 1) THEN
4001 IF (this%vcoord_out(k) >= &
4002 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
4003 this%vcoord_out(k) < &
4004 max(coord_in(kkdown), coord_in(kkdown+1))) THEN
4011 IF (kkup < inused) THEN
4012 IF (this%vcoord_out(k) >= &
4013 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
4014 this%vcoord_out(k) < &
4015 max(coord_in(kkup), coord_in(kkup+1))) THEN
4025 IF (c_e(kfound)) THEN
4026 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
4027 (coord_in(kfound) - coord_in(kfound-1)))
4029 IF (vartype == var_dir360) THEN
4030 field_out(i,j,k) = &
4031 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
4032 ELSE IF (vartype == var_press) THEN
4033 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
4035 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
4045 DEALLOCATE(coord_in,val_in)
4050 ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN
4052 field_out(:,:,:) = field_in(:,:,:)
4062 FUNCTION interp_var_360(v1, v2, w1, w2)
4063 REAL, INTENT(in) :: v1, v2, w1, w2
4064 REAL :: interp_var_360
4068 IF (abs(v1 - v2) > 180.) THEN
4076 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4078 interp_var_360 = v1*w2 + v2*w1
4081 END FUNCTION interp_var_360
4084 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
4085 REAL, INTENT(in) :: v1(:,:)
4086 REAL, INTENT(in) :: l, h, res
4087 LOGICAL, INTENT(in), OPTIONAL :: mask(:,:)
4095 IF ((h - l) <= res) THEN
4100 IF ( PRESENT(mask)) THEN
4101 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4102 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4104 nl = count(v1 >= l .AND. v1 < m)
4105 nh = count(v1 >= m .AND. v1 < h)
4108 prev = find_prevailing_direction(v1, m, h, res)
4109 ELSE IF (nl > nh) THEN
4110 prev = find_prevailing_direction(v1, l, m, res)
4111 ELSE IF (nl == 0 .AND. nh == 0) THEN
4117 END FUNCTION find_prevailing_direction
4120 END SUBROUTINE grid_transform_compute
4128 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4130 TYPE(grid_transform), INTENT(in) :: this
4131 REAL, INTENT(in) :: field_in(:,:)
4132 REAL, INTENT(out):: field_out(:,:,:)
4133 TYPE(vol7d_var), INTENT(in), OPTIONAL :: var
4134 REAL, INTENT(in), OPTIONAL, TARGET :: coord_3d_in(:,:,:)
4136 real, allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4137 INTEGER :: inn_p, ier, k
4138 INTEGER :: innx, inny, innz, outnx, outny, outnz
4141 call l4f_category_log(this%category,l4f_debug, "start v7d_grid_transform_compute")
4144 field_out(:,:,:) = rmiss
4146 IF (.NOT.this%valid) THEN
4147 call l4f_category_log(this%category,l4f_error, &
4148 "refusing to perform a non valid transformation")
4153 innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
4154 outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
4157 IF (this%trans%trans_type == 'vertint') THEN
4159 IF (innz /= this%innz) THEN
4160 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation")
4161 CALL l4f_category_log(this%category,l4f_error, "inconsistent input shape: "//&
4162 t2c(this%innz)// " /= "//t2c(innz))
4167 IF (outnz /= this%outnz) THEN
4168 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation")
4169 CALL l4f_category_log(this%category,l4f_error, "inconsistent output shape: "//&
4170 t2c(this%outnz)// " /= "//t2c(outnz))
4175 IF (innx /= outnx .OR. inny /= outny) THEN
4176 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation")
4177 CALL l4f_category_log(this%category,l4f_error, "inconsistent hor. sizes: "//&
4178 t2c(innx)// ","//t2c(inny)// " /= "//&
4179 t2c(outnx)// ","//t2c(outny))
4186 IF (innx /= this%innx .OR. inny /= this%inny) THEN
4187 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation")
4188 CALL l4f_category_log(this%category,l4f_error, "inconsistent input shape: "//&
4189 t2c(this%innx)// ","//t2c(this%inny)// " /= "//&
4190 t2c(innx)// ","//t2c(inny))
4195 IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
4196 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation")
4197 CALL l4f_category_log(this%category,l4f_error, "inconsistent output shape: "//&
4198 t2c(this%outnx)// ","//t2c(this%outny)// " /= "//&
4199 t2c(outnx)// ","//t2c(outny))
4204 IF (innz /= outnz) THEN
4205 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation")
4206 CALL l4f_category_log(this%category,l4f_error, "inconsistent vert. sizes: "//&
4207 t2c(innz)// " /= "//t2c(outnz))
4215 call l4f_category_log(this%category,l4f_debug, &
4216 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)// ':'// &
4217 trim(this%trans%sub_type))
4220 IF (this%trans%trans_type == 'inter') THEN
4222 IF (this%trans%sub_type == 'linear') THEN
4224 #ifdef HAVE_LIBNGMATH
4226 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4228 inn_p = count(c_e(field_in(:,k)))
4230 CALL l4f_category_log(this%category,l4f_info, &
4231 "Number of sparse data points: "//t2c(inn_p)// ','//t2c( SIZE(field_in(:,k))))
4235 field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k)))
4236 x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k)))
4237 y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k)))
4239 IF (.NOT.this%trans%extrap) THEN
4240 CALL nnseti( 'ext', 0)
4241 CALL nnsetr( 'nul', rmiss)
4244 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4245 this%outnx, this%outny, real(this%inter_x(:,1)), &
4246 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4249 CALL l4f_category_log(this%category,l4f_error, &
4250 "Error code from NCAR natgrids: "//t2c(ier))
4256 CALL l4f_category_log(this%category,l4f_info, &
4257 "insufficient data in gridded region to triangulate")
4261 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4264 CALL l4f_category_log(this%category,l4f_error, &
4265 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4272 ELSE IF (this%trans%trans_type == 'boxinter' .OR. &
4273 this%trans%trans_type == 'polyinter' .OR. &
4274 this%trans%trans_type == 'vertint' .OR. &
4275 this%trans%trans_type == 'metamorphosis') THEN
4277 CALL compute(this, &
4278 reshape(field_in, (/ SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
4283 END SUBROUTINE grid_transform_v7d_grid_compute
4298 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
4299 REAL, INTENT(in) :: z1,z2,z3,z4
4300 DOUBLE PRECISION, INTENT(in):: x1,y1
4301 DOUBLE PRECISION, INTENT(in):: x3,y3
4302 DOUBLE PRECISION, INTENT(in):: xp,yp
4309 p2 = real((yp-y1)/(y3-y1))
4310 p1 = real((xp-x1)/(x3-x1))
|