|
◆ 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 3067 del file grid_transform_class.F90.
3074 navg = this%trans%box_info%npx*this%trans%box_info%npy
3077 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3078 je = j+this%trans%box_info%npy-1
3081 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3082 ie = i+this%trans%box_info%npx-1
3084 field_out(ii,jj,k) = stat_stddev( &
3085 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3090 ELSE IF (this%trans%sub_type == 'max') THEN
3093 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3094 je = j+this%trans%box_info%npy-1
3097 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3098 ie = i+this%trans%box_info%npx-1
3100 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3102 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3103 mask=(field_in(i:ie,j:je,k) /= rmiss))
3109 ELSE IF (this%trans%sub_type == 'min') THEN
3112 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3113 je = j+this%trans%box_info%npy-1
3116 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3117 ie = i+this%trans%box_info%npx-1
3119 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3121 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3122 mask=(field_in(i:ie,j:je,k) /= rmiss))
3128 ELSE IF (this%trans%sub_type == 'percentile') THEN
3130 navg = this%trans%box_info%npx*this%trans%box_info%npy
3133 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3134 je = j+this%trans%box_info%npy-1
3137 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3138 ie = i+this%trans%box_info%npx-1
3140 field_out(ii:ii,jj,k) = stat_percentile( &
3141 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3142 (/real(this%trans%stat_info%percentile)/))
3147 ELSE IF (this%trans%sub_type == 'frequency') THEN
3151 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3152 je = j+this%trans%box_info%npy-1
3155 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3156 ie = i+this%trans%box_info%npx-1
3158 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3160 field_out(ii,jj,k) = sum(interval_info_valid( &
3161 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3162 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3170 ELSE IF (this%trans%trans_type == 'inter') THEN
3172 IF (this%trans%sub_type == 'near') THEN
3175 DO j = 1, this%outny
3176 DO i = 1, this%outnx
3178 IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3179 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3185 ELSE IF (this%trans%sub_type == 'bilin') THEN
3188 DO j = 1, this%outny
3189 DO i = 1, this%outnx
3191 IF (c_e(this%inter_index_x(i,j))) THEN
3193 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3194 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j
3195 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j
3196 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+
3198 IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN
3200 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y
3201 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y
3202 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y
3203 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y
3205 xp=this%inter_xp(i,j)
3206 yp=this%inter_yp(i,j)
3208 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3216 ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
3218 DO j = 1, this%outny
3219 DO i = 1, this%outnx
3221 IF (c_e(this%inter_index_x(i,j))) THEN
3223 IF(this%inter_index_x(i,j)-1>0) THEN
3224 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y
3228 IF(this%inter_index_x(i,j)+1<=this%outnx) THEN
3229 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y
3233 IF(this%inter_index_y(i,j)+1<=this%outny) THEN
3234 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y
3238 IF(this%inter_index_y(i,j)-1>0) THEN
3239 z(4)=field_in(this%inter_index_x(i,j) ,this%inter_index_y
3243 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j
3252 ELSE IF (this%trans%trans_type == 'boxinter' &
3253 .OR. this%trans%trans_type == 'polyinter' &
3254 .OR. this%trans%trans_type == 'maskinter') THEN
3256 IF (this%trans%sub_type == 'average') THEN
3258 IF (vartype == var_dir360) THEN
3260 DO j = 1, this%outny
3261 DO i = 1, this%outnx
3262 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k)
3264 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(
3270 ALLOCATE(nval(this%outnx, this%outny))
3271 field_out(:,:,:) = 0.0
3276 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3277 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3278 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j)
3280 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3281 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3285 WHERE (nval(:,:) /= 0)
3286 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3288 field_out(:,:,k) = rmiss
3294 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3295 this%trans%sub_type == 'stddevnm1') THEN
3297 IF (this%trans%sub_type == 'stddev') THEN
3303 DO j = 1, this%outny
3304 DO i = 1, this%outnx
3306 field_out(i:i,j,k) = stat_stddev( &
3307 reshape(field_in(:,:,k), (/ SIZE(field_in(:,:,k))/)), &
3308 mask=reshape((this%inter_index_x == i .AND. &
3309 this%inter_index_y == j), (/ SIZE(field_in(:,:,k))/)), nm1=nm1
3314 ELSE IF (this%trans%sub_type == 'max') THEN
3319 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3320 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y THEN
3321 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3322 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i
3325 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3334 ELSE IF (this%trans%sub_type == 'min') THEN
3339 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3340 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y THEN
3341 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3342 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i
3345 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3353 ELSE IF (this%trans%sub_type == 'percentile') THEN
3356 DO j = 1, this%outny
3357 DO i = 1, this%outnx
3359 field_out(i:i,j,k) = stat_percentile( &
3360 reshape(field_in(:,:,k), (/ SIZE(field_in(:,:,k))/)), &
3361 (/real(this%trans%stat_info%percentile)/), &
3362 mask=reshape((this%inter_index_x == i .AND. &
3363 this%inter_index_y == j), (/ SIZE(field_in(:,:,k))/)))
3368 ELSE IF (this%trans%sub_type == 'frequency') THEN
3370 ALLOCATE(nval(this%outnx, this%outny))
3371 field_out(:,:,:) = 0.0
3376 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3377 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3378 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k
3379 interval_info_valid(this%trans%interval_info, field_in(i,j,k
3380 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3381 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3385 WHERE (nval(:,:) /= 0)
3386 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3388 field_out(:,:,k) = rmiss
3395 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
3396 np = SIZE(this%stencil,1)/2
3397 ns = SIZE(this%stencil)
3399 IF (this%trans%sub_type == 'average') THEN
3401 IF (vartype == var_dir360) THEN
3403 DO j = 1, this%outny
3404 DO i = 1, this%outnx
3405 IF (c_e(this%inter_index_x(i,j))) THEN
3406 i1 = this%inter_index_x(i,j) - np
3407 i2 = this%inter_index_x(i,j) + np
3408 j1 = this%inter_index_y(i,j) - np
3409 j2 = this%inter_index_y(i,j) + np
3410 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2
3412 mask=this%stencil(:,:))
3423 DO j = 1, this%outny
3424 DO i = 1, this%outnx
3425 IF (c_e(this%inter_index_x(i,j))) THEN
3426 i1 = this%inter_index_x(i,j) - np
3427 i2 = this%inter_index_x(i,j) + np
3428 j1 = this%inter_index_y(i,j) - np
3429 j2 = this%inter_index_y(i,j) + np
3430 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil
3432 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3433 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil
3442 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3443 this%trans%sub_type == 'stddevnm1') THEN
3445 IF (this%trans%sub_type == 'stddev') THEN
3454 DO j = 1, this%outny
3455 DO i = 1, this%outnx
3456 IF (c_e(this%inter_index_x(i,j))) THEN
3457 i1 = this%inter_index_x(i,j) - np
3458 i2 = this%inter_index_x(i,j) + np
3459 j1 = this%inter_index_y(i,j) - np
3460 j2 = this%inter_index_y(i,j) + np
3462 field_out(i:i,j,k) = stat_stddev( &
3463 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3464 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3465 this%stencil(:,:), (/ns/)), nm1=nm1)
3472 ELSE IF (this%trans%sub_type == 'max') THEN
3477 DO j = 1, this%outny
3478 DO i = 1, this%outnx
3479 IF (c_e(this%inter_index_x(i,j))) THEN
3480 i1 = this%inter_index_x(i,j) - np
3481 i2 = this%inter_index_x(i,j) + np
3482 j1 = this%inter_index_y(i,j) - np
3483 j2 = this%inter_index_y(i,j) + np
3484 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil
3486 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3487 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(
3495 ELSE IF (this%trans%sub_type == 'min') THEN
3500 DO j = 1, this%outny
3501 DO i = 1, this%outnx
3502 IF (c_e(this%inter_index_x(i,j))) THEN
3503 i1 = this%inter_index_x(i,j) - np
3504 i2 = this%inter_index_x(i,j) + np
3505 j1 = this%inter_index_y(i,j) - np
3506 j2 = this%inter_index_y(i,j) + np
3507 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil
3509 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3510 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(
3518 ELSE IF (this%trans%sub_type == 'percentile') THEN
3523 DO j = 1, this%outny
3524 DO i = 1, this%outnx
3525 IF (c_e(this%inter_index_x(i,j))) THEN
3526 i1 = this%inter_index_x(i,j) - np
3527 i2 = this%inter_index_x(i,j) + np
3528 j1 = this%inter_index_y(i,j) - np
3529 j2 = this%inter_index_y(i,j) + np
3531 field_out(i:i,j,k) = stat_percentile( &
3532 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3533 (/real(this%trans%stat_info%percentile)/), &
3534 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3535 this%stencil(:,:), (/ns/)))
3542 ELSE IF (this%trans%sub_type == 'frequency') THEN
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(interval_info_valid( &
3557 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3558 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(
3568 ELSE IF (this%trans%trans_type == 'maskgen') THEN
3571 WHERE(c_e(this%inter_index_x(:,:)))
3572 field_out(:,:,k) = real(this%inter_index_x(:,:))
3576 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
3578 IF (this%trans%sub_type == 'all') THEN
3580 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz
3582 ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly'
3583 .OR. this%trans%sub_type == 'mask') THEN
3587 field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)
3590 ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
3591 this%trans%sub_type == 'maskinvalid') THEN
3594 WHERE (this%point_mask(:,:))
3595 field_out(:,:,k) = field_in(:,:,k)
3599 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
3602 WHERE (c_e(field_in(:,:,k)))
3603 field_out(:,:,k) = field_in(:,:,k)
3605 field_out(:,:,k) = this%val1
3609 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
3610 IF (c_e(this%val1) .AND. c_e(this%val2)) THEN
3611 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3612 .AND. field_in(:,:,:) <= this%val2)
3613 field_out(:,:,:) = rmiss
3615 field_out(:,:,:) = field_in(:,:,:)
3617 ELSE IF (c_e(this%val1)) THEN
3618 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3619 field_out(:,:,:) = rmiss
3621 field_out(:,:,:) = field_in(:,:,:)
3623 ELSE IF (c_e(this%val2)) THEN
3624 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3625 field_out(:,:,:) = rmiss
3627 field_out(:,:,:) = field_in(:,:,:)
3632 ELSE IF (this%trans%trans_type == 'vertint') THEN
3634 IF (this%trans%sub_type == 'linear') THEN
3636 alloc_coord_3d_in_act = .false.
3637 IF ( ASSOCIATED(this%inter_index_z)) THEN
3640 IF (c_e(this%inter_index_z(k))) THEN
3641 z1 = real(this%inter_zp(k))
3642 z2 = real(1.0d0 - this%inter_zp(k))
3643 SELECT CASE(vartype)
3648 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3649 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3650 field_out(i,j,k) = &
3651 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3652 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3660 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3661 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3662 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3663 field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN
3664 field_out(i,j,k) = exp( &
3665 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3666 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3674 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3675 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3676 field_out(i,j,k) = &
3677 field_in(i,j,this%inter_index_z(k))*z2 + &
3678 field_in(i,j,this%inter_index_z(k)+1)*z1
3690 IF ( ALLOCATED(this%coord_3d_in)) THEN
3691 coord_3d_in_act => this%coord_3d_in
3693 CALL l4f_category_log(this%category,l4f_debug, &
3694 "Using external vertical coordinate for vertint")
3697 IF ( PRESENT(coord_3d_in)) THEN
3699 CALL l4f_category_log(this%category,l4f_debug, &
3700 "Using internal vertical coordinate for vertint")
3702 IF (this%dolog) THEN
3703 ALLOCATE(coord_3d_in_act( SIZE(coord_3d_in,1), &
3704 SIZE(coord_3d_in,2), SIZE(coord_3d_in,3)))
3705 alloc_coord_3d_in_act = .true.
3706 WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3707 coord_3d_in_act = log(coord_3d_in)
3709 coord_3d_in_act = rmiss
3712 coord_3d_in_act => coord_3d_in
3715 CALL l4f_category_log(this%category,l4f_error, &
3716 "Missing internal and external vertical coordinate for vertint"
3722 inused = SIZE(coord_3d_in_act,3)
3723 IF (inused < 2) RETURN
3727 IF (.NOT.c_e(this%vcoord_out(k))) cycle
3731 DO kk = 1, max(inused-kkcache-1, kkcache)
3733 kkdown = kkcache - kk + 1
3735 IF (kkdown >= 1) THEN
3736 IF (this%vcoord_out(k) >= &
3737 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown
3738 this%vcoord_out(k) < &
3739 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown THEN
3742 kfound = kkcache + this%levshift
3746 IF (kkup < inused) THEN
3747 IF (this%vcoord_out(k) >= &
3748 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup
3749 this%vcoord_out(k) < &
3750 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup THEN
3753 kfound = kkcache + this%levshift
3760 IF (c_e(kfound)) THEN
3761 IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound THEN
3762 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin
3763 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin
3765 SELECT CASE(vartype)
3768 field_out(i,j,k) = &
3769 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound
3771 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 +
3772 log(field_in(i,j,kfound+1))*z1)
3775 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i
3784 IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
3787 ELSE IF (this%trans%sub_type == 'linearsparse') THEN
3790 IF (.NOT. ASSOCIATED(this%vcoord_in) .AND. .NOT. PRESENT(coord_3d_in)) THEN
3791 CALL l4f_category_log(this%category,l4f_error, &
3792 "linearsparse interpolation, no input vert coord available")
3796 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3800 IF ( ASSOCIATED(this%vcoord_in)) THEN
3801 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused
3802 .AND. c_e(this%vcoord_in(:))
3804 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused
3805 .AND. c_e(coord_3d_in(i,j,:))
3808 IF (vartype == var_press) THEN
3809 mask_in(:) = mask_in(:) .AND. &
3810 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0
3812 inused = count(mask_in)
3813 IF (inused > 1) THEN
3814 IF ( ASSOCIATED(this%vcoord_in)) THEN
3815 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3817 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3819 IF (vartype == var_press) THEN
3820 val_in(1:inused) = log(pack( &
3821 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3824 val_in(1:inused) = pack( &
3825 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3832 DO kk = 1, max(inused-kkcache-1, kkcache)
3834 kkdown = kkcache - kk + 1
3836 IF (kkdown >= 1) THEN
3837 IF (this%vcoord_out(k) >= &
3838 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
3839 this%vcoord_out(k) < &
3840 max(coord_in(kkdown), coord_in(kkdown+1))) THEN
3847 IF (kkup < inused) THEN
3848 IF (this%vcoord_out(k) >= &
3849 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
3850 this%vcoord_out(k) < &
3851 max(coord_in(kkup), coord_in(kkup+1))) THEN
3861 IF (c_e(kfound)) THEN
3862 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
3863 (coord_in(kfound) - coord_in(kfound-1)))
3865 IF (vartype == var_dir360) THEN
3866 field_out(i,j,k) = &
3867 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2
3868 ELSE IF (vartype == var_press) THEN
3869 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound
3871 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
3881 DEALLOCATE(coord_in,val_in)
3886 ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none'THEN
3888 field_out(:,:,:) = field_in(:,:,:)
3898 FUNCTION interp_var_360(v1, v2, w1, w2)
3899 REAL, INTENT(in) :: v1, v2, w1, w2
3900 REAL :: interp_var_360
3904 IF (abs(v1 - v2) > 180.) THEN
3912 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
3914 interp_var_360 = v1*w2 + v2*w1
3917 END FUNCTION interp_var_360
3920 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
3921 REAL, INTENT(in) :: v1(:,:)
3922 REAL, INTENT(in) :: l, h, res
3923 LOGICAL, INTENT(in), OPTIONAL :: mask(:,:)
3931 IF ((h - l) <= res) THEN
3936 IF ( PRESENT(mask)) THEN
3937 nl = count(v1 >= l .AND. v1 < m .AND. mask)
3938 nh = count(v1 >= m .AND. v1 < h .AND. mask)
3940 nl = count(v1 >= l .AND. v1 < m)
3941 nh = count(v1 >= m .AND. v1 < h)
3944 prev = find_prevailing_direction(v1, m, h, res)
3945 ELSE IF (nl > nh) THEN
3946 prev = find_prevailing_direction(v1, l, m, res)
3947 ELSE IF (nl == 0 .AND. nh == 0) THEN
3953 END FUNCTION find_prevailing_direction
3956 END SUBROUTINE grid_transform_compute
3964 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
3966 TYPE(grid_transform), INTENT(in) :: this
3967 REAL, INTENT(in) :: field_in(:,:)
3968 REAL, INTENT(out):: field_out(:,:,:)
3969 TYPE(vol7d_var), INTENT(in), OPTIONAL :: var
3970 REAL, INTENT(in), OPTIONAL, TARGET :: coord_3d_in(:,:,:)
3972 real, allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
3973 INTEGER :: inn_p, ier, k
3974 INTEGER :: innx, inny, innz, outnx, outny, outnz
3980 field_out(:,:,:) = rmiss
3982 IF (.NOT.this%valid) THEN
3983 call l4f_category_log(this%category,l4f_error, &
3984 "refusing to perform a non valid transformation")
3989 innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
3990 outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out
3993 IF (this%trans%trans_type == 'vertint') THEN
3995 IF (innz /= this%innz) THEN
3996 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation"
3997 CALL l4f_category_log(this%category,l4f_error, "inconsistent input shape: "
3998 t2c(this%innz)// " /= "//t2c(innz))
4003 IF (outnz /= this%outnz) THEN
4004 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation"
4005 CALL l4f_category_log(this%category,l4f_error, "inconsistent output shape: "
4006 t2c(this%outnz)// " /= "//t2c(outnz))
4011 IF (innx /= outnx .OR. inny /= outny) THEN
4012 CALL l4f_category_log(this%category,l4f_error, "vertical interpolation"
4013 CALL l4f_category_log(this%category,l4f_error, "inconsistent hor. sizes: "
4014 t2c(innx)// ","//t2c(inny)// " /= "//&
4015 t2c(outnx)// ","//t2c(outny))
4020 ELSE ! horizontal interpolation
4022 IF (innx /= this%innx .OR. inny /= this%inny) THEN
4023 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation"
4024 CALL l4f_category_log(this%category,l4f_error, "inconsistent input shape: "
4025 t2c(this%innx)// ","//t2c(this%inny)// " /= "//&
4026 t2c(innx)// ","//t2c(inny))
4031 IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
4032 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation"
4033 CALL l4f_category_log(this%category,l4f_error, "inconsistent output shape: "
4034 t2c(this%outnx)// ","//t2c(this%outny)// " /= "//&
4035 t2c(outnx)// ","//t2c(outny))
4040 IF (innz /= outnz) THEN
4041 CALL l4f_category_log(this%category,l4f_error, "horizontal interpolation"
4042 CALL l4f_category_log(this%category,l4f_error, "inconsistent vert. sizes: "
4043 t2c(innz)// " /= "//t2c(outnz))
4051 call l4f_category_log(this%category,l4f_debug, &
4052 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type) ':'
4053 trim(this%trans%sub_type))
4056 IF (this%trans%trans_type == 'inter') THEN
4058 IF (this%trans%sub_type == 'linear') THEN
4060 #ifdef HAVE_LIBNGMATH
4062 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny)
4064 inn_p = count(c_e(field_in(:,k)))
4066 CALL l4f_category_log(this%category,l4f_info, &
4067 "Number of sparse data points: "//t2c(inn_p)// ','//t2c( SIZE(field_in
4071 field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k)))
4072 x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k)))
4073 y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k)))
4075 IF (.NOT.this%trans%extrap) THEN
4076 CALL nnseti( 'ext', 0)
4077 CALL nnsetr( 'nul', rmiss)
4080 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4081 this%outnx, this%outny, real(this%inter_x(:,1)), &
4082 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4085 CALL l4f_category_log(this%category,l4f_error, &
4086 "Error code from NCAR natgrids: "//t2c(ier))
4092 CALL l4f_category_log(this%category,l4f_info, &
4093 "insufficient data in gridded region to triangulate")
4097 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4100 CALL l4f_category_log(this%category,l4f_error, &
4101 "libsim compiled without natgridd(ngmath ncarg library) ")
4108 .OR. ELSE IF (this%trans%trans_type == 'boxinter' &
4109 .OR. this%trans%trans_type == 'polyinter' &
4110 .OR. this%trans%trans_type == 'vertint' &
4111 this%trans%trans_type == 'metamorphosis') THEN ! use the grid-to-grid method
4113 CALL compute(this, &
4114 RESHAPE(field_in, (/SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
4119 END SUBROUTINE grid_transform_v7d_grid_compute
4122 ! Bilinear interpolation
4123 ! Effettua interpolazione bilineare dati i valori nei punti 1,2,3,4 e
4124 ! le coordinate dei punti 1 e 3 oltre a quelle del punto p dove viene
4125 ! valutato il campo.
4126 !_____________________________________________________________
4127 ! disposizione punti
4133 ! _____________________________________________________________
4134 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
4135 REAL,INTENT(in) :: z1,z2,z3,z4 ! Z values on the four points
4136 DOUBLE PRECISION,INTENT(in):: x1,y1 ! coordinate of the lower left point
4137 DOUBLE PRECISION,INTENT(in):: x3,y3 ! coordinate of the upper right point
4138 DOUBLE PRECISION,INTENT(in):: xp,yp ! coordinate of point where interpolate
|