libsim Versione 7.2.6

◆ 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 )
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]thisgrid_transformation object
[in]field_ininput array
[out]field_outoutput array
[in]varphysical variable to be interpolated, if provided, some ad-hoc algorithms may be used where possible
[in]coord_3d_ininput vertical coordinate for vertical interpolation, if not provided by other means

Definizione alla linea 3120 del file grid_transform_class.F90.

3122
3123 ENDIF
3124
3125 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3126 this%trans%sub_type == 'stddevnm1') THEN
3127
3128 IF (this%trans%sub_type == 'stddev') THEN
3129 nm1 = .false.
3130 ELSE
3131 nm1 = .true.
3132 ENDIF
3133
3134 navg = this%trans%box_info%npx*this%trans%box_info%npy
3135 DO k = 1, innz
3136 jj = 0
3137 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3138 je = j+this%trans%box_info%npy-1
3139 jj = jj+1
3140 ii = 0
3141 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3142 ie = i+this%trans%box_info%npx-1
3143 ii = ii+1
3144 field_out(ii,jj,k) = stat_stddev( &
3145 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3146 ENDDO
3147 ENDDO
3148 ENDDO
3149
3150 ELSE IF (this%trans%sub_type == 'max') THEN
3151 DO k = 1, innz
3152 jj = 0
3153 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3154 je = j+this%trans%box_info%npy-1
3155 jj = jj+1
3156 ii = 0
3157 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3158 ie = i+this%trans%box_info%npx-1
3159 ii = ii+1
3160 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3161 IF (navg > 0) THEN
3162 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3163 mask=(field_in(i:ie,j:je,k) /= rmiss))
3164 ENDIF
3165 ENDDO
3166 ENDDO
3167 ENDDO
3168
3169 ELSE IF (this%trans%sub_type == 'min') THEN
3170 DO k = 1, innz
3171 jj = 0
3172 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3173 je = j+this%trans%box_info%npy-1
3174 jj = jj+1
3175 ii = 0
3176 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3177 ie = i+this%trans%box_info%npx-1
3178 ii = ii+1
3179 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3180 IF (navg > 0) THEN
3181 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3182 mask=(field_in(i:ie,j:je,k) /= rmiss))
3183 ENDIF
3184 ENDDO
3185 ENDDO
3186 ENDDO
3187
3188 ELSE IF (this%trans%sub_type == 'percentile') THEN
3189
3190 navg = this%trans%box_info%npx*this%trans%box_info%npy
3191 DO k = 1, innz
3192 jj = 0
3193 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3194 je = j+this%trans%box_info%npy-1
3195 jj = jj+1
3196 ii = 0
3197 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3198 ie = i+this%trans%box_info%npx-1
3199 ii = ii+1
3200 field_out(ii:ii,jj,k) = stat_percentile( &
3201 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3202 (/real(this%trans%stat_info%percentile)/))
3203 ENDDO
3204 ENDDO
3205 ENDDO
3206
3207 ELSE IF (this%trans%sub_type == 'frequency') THEN
3208
3209 DO k = 1, innz
3210 jj = 0
3211 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3212 je = j+this%trans%box_info%npy-1
3213 jj = jj+1
3214 ii = 0
3215 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3216 ie = i+this%trans%box_info%npx-1
3217 ii = ii+1
3218 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3219 IF (navg > 0) THEN
3220 field_out(ii,jj,k) = sum(interval_info_valid( &
3221 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3222 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3223 ENDIF
3224 ENDDO
3225 ENDDO
3226 ENDDO
3227
3228 ENDIF
3229
3230ELSE IF (this%trans%trans_type == 'inter') THEN
3231
3232 IF (this%trans%sub_type == 'near') THEN
3233
3234 DO k = 1, innz
3235 DO j = 1, this%outny
3236 DO i = 1, this%outnx
3237
3238 IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3239 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3240
3241 ENDDO
3242 ENDDO
3243 ENDDO
3244
3245 ELSE IF (this%trans%sub_type == 'bilin') THEN
3246
3247 DO k = 1, innz
3248 DO j = 1, this%outny
3249 DO i = 1, this%outnx
3250
3251 IF (c_e(this%inter_index_x(i,j))) THEN
3252
3253 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3254 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3255 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3256 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3257
3258 IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN
3259
3260 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3261 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3262 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3263 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3264
3265 xp=this%inter_xp(i,j)
3266 yp=this%inter_yp(i,j)
3267
3268 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3269
3270 ENDIF
3271 ENDIF
3272
3273 ENDDO
3274 ENDDO
3275 ENDDO
3276 ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
3277 DO k = 1, innz
3278 DO j = 1, this%outny
3279 DO i = 1, this%outnx
3280
3281 IF (c_e(this%inter_index_x(i,j))) THEN
3282
3283 IF(this%inter_index_x(i,j)-1>0)THEN
3284 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3285 ELSE
3286 z(1) = rmiss
3287 END IF
3288 IF(this%inter_index_x(i,j)+1<=this%outnx)THEN
3289 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3290 ELSE
3291 z(3) = rmiss
3292 END IF
3293 IF(this%inter_index_y(i,j)+1<=this%outny)THEN
3294 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3295 ELSE
3296 z(2) = rmiss
3297 END IF
3298 IF(this%inter_index_y(i,j)-1>0)THEN
3299 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3300 ELSE
3301 z(4) = rmiss
3302 END IF
3303 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3304
3305 ENDIF
3306
3307 ENDDO
3308 ENDDO
3309 ENDDO
3310
3311 ENDIF
3312ELSE IF (this%trans%trans_type == 'intersearch') THEN
3313
3314 likethis = this
3315 likethis%trans%trans_type = 'inter' ! fake type and make a recursive call to compute base field
3316 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3317 CALL getenv('LIBSIM_DISABLEOPTSEARCH', env_var)
3318 optsearch = len_trim(env_var) == 0
3319
3320 DO k = 1, innz
3321 IF ((.NOT.all(c_e(field_out(:,:,k)))) .AND. (any(c_e(field_in(:,:,k))))) THEN ! must fill some values
3322 DO j = 1, this%outny
3323 DO i = 1, this%outnx
3324 IF (.NOT.c_e(field_out(i,j,k))) THEN
3325 dist = huge(dist)
3326 nearcount = 0
3327 IF (optsearch) THEN ! optimized, error-prone algorithm
3328 ix = this%inter_index_x(i,j)
3329 iy = this%inter_index_y(i,j)
3330 DO s = 0, max(this%innx, this%inny)
3331 farenough = .true.
3332 DO m = iy-s, iy+s, max(2*s, 1) ! y loop on upper and lower frames
3333 IF (m < 1 .OR. m > this%inny) cycle
3334 DO l = max(1, ix-s), min(this%innx, ix+s) ! x loop on upper and lower frames
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
3338 dist = disttmp
3339 field_out(i,j,k) = field_in(l,m,k)
3340 nearcount = 1
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
3344 ENDIF
3345 ENDIF
3346 IF (disttmp < dist) farenough = .false.
3347 ENDDO
3348 ENDDO
3349 DO m = max(1, iy-s+1), min(this%inny, iy+s-1) ! y loop on left and right frames (avoid corners)
3350 DO l = ix-s, ix+s, 2*s ! x loop on left and right frames (exchange loops?)
3351 IF (l < 1 .OR. l > this%innx) cycle
3352 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3353 IF (c_e(field_in(l,m,k))) THEN
3354 IF (disttmp < dist) THEN
3355 dist = disttmp
3356 field_out(i,j,k) = field_in(l,m,k)
3357 nearcount = 1
3358 ELSE IF (disttmp == dist) THEN
3359 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3360 nearcount = nearcount + 1
3361 ENDIF
3362 ENDIF
3363 IF (disttmp < dist) farenough = .false.
3364 ENDDO
3365 ENDDO
3366 IF (s > 0 .AND. farenough) EXIT ! nearest point found, do not trust the same point, in case of bilin it could be not the nearest
3367 ENDDO
3368 ELSE ! linear, simple, slow algorithm
3369 DO m = 1, this%inny
3370 DO l = 1, this%innx
3371 IF (c_e(field_in(l,m,k))) THEN
3372 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3373 IF (disttmp < dist) THEN
3374 dist = disttmp
3375 field_out(i,j,k) = field_in(l,m,k)
3376 nearcount = 1
3377 ELSE IF (disttmp == dist) THEN
3378 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3379 nearcount = nearcount + 1
3380 ENDIF
3381 ENDIF
3382 ENDDO
3383 ENDDO
3384 ENDIF
3385! average points with same minimum distance
3386 IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount
3387 ENDIF
3388 ENDDO
3389 ENDDO
3390 ENDIF
3391 ENDDO
3392
3393ELSE IF (this%trans%trans_type == 'boxinter' &
3394 .OR. this%trans%trans_type == 'polyinter' &
3395 .OR. this%trans%trans_type == 'maskinter') THEN
3396
3397 IF (this%trans%sub_type == 'average') THEN
3398
3399 IF (vartype == var_dir360) THEN
3400 DO k = 1, innz
3401 DO j = 1, this%outny
3402 DO i = 1, this%outnx
3403 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3404 0.0, 360.0, 5.0, &
3405 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3406 ENDDO
3407 ENDDO
3408 ENDDO
3409
3410 ELSE
3411 ALLOCATE(nval(this%outnx, this%outny))
3412 field_out(:,:,:) = 0.0
3413 DO k = 1, innz
3414 nval(:,:) = 0
3415 DO j = 1, this%inny
3416 DO i = 1, this%innx
3417 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3418 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3419 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3420 field_in(i,j,k)
3421 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3422 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3423 ENDIF
3424 ENDDO
3425 ENDDO
3426 WHERE (nval(:,:) /= 0)
3427 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3428 ELSEWHERE
3429 field_out(:,:,k) = rmiss
3430 END WHERE
3431 ENDDO
3432 DEALLOCATE(nval)
3433 ENDIF
3434
3435 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3436 this%trans%sub_type == 'stddevnm1') THEN
3437
3438 IF (this%trans%sub_type == 'stddev') THEN
3439 nm1 = .false.
3440 ELSE
3441 nm1 = .true.
3442 ENDIF
3443 DO k = 1, innz
3444 DO j = 1, this%outny
3445 DO i = 1, this%outnx
3446! da paura
3447 field_out(i:i,j,k) = stat_stddev( &
3448 reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
3449 mask=reshape((this%inter_index_x == i .AND. &
3450 this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)), nm1=nm1)
3451 ENDDO
3452 ENDDO
3453 ENDDO
3454
3455 ELSE IF (this%trans%sub_type == 'max') THEN
3456
3457 DO k = 1, innz
3458 DO j = 1, this%inny
3459 DO i = 1, this%innx
3460 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3461 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3462 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3463 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3464 field_in(i,j,k))
3465 ELSE
3466 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3467 field_in(i,j,k)
3468 ENDIF
3469 ENDIF
3470 ENDDO
3471 ENDDO
3472 ENDDO
3473
3474
3475 ELSE IF (this%trans%sub_type == 'min') THEN
3476
3477 DO k = 1, innz
3478 DO j = 1, this%inny
3479 DO i = 1, this%innx
3480 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3481 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3482 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3483 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3484 field_in(i,j,k))
3485 ELSE
3486 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3487 field_in(i,j,k)
3488 ENDIF
3489 ENDIF
3490 ENDDO
3491 ENDDO
3492 ENDDO
3493
3494 ELSE IF (this%trans%sub_type == 'percentile') THEN
3495
3496 DO k = 1, innz
3497 DO j = 1, this%outny
3498 DO i = 1, this%outnx
3499! da paura
3500 field_out(i:i,j,k) = stat_percentile( &
3501 reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
3502 (/real(this%trans%stat_info%percentile)/), &
3503 mask=reshape((this%inter_index_x == i .AND. &
3504 this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)))
3505 ENDDO
3506 ENDDO
3507 ENDDO
3508
3509 ELSE IF (this%trans%sub_type == 'frequency') THEN
3510
3511 ALLOCATE(nval(this%outnx, this%outny))
3512 field_out(:,:,:) = 0.0
3513 DO k = 1, innz
3514 nval(:,:) = 0
3515 DO j = 1, this%inny
3516 DO i = 1, this%innx
3517 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3518 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3519 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3520 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3521 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3522 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3523 ENDIF
3524 ENDDO
3525 ENDDO
3526 WHERE (nval(:,:) /= 0)
3527 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3528 ELSEWHERE
3529 field_out(:,:,k) = rmiss
3530 END WHERE
3531 ENDDO
3532 DEALLOCATE(nval)
3533
3534 ENDIF
3535
3536ELSE IF (this%trans%trans_type == 'stencilinter') THEN
3537 np = SIZE(this%stencil,1)/2
3538 ns = SIZE(this%stencil)
3539
3540 IF (this%trans%sub_type == 'average') THEN
3541
3542 IF (vartype == var_dir360) THEN
3543 DO k = 1, innz
3544 DO j = 1, this%outny
3545 DO i = 1, this%outnx
3546 IF (c_e(this%inter_index_x(i,j))) THEN
3547 i1 = this%inter_index_x(i,j) - np
3548 i2 = this%inter_index_x(i,j) + np
3549 j1 = this%inter_index_y(i,j) - np
3550 j2 = this%inter_index_y(i,j) + np
3551 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3552 0.0, 360.0, 5.0, &
3553 mask=this%stencil(:,:)) ! simpler and equivalent form
3554! mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3555 ENDIF
3556 ENDDO
3557 ENDDO
3558 ENDDO
3559
3560 ELSE
3561!$OMP PARALLEL DEFAULT(SHARED)
3562!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3563 DO k = 1, innz
3564 DO j = 1, this%outny
3565 DO i = 1, this%outnx
3566 IF (c_e(this%inter_index_x(i,j))) THEN
3567 i1 = this%inter_index_x(i,j) - np
3568 i2 = this%inter_index_x(i,j) + np
3569 j1 = this%inter_index_y(i,j) - np
3570 j2 = this%inter_index_y(i,j) + np
3571 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3572 IF (n > 0) THEN
3573 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3574 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3575 ENDIF
3576 ENDIF
3577 ENDDO
3578 ENDDO
3579 ENDDO
3580!$OMP END PARALLEL
3581 ENDIF
3582
3583 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3584 this%trans%sub_type == 'stddevnm1') THEN
3585
3586 IF (this%trans%sub_type == 'stddev') THEN
3587 nm1 = .false.
3588 ELSE
3589 nm1 = .true.
3590 ENDIF
3591
3592!$OMP PARALLEL DEFAULT(SHARED)
3593!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3594 DO k = 1, innz
3595 DO j = 1, this%outny
3596 DO i = 1, this%outnx
3597 IF (c_e(this%inter_index_x(i,j))) THEN
3598 i1 = this%inter_index_x(i,j) - np
3599 i2 = this%inter_index_x(i,j) + np
3600 j1 = this%inter_index_y(i,j) - np
3601 j2 = this%inter_index_y(i,j) + np
3602! da paura
3603 field_out(i:i,j,k) = stat_stddev( &
3604 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3605 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3606 this%stencil(:,:), (/ns/)), nm1=nm1)
3607 ENDIF
3608 ENDDO
3609 ENDDO
3610 ENDDO
3611!$OMP END PARALLEL
3612
3613 ELSE IF (this%trans%sub_type == 'max') THEN
3614
3615!$OMP PARALLEL DEFAULT(SHARED)
3616!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3617 DO k = 1, innz
3618 DO j = 1, this%outny
3619 DO i = 1, this%outnx
3620 IF (c_e(this%inter_index_x(i,j))) THEN
3621 i1 = this%inter_index_x(i,j) - np
3622 i2 = this%inter_index_x(i,j) + np
3623 j1 = this%inter_index_y(i,j) - np
3624 j2 = this%inter_index_y(i,j) + np
3625 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3626 IF (n > 0) THEN
3627 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3628 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3629 ENDIF
3630 ENDIF
3631 ENDDO
3632 ENDDO
3633 ENDDO
3634!$OMP END PARALLEL
3635
3636 ELSE IF (this%trans%sub_type == 'min') THEN
3637
3638!$OMP PARALLEL DEFAULT(SHARED)
3639!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3640 DO k = 1, innz
3641 DO j = 1, this%outny
3642 DO i = 1, this%outnx
3643 IF (c_e(this%inter_index_x(i,j))) THEN
3644 i1 = this%inter_index_x(i,j) - np
3645 i2 = this%inter_index_x(i,j) + np
3646 j1 = this%inter_index_y(i,j) - np
3647 j2 = this%inter_index_y(i,j) + np
3648 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3649 IF (n > 0) THEN
3650 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3651 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3652 ENDIF
3653 ENDIF
3654 ENDDO
3655 ENDDO
3656 ENDDO
3657!$OMP END PARALLEL
3658
3659 ELSE IF (this%trans%sub_type == 'percentile') THEN
3660
3661!$OMP PARALLEL DEFAULT(SHARED)
3662!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3663 DO k = 1, innz
3664 DO j = 1, this%outny
3665 DO i = 1, this%outnx
3666 IF (c_e(this%inter_index_x(i,j))) THEN
3667 i1 = this%inter_index_x(i,j) - np
3668 i2 = this%inter_index_x(i,j) + np
3669 j1 = this%inter_index_y(i,j) - np
3670 j2 = this%inter_index_y(i,j) + np
3671! da paura
3672 field_out(i:i,j,k) = stat_percentile( &
3673 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3674 (/real(this%trans%stat_info%percentile)/), &
3675 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3676 this%stencil(:,:), (/ns/)))
3677 ENDIF
3678 ENDDO
3679 ENDDO
3680 ENDDO
3681!$OMP END PARALLEL
3682
3683 ELSE IF (this%trans%sub_type == 'frequency') THEN
3684
3685!$OMP PARALLEL DEFAULT(SHARED)
3686!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3687 DO k = 1, innz
3688 DO j = 1, this%outny
3689 DO i = 1, this%outnx
3690 IF (c_e(this%inter_index_x(i,j))) THEN
3691 i1 = this%inter_index_x(i,j) - np
3692 i2 = this%inter_index_x(i,j) + np
3693 j1 = this%inter_index_y(i,j) - np
3694 j2 = this%inter_index_y(i,j) + np
3695 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3696 IF (n > 0) THEN
3697 field_out(i,j,k) = sum(interval_info_valid( &
3698 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3699 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3700 ENDIF
3701 ENDIF
3702 ENDDO
3703 ENDDO
3704 ENDDO
3705!$OMP END PARALLEL
3706
3707 ENDIF
3708
3709ELSE IF (this%trans%trans_type == 'maskgen') THEN
3710
3711 DO k = 1, innz
3712 WHERE(c_e(this%inter_index_x(:,:)))
3713 field_out(:,:,k) = real(this%inter_index_x(:,:))
3714 END WHERE
3715 ENDDO
3716
3717ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
3718
3719 IF (this%trans%sub_type == 'all') THEN
3720
3721 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3722
3723 ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' &
3724 .OR. this%trans%sub_type == 'mask') THEN
3725
3726 DO k = 1, innz
3727! this is to sparse-points only, so field_out(:,1,k) is acceptable
3728 field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)))
3729 ENDDO
3730
3731 ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
3732 this%trans%sub_type == 'maskinvalid') THEN
3733
3734 DO k = 1, innz
3735 WHERE (this%point_mask(:,:))
3736 field_out(:,:,k) = field_in(:,:,k)
3737 END WHERE
3738 ENDDO
3739
3740 ELSE IF (this%trans%sub_type == 'lemaskinvalid') THEN
3741
3742 DO k = 1, innz
3743 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3744 field_out(:,:,k) = field_in(:,:,k)
3745 ELSEWHERE
3746 field_out(:,:,k) = rmiss
3747 END WHERE
3748 ENDDO
3749
3750 ELSE IF (this%trans%sub_type == 'ltmaskinvalid') THEN
3751
3752 DO k = 1, innz
3753 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3754 field_out(:,:,k) = field_in(:,:,k)
3755 ELSEWHERE
3756 field_out(:,:,k) = rmiss
3757 END WHERE
3758 ENDDO
3759
3760 ELSE IF (this%trans%sub_type == 'gemaskinvalid') THEN
3761
3762 DO k = 1, innz
3763 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3764 field_out(:,:,k) = field_in(:,:,k)
3765 ELSEWHERE
3766 field_out(:,:,k) = rmiss
3767 END WHERE
3768 ENDDO
3769
3770 ELSE IF (this%trans%sub_type == 'gtmaskinvalid') THEN
3771
3772 DO k = 1, innz
3773 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3774 field_out(:,:,k) = field_in(:,:,k)
3775 ELSEWHERE
3776 field_out(:,:,k) = rmiss
3777 END WHERE
3778 ENDDO
3779
3780 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
3781
3782 DO k = 1, innz
3783 WHERE (c_e(field_in(:,:,k)))
3784 field_out(:,:,k) = field_in(:,:,k)
3785 ELSE WHERE
3786 field_out(:,:,k) = this%val1
3787 END WHERE
3788 ENDDO
3789
3790 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
3791 IF (c_e(this%val1) .AND. c_e(this%val2)) THEN
3792 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3793 .AND. field_in(:,:,:) <= this%val2)
3794 field_out(:,:,:) = rmiss
3795 ELSE WHERE
3796 field_out(:,:,:) = field_in(:,:,:)
3797 END WHERE
3798 ELSE IF (c_e(this%val1)) THEN
3799 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3800 field_out(:,:,:) = rmiss
3801 ELSE WHERE
3802 field_out(:,:,:) = field_in(:,:,:)
3803 END WHERE
3804 ELSE IF (c_e(this%val2)) THEN
3805 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3806 field_out(:,:,:) = rmiss
3807 ELSE WHERE
3808 field_out(:,:,:) = field_in(:,:,:)
3809 END WHERE
3810 ENDIF
3811 ENDIF
3812
3813ELSE IF (this%trans%trans_type == 'vertint') THEN
3814
3815 IF (this%trans%sub_type == 'linear') THEN
3816
3817 alloc_coord_3d_in_act = .false.
3818 IF (ASSOCIATED(this%inter_index_z)) THEN
3819
3820 DO k = 1, outnz
3821 IF (c_e(this%inter_index_z(k))) THEN
3822 z1 = real(this%inter_zp(k)) ! weight for k+1
3823 z2 = real(1.0d0 - this%inter_zp(k)) ! weight for k
3824 SELECT CASE(vartype)
3825
3826 CASE(var_dir360)
3827 DO j = 1, inny
3828 DO i = 1, innx
3829 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3830 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3831 field_out(i,j,k) = &
3832 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3833 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3834 ENDIF
3835 ENDDO
3836 ENDDO
3837
3838 CASE(var_press)
3839 DO j = 1, inny
3840 DO i = 1, innx
3841 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3842 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3843 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3844 field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN
3845 field_out(i,j,k) = exp( &
3846 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3847 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3848 ENDIF
3849 ENDDO
3850 ENDDO
3851
3852 CASE default
3853 DO j = 1, inny
3854 DO i = 1, innx
3855 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3856 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3857 field_out(i,j,k) = &
3858 field_in(i,j,this%inter_index_z(k))*z2 + &
3859 field_in(i,j,this%inter_index_z(k)+1)*z1
3860 ENDIF
3861 ENDDO
3862 ENDDO
3863
3864 END SELECT
3865
3866 ENDIF
3867 ENDDO
3868
3869 ELSE ! use coord_3d_in
3870
3871 IF (ALLOCATED(this%coord_3d_in)) THEN
3872 coord_3d_in_act => this%coord_3d_in
3873#ifdef DEBUG
3874 CALL l4f_category_log(this%category,l4f_debug, &
3875 "Using external vertical coordinate for vertint")
3876#endif
3877 ELSE
3878 IF (PRESENT(coord_3d_in)) THEN
3879#ifdef DEBUG
3880 CALL l4f_category_log(this%category,l4f_debug, &
3881 "Using internal vertical coordinate for vertint")
3882#endif
3883 IF (this%dolog) THEN
3884 ALLOCATE(coord_3d_in_act(SIZE(coord_3d_in,1), &
3885 SIZE(coord_3d_in,2), SIZE(coord_3d_in,3)))
3886 alloc_coord_3d_in_act = .true.
3887 WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3888 coord_3d_in_act = log(coord_3d_in)
3889 ELSEWHERE
3890 coord_3d_in_act = rmiss
3891 END WHERE
3892 ELSE
3893 coord_3d_in_act => coord_3d_in
3894 ENDIF
3895 ELSE
3896 CALL l4f_category_log(this%category,l4f_error, &
3897 "Missing internal and external vertical coordinate for vertint")
3898 CALL raise_error()
3899 RETURN
3900 ENDIF
3901 ENDIF
3902
3903 inused = SIZE(coord_3d_in_act,3)
3904 IF (inused < 2) RETURN ! to avoid algorithm failure
3905 kkcache = 1
3906
3907 DO k = 1, outnz
3908 IF (.NOT.c_e(this%vcoord_out(k))) cycle
3909 DO j = 1, inny
3910 DO i = 1, innx
3911 kfound = imiss
3912 DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
3913 kkup = kkcache + kk
3914 kkdown = kkcache - kk + 1
3915
3916 IF (kkdown >= 1) THEN ! search down
3917 IF (this%vcoord_out(k) >= &
3918 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3919 this%vcoord_out(k) < &
3920 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1))) THEN
3921 kkcache = kkdown
3922 kfoundin = kkcache
3923 kfound = kkcache + this%levshift
3924 EXIT ! kk
3925 ENDIF
3926 ENDIF
3927 IF (kkup < inused) THEN ! search up
3928 IF (this%vcoord_out(k) >= &
3929 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3930 this%vcoord_out(k) < &
3931 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1))) THEN
3932 kkcache = kkup
3933 kfoundin = kkcache
3934 kfound = kkcache + this%levshift
3935 EXIT ! kk
3936 ENDIF
3937 ENDIF
3938
3939 ENDDO
3940
3941 IF (c_e(kfound)) THEN
3942 IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound+1))) THEN
3943 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3944 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3945 z2 = 1.0 - z1
3946 SELECT CASE(vartype)
3947
3948 CASE(var_dir360)
3949 field_out(i,j,k) = &
3950 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3951 CASE(var_press)
3952 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3953 log(field_in(i,j,kfound+1))*z1)
3954
3955 CASE default
3956 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3957 END SELECT
3958
3959 ENDIF
3960 ELSE
3961 ENDIF
3962 ENDDO ! i
3963 ENDDO ! j
3964 ENDDO ! k
3965 IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
3966 ENDIF
3967
3968 ELSE IF (this%trans%sub_type == 'linearsparse') THEN
3969
3970
3971 IF (.NOT.ASSOCIATED(this%vcoord_in) .AND. .NOT.PRESENT(coord_3d_in)) THEN
3972 CALL l4f_category_log(this%category,l4f_error, &
3973 "linearsparse interpolation, no input vert coord available")
3974 RETURN
3975 ENDIF
3976
3977 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3978 DO j = 1, inny
3979 DO i = 1, innx
3980
3981 IF (ASSOCIATED(this%vcoord_in)) THEN
3982 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3983 .AND. c_e(this%vcoord_in(:))
3984 ELSE
3985 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3986 .AND. c_e(coord_3d_in(i,j,:))
3987 ENDIF
3988
3989 IF (vartype == var_press) THEN
3990 mask_in(:) = mask_in(:) .AND. &
3991 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3992 ENDIF
3993 inused = count(mask_in)
3994 IF (inused > 1) THEN
3995 IF (ASSOCIATED(this%vcoord_in)) THEN
3996 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3997 ELSE
3998 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3999 ENDIF
4000 IF (vartype == var_press) THEN
4001 val_in(1:inused) = log(pack( &
4002 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4003 mask=mask_in))
4004 ELSE
4005 val_in(1:inused) = pack( &
4006 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4007 mask=mask_in)
4008 ENDIF
4009 kkcache = 1
4010 DO k = 1, outnz
4011
4012 kfound = imiss
4013 DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
4014 kkup = kkcache + kk
4015 kkdown = kkcache - kk + 1
4016
4017 IF (kkdown >= 1) THEN ! search down
4018 IF (this%vcoord_out(k) >= &
4019 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
4020 this%vcoord_out(k) < &
4021 max(coord_in(kkdown), coord_in(kkdown+1))) THEN
4022 kkcache = kkdown
4023 kfoundin = kkcache
4024 kfound = kkcache
4025 EXIT ! kk
4026 ENDIF
4027 ENDIF
4028 IF (kkup < inused) THEN ! search up
4029 IF (this%vcoord_out(k) >= &
4030 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
4031 this%vcoord_out(k) < &
4032 max(coord_in(kkup), coord_in(kkup+1))) THEN
4033 kkcache = kkup
4034 kfoundin = kkcache
4035 kfound = kkcache
4036 EXIT ! kk
4037 ENDIF
4038 ENDIF
4039
4040 ENDDO
4041
4042 IF (c_e(kfound)) THEN
4043 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
4044 (coord_in(kfound) - coord_in(kfound-1)))
4045 z2 = 1.0 - z1
4046 IF (vartype == var_dir360) THEN
4047 field_out(i,j,k) = &
4048 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
4049 ELSE IF (vartype == var_press) THEN
4050 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
4051 ELSE
4052 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
4053 ENDIF
4054 ENDIF
4055
4056 ENDDO
4057
4058 ENDIF
4059
4060 ENDDO
4061 ENDDO
4062 DEALLOCATE(coord_in,val_in)
4063
4064
4065 ENDIF
4066
4067ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN
4068
4069 field_out(:,:,:) = field_in(:,:,:)
4070
4071ENDIF
4072
4073
4074CONTAINS
4075
4076
4077! internal function for interpolating directions from 0 to 360 degree
4078! hope it is inlined by the compiler
4079FUNCTION interp_var_360(v1, v2, w1, w2)
4080REAL,INTENT(in) :: v1, v2, w1, w2
4081REAL :: interp_var_360
4082
4083REAL :: lv1, lv2
4084
4085IF (abs(v1 - v2) > 180.) THEN
4086 IF (v1 > v2) THEN
4087 lv1 = v1 - 360.
4088 lv2 = v2
4089 ELSE
4090 lv1 = v1
4091 lv2 = v2 - 360.
4092 ENDIF
4093 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4094ELSE
4095 interp_var_360 = v1*w2 + v2*w1
4096ENDIF
4097
4098END FUNCTION interp_var_360
4099
4100
4101RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
4102REAL,INTENT(in) :: v1(:,:)
4103REAL,INTENT(in) :: l, h, res
4104LOGICAL,INTENT(in),OPTIONAL :: mask(:,:)
4105REAL :: prev
4106
4107REAL :: m
4108INTEGER :: nh, nl
4109!REAL,PARAMETER :: res = 1.0
4110
4111m = (l + h)/2.
4112IF ((h - l) <= res) THEN
4113 prev = m
4114 RETURN
4115ENDIF
4116
4117IF (PRESENT(mask)) THEN
4118 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4119 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4120ELSE
4121 nl = count(v1 >= l .AND. v1 < m)
4122 nh = count(v1 >= m .AND. v1 < h)
4123ENDIF
4124IF (nh > nl) THEN
4125 prev = find_prevailing_direction(v1, m, h, res)
4126ELSE IF (nl > nh) THEN
4127 prev = find_prevailing_direction(v1, l, m, res)
4128ELSE IF (nl == 0 .AND. nh == 0) THEN
4129 prev = rmiss
4130ELSE
4131 prev = m
4132ENDIF
4133
4134END FUNCTION find_prevailing_direction
4135
4136
4137END SUBROUTINE grid_transform_compute
4138
4139
4140!> Compute the output data array from input data array according to
4141!! the defined transformation. The \a grid_transform object \a this
4142!! must have been properly initialised, so that it contains all the
4143!! information needed for computing the transformation. This is the
4144!! sparse points-to-grid and sparse points-to-sparse points version.
4145SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4146 coord_3d_in)
4147TYPE(grid_transform),INTENT(in) :: this !< grid_tranform object
4148REAL, INTENT(in) :: field_in(:,:) !< input array
4149REAL, INTENT(out):: field_out(:,:,:) !< output array
4150TYPE(vol7d_var),INTENT(in),OPTIONAL :: var !< physical variable to be interpolated, if provided, some ad-hoc algorithms may be used where possible
4151REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:) !< input vertical coordinate for vertical interpolation, if not provided by other means
4152
4153real,allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4154INTEGER :: inn_p, ier, k
4155INTEGER :: innx, inny, innz, outnx, outny, outnz
4156
4157#ifdef DEBUG
4158call l4f_category_log(this%category,l4f_debug,"start v7d_grid_transform_compute")
4159#endif
4160
4161field_out(:,:,:) = rmiss
4162
4163IF (.NOT.this%valid) THEN
4164 call l4f_category_log(this%category,l4f_error, &
4165 "refusing to perform a non valid transformation")
4166 call raise_error()
4167 RETURN
4168ENDIF
4169
4170innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
4171outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
4172
4173! check size of field_in, field_out
4174IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation
4175
4176 IF (innz /= this%innz) THEN
4177 CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4178 CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
4179 t2c(this%innz)//" /= "//t2c(innz))
4180 CALL raise_error()
4181 RETURN
4182 ENDIF
4183
4184 IF (outnz /= this%outnz) THEN
4185 CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4186 CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
4187 t2c(this%outnz)//" /= "//t2c(outnz))
4188 CALL raise_error()
4189 RETURN
4190 ENDIF
4191
4192 IF (innx /= outnx .OR. inny /= outny) THEN
4193 CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4194 CALL l4f_category_log(this%category,l4f_error,"inconsistent hor. sizes: "//&
4195 t2c(innx)//","//t2c(inny)//" /= "//&
4196 t2c(outnx)//","//t2c(outny))
4197 CALL raise_error()
4198 RETURN
4199 ENDIF
4200
4201ELSE ! horizontal interpolation
4202
4203 IF (innx /= this%innx .OR. inny /= this%inny) THEN
4204 CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4205 CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
4206 t2c(this%innx)//","//t2c(this%inny)//" /= "//&
4207 t2c(innx)//","//t2c(inny))
4208 CALL raise_error()
4209 RETURN
4210 ENDIF
4211
4212 IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
4213 CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4214 CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
4215 t2c(this%outnx)//","//t2c(this%outny)//" /= "//&
4216 t2c(outnx)//","//t2c(outny))
4217 CALL raise_error()
4218 RETURN
4219 ENDIF
4220
4221 IF (innz /= outnz) THEN
4222 CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4223 CALL l4f_category_log(this%category,l4f_error,"inconsistent vert. sizes: "//&
4224 t2c(innz)//" /= "//t2c(outnz))
4225 CALL raise_error()
4226 RETURN
4227 ENDIF
4228
4229ENDIF
4230
4231#ifdef DEBUG
4232 call l4f_category_log(this%category,l4f_debug, &
4233 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//':'// &
4234 trim(this%trans%sub_type))
4235#endif
4236
4237IF (this%trans%trans_type == 'inter') THEN
4238
4239 IF (this%trans%sub_type == 'linear') THEN
4240
4241#ifdef HAVE_LIBNGMATH
4242! optimization, allocate only once with a safe size
4243 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4244 DO k = 1, innz
4245 inn_p = count(c_e(field_in(:,k)))
4246
4247 CALL l4f_category_log(this%category,l4f_info, &
4248 "Number of sparse data points: "//t2c(inn_p)//','//t2c(SIZE(field_in(:,k))))
4249
4250 IF (inn_p > 2) THEN
4251
4252 field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k)))
4253 x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k)))
4254 y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k)))
4255
4256 IF (.NOT.this%trans%extrap) THEN
4257 CALL nnseti('ext', 0) ! 0 = no extrapolation
4258 CALL nnsetr('nul', rmiss)
4259 ENDIF
4260
4261 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, & ! (1:inn_p) omitted
4262 this%outnx, this%outny, real(this%inter_x(:,1)), & ! no f90 interface
4263 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4264
4265 IF (ier /= 0) THEN
4266 CALL l4f_category_log(this%category,l4f_error, &
4267 "Error code from NCAR natgrids: "//t2c(ier))
4268 CALL raise_error()
4269 EXIT
4270 ENDIF ! exit loop to deallocate
4271 ELSE
4272
4273 CALL l4f_category_log(this%category,l4f_info, &
4274 "insufficient data in gridded region to triangulate")
4275
4276 ENDIF
4277 ENDDO
4278 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4279
4280#else
4281 CALL l4f_category_log(this%category,l4f_error, &
4282 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4283 CALL raise_error()
4284 RETURN
4285#endif
4286
4287 ENDIF
4288
4289ELSE IF (this%trans%trans_type == 'boxinter' .OR. &
4290 this%trans%trans_type == 'polyinter' .OR. &
4291 this%trans%trans_type == 'vertint' .OR. &
4292 this%trans%trans_type == 'metamorphosis') THEN ! use the grid-to-grid method
4293
4294 CALL compute(this, &
4295 reshape(field_in, (/SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
4296 coord_3d_in)
4297
4298ENDIF
4299
4300END SUBROUTINE grid_transform_v7d_grid_compute
4301
4302
4303! Bilinear interpolation
4304! Effettua interpolazione bilineare dati i valori nei punti 1,2,3,4 e
4305! le coordinate dei punti 1 e 3 oltre a quelle del punto p dove viene
4306! valutato il campo.
4307!_____________________________________________________________
4308! disposizione punti
4309! 4 3
4310!
4311! p
4312!
4313! 1 2
4314! _____________________________________________________________
4315ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
4316REAL,INTENT(in) :: z1,z2,z3,z4 ! Z values on the four points
4317DOUBLE PRECISION,INTENT(in):: x1,y1 ! coordinate of the lower left point

Generated with Doxygen.