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