libsim Versione 7.2.1

◆ 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]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 3095 del file grid_transform_class.F90.

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

Generated with Doxygen.