libsim Versione 7.2.4
|
◆ grid_transform_compute()
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.
Definizione alla linea 3103 del file grid_transform_class.F90. 3105
3106 ENDIF
3107
3108 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3109 this%trans%sub_type == 'stddevnm1') THEN
3110
3111 IF (this%trans%sub_type == 'stddev') THEN
3112 nm1 = .false.
3113 ELSE
3114 nm1 = .true.
3115 ENDIF
3116
3117 navg = this%trans%box_info%npx*this%trans%box_info%npy
3118 DO k = 1, innz
3119 jj = 0
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
3122 jj = jj+1
3123 ii = 0
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
3126 ii = ii+1
3127 field_out(ii,jj,k) = stat_stddev( &
3128 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3129 ENDDO
3130 ENDDO
3131 ENDDO
3132
3133 ELSE IF (this%trans%sub_type == 'max') THEN
3134 DO k = 1, innz
3135 jj = 0
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
3138 jj = jj+1
3139 ii = 0
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
3142 ii = ii+1
3143 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3144 IF (navg > 0) THEN
3145 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3146 mask=(field_in(i:ie,j:je,k) /= rmiss))
3147 ENDIF
3148 ENDDO
3149 ENDDO
3150 ENDDO
3151
3152 ELSE IF (this%trans%sub_type == 'min') THEN
3153 DO k = 1, innz
3154 jj = 0
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
3157 jj = jj+1
3158 ii = 0
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
3161 ii = ii+1
3162 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3163 IF (navg > 0) THEN
3164 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3165 mask=(field_in(i:ie,j:je,k) /= rmiss))
3166 ENDIF
3167 ENDDO
3168 ENDDO
3169 ENDDO
3170
3171 ELSE IF (this%trans%sub_type == 'percentile') THEN
3172
3173 navg = this%trans%box_info%npx*this%trans%box_info%npy
3174 DO k = 1, innz
3175 jj = 0
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
3178 jj = jj+1
3179 ii = 0
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
3182 ii = ii+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)/))
3186 ENDDO
3187 ENDDO
3188 ENDDO
3189
3190 ELSE IF (this%trans%sub_type == 'frequency') THEN
3191
3192 DO k = 1, innz
3193 jj = 0
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
3196 jj = jj+1
3197 ii = 0
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
3200 ii = ii+1
3201 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3202 IF (navg > 0) THEN
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
3206 ENDIF
3207 ENDDO
3208 ENDDO
3209 ENDDO
3210
3211 ENDIF
3212
3213ELSE IF (this%trans%trans_type == 'inter') THEN
3214
3215 IF (this%trans%sub_type == 'near') THEN
3216
3217 DO k = 1, innz
3218 DO j = 1, this%outny
3219 DO i = 1, this%outnx
3220
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)
3223
3224 ENDDO
3225 ENDDO
3226 ENDDO
3227
3228 ELSE IF (this%trans%sub_type == 'bilin') THEN
3229
3230 DO k = 1, innz
3231 DO j = 1, this%outny
3232 DO i = 1, this%outnx
3233
3234 IF (c_e(this%inter_index_x(i,j))) THEN
3235
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)
3240
3241 IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN
3242
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)
3247
3248 xp=this%inter_xp(i,j)
3249 yp=this%inter_yp(i,j)
3250
3251 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3252
3253 ENDIF
3254 ENDIF
3255
3256 ENDDO
3257 ENDDO
3258 ENDDO
3259 ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
3260 DO k = 1, innz
3261 DO j = 1, this%outny
3262 DO i = 1, this%outnx
3263
3264 IF (c_e(this%inter_index_x(i,j))) THEN
3265
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)
3268 ELSE
3269 z(1) = rmiss
3270 END IF
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)
3273 ELSE
3274 z(3) = rmiss
3275 END IF
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)
3278 ELSE
3279 z(2) = rmiss
3280 END IF
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)
3283 ELSE
3284 z(4) = rmiss
3285 END IF
3286 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3287
3288 ENDIF
3289
3290 ENDDO
3291 ENDDO
3292 ENDDO
3293
3294 ENDIF
3295ELSE IF (this%trans%trans_type == 'intersearch') THEN
3296
3297 likethis = this
3298 likethis%trans%trans_type = 'inter' ! fake type and make a recursive call to compute base field
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
3302
3303 DO k = 1, innz
3304 IF ((.NOT.all(c_e(field_out(:,:,k)))) .AND. (any(c_e(field_in(:,:,k))))) THEN ! must fill some values
3305 DO j = 1, this%outny
3306 DO i = 1, this%outnx
3307 IF (.NOT.c_e(field_out(i,j,k))) THEN
3308 dist = huge(dist)
3309 nearcount = 0
3310 IF (optsearch) THEN ! optimized, error-prone algorithm
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)
3314 farenough = .true.
3315 DO m = iy-s, iy+s, max(2*s, 1) ! y loop on upper and lower frames
3316 IF (m < 1 .OR. m > this%inny) cycle
3317 DO l = max(1, ix-s), min(this%innx, ix+s) ! x loop on upper and lower frames
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
3321 dist = disttmp
3322 field_out(i,j,k) = field_in(l,m,k)
3323 nearcount = 1
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
3327 ENDIF
3328 ENDIF
3329 IF (disttmp < dist) farenough = .false.
3330 ENDDO
3331 ENDDO
3332 DO m = max(1, iy-s+1), min(this%inny, iy+s-1) ! y loop on left and right frames (avoid corners)
3333 DO l = ix-s, ix+s, 2*s ! x loop on left and right frames (exchange loops?)
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
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 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
3350 ENDDO
3351 ELSE ! linear, simple, slow algorithm
3352 DO m = 1, this%inny
3353 DO l = 1, this%innx
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
3357 dist = disttmp
3358 field_out(i,j,k) = field_in(l,m,k)
3359 nearcount = 1
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
3363 ENDIF
3364 ENDIF
3365 ENDDO
3366 ENDDO
3367 ENDIF
3368! average points with same minimum distance
3369 IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount
3370 ENDIF
3371 ENDDO
3372 ENDDO
3373 ENDIF
3374 ENDDO
3375
3376ELSE IF (this%trans%trans_type == 'boxinter' &
3377 .OR. this%trans%trans_type == 'polyinter' &
3378 .OR. this%trans%trans_type == 'maskinter') THEN
3379
3380 IF (this%trans%sub_type == 'average') THEN
3381
3382 IF (vartype == var_dir360) THEN
3383 DO k = 1, innz
3384 DO j = 1, this%outny
3385 DO i = 1, this%outnx
3386 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3387 0.0, 360.0, 5.0, &
3388 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3389 ENDDO
3390 ENDDO
3391 ENDDO
3392
3393 ELSE
3394 ALLOCATE(nval(this%outnx, this%outny))
3395 field_out(:,:,:) = 0.0
3396 DO k = 1, innz
3397 nval(:,:) = 0
3398 DO j = 1, this%inny
3399 DO i = 1, this%innx
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) + &
3403 field_in(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
3406 ENDIF
3407 ENDDO
3408 ENDDO
3409 WHERE (nval(:,:) /= 0)
3410 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3411 ELSEWHERE
3412 field_out(:,:,k) = rmiss
3413 END WHERE
3414 ENDDO
3415 DEALLOCATE(nval)
3416 ENDIF
3417
3418 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3419 this%trans%sub_type == 'stddevnm1') THEN
3420
3421 IF (this%trans%sub_type == 'stddev') THEN
3422 nm1 = .false.
3423 ELSE
3424 nm1 = .true.
3425 ENDIF
3426 DO k = 1, innz
3427 DO j = 1, this%outny
3428 DO i = 1, this%outnx
3429! da paura
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)
3434 ENDDO
3435 ENDDO
3436 ENDDO
3437
3438 ELSE IF (this%trans%sub_type == 'max') THEN
3439
3440 DO k = 1, innz
3441 DO j = 1, this%inny
3442 DO i = 1, this%innx
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), &
3447 field_in(i,j,k))
3448 ELSE
3449 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3450 field_in(i,j,k)
3451 ENDIF
3452 ENDIF
3453 ENDDO
3454 ENDDO
3455 ENDDO
3456
3457
3458 ELSE IF (this%trans%sub_type == 'min') THEN
3459
3460 DO k = 1, innz
3461 DO j = 1, this%inny
3462 DO i = 1, this%innx
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), &
3467 field_in(i,j,k))
3468 ELSE
3469 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3470 field_in(i,j,k)
3471 ENDIF
3472 ENDIF
3473 ENDDO
3474 ENDDO
3475 ENDDO
3476
3477 ELSE IF (this%trans%sub_type == 'percentile') THEN
3478
3479 DO k = 1, innz
3480 DO j = 1, this%outny
3481 DO i = 1, this%outnx
3482! da paura
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))/)))
3488 ENDDO
3489 ENDDO
3490 ENDDO
3491
3492 ELSE IF (this%trans%sub_type == 'frequency') THEN
3493
3494 ALLOCATE(nval(this%outnx, this%outny))
3495 field_out(:,:,:) = 0.0
3496 DO k = 1, innz
3497 nval(:,:) = 0
3498 DO j = 1, this%inny
3499 DO i = 1, this%innx
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
3506 ENDIF
3507 ENDDO
3508 ENDDO
3509 WHERE (nval(:,:) /= 0)
3510 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3511 ELSEWHERE
3512 field_out(:,:,k) = rmiss
3513 END WHERE
3514 ENDDO
3515 DEALLOCATE(nval)
3516
3517 ENDIF
3518
3519ELSE IF (this%trans%trans_type == 'stencilinter') THEN
3520 np = SIZE(this%stencil,1)/2
3521 ns = SIZE(this%stencil)
3522
3523 IF (this%trans%sub_type == 'average') THEN
3524
3525 IF (vartype == var_dir360) THEN
3526 DO k = 1, innz
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), &
3535 0.0, 360.0, 5.0, &
3536 mask=this%stencil(:,:)) ! simpler and equivalent form
3537! mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3538 ENDIF
3539 ENDDO
3540 ENDDO
3541 ENDDO
3542
3543 ELSE
3544!$OMP PARALLEL DEFAULT(SHARED)
3545!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3546 DO k = 1, innz
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(:,:))
3555 IF (n > 0) THEN
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
3558 ENDIF
3559 ENDIF
3560 ENDDO
3561 ENDDO
3562 ENDDO
3563!$OMP END PARALLEL
3564 ENDIF
3565
3566 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3567 this%trans%sub_type == 'stddevnm1') THEN
3568
3569 IF (this%trans%sub_type == 'stddev') THEN
3570 nm1 = .false.
3571 ELSE
3572 nm1 = .true.
3573 ENDIF
3574
3575!$OMP PARALLEL DEFAULT(SHARED)
3576!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3577 DO k = 1, innz
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
3585! da paura
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)
3590 ENDIF
3591 ENDDO
3592 ENDDO
3593 ENDDO
3594!$OMP END PARALLEL
3595
3596 ELSE IF (this%trans%sub_type == 'max') THEN
3597
3598!$OMP PARALLEL DEFAULT(SHARED)
3599!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3600 DO k = 1, innz
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(:,:))
3609 IF (n > 0) THEN
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(:,:))
3612 ENDIF
3613 ENDIF
3614 ENDDO
3615 ENDDO
3616 ENDDO
3617!$OMP END PARALLEL
3618
3619 ELSE IF (this%trans%sub_type == 'min') THEN
3620
3621!$OMP PARALLEL DEFAULT(SHARED)
3622!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3623 DO k = 1, innz
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(:,:))
3632 IF (n > 0) THEN
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(:,:))
3635 ENDIF
3636 ENDIF
3637 ENDDO
3638 ENDDO
3639 ENDDO
3640!$OMP END PARALLEL
3641
3642 ELSE IF (this%trans%sub_type == 'percentile') THEN
3643
3644!$OMP PARALLEL DEFAULT(SHARED)
3645!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3646 DO k = 1, innz
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
3654! da paura
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/)))
3660 ENDIF
3661 ENDDO
3662 ENDDO
3663 ENDDO
3664!$OMP END PARALLEL
3665
3666 ELSE IF (this%trans%sub_type == 'frequency') THEN
3667
3668!$OMP PARALLEL DEFAULT(SHARED)
3669!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3670 DO k = 1, innz
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(:,:))
3679 IF (n > 0) THEN
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
3683 ENDIF
3684 ENDIF
3685 ENDDO
3686 ENDDO
3687 ENDDO
3688!$OMP END PARALLEL
3689
3690 ENDIF
3691
3692ELSE IF (this%trans%trans_type == 'maskgen') THEN
3693
3694 DO k = 1, innz
3695 WHERE(c_e(this%inter_index_x(:,:)))
3696 field_out(:,:,k) = real(this%inter_index_x(:,:))
3697 END WHERE
3698 ENDDO
3699
3700ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
3701
3702 IF (this%trans%sub_type == 'all') THEN
3703
3704 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3705
3706 ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' &
3707 .OR. this%trans%sub_type == 'mask') THEN
3708
3709 DO k = 1, innz
3710! this is to sparse-points only, so field_out(:,1,k) is acceptable
3711 field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)))
3712 ENDDO
3713
3714 ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
3715 this%trans%sub_type == 'maskinvalid') THEN
3716
3717 DO k = 1, innz
3718 WHERE (this%point_mask(:,:))
3719 field_out(:,:,k) = field_in(:,:,k)
3720 END WHERE
3721 ENDDO
3722
3723 ELSE IF (this%trans%sub_type == 'lemaskinvalid') THEN
3724
3725 DO k = 1, innz
3726 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3727 field_out(:,:,k) = field_in(:,:,k)
3728 ELSEWHERE
3729 field_out(:,:,k) = rmiss
3730 END WHERE
3731 ENDDO
3732
3733 ELSE IF (this%trans%sub_type == 'ltmaskinvalid') THEN
3734
3735 DO k = 1, innz
3736 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3737 field_out(:,:,k) = field_in(:,:,k)
3738 ELSEWHERE
3739 field_out(:,:,k) = rmiss
3740 END WHERE
3741 ENDDO
3742
3743 ELSE IF (this%trans%sub_type == 'gemaskinvalid') THEN
3744
3745 DO k = 1, innz
3746 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3747 field_out(:,:,k) = field_in(:,:,k)
3748 ELSEWHERE
3749 field_out(:,:,k) = rmiss
3750 END WHERE
3751 ENDDO
3752
3753 ELSE IF (this%trans%sub_type == 'gtmaskinvalid') THEN
3754
3755 DO k = 1, innz
3756 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3757 field_out(:,:,k) = field_in(:,:,k)
3758 ELSEWHERE
3759 field_out(:,:,k) = rmiss
3760 END WHERE
3761 ENDDO
3762
3763 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
3764
3765 DO k = 1, innz
3766 WHERE (c_e(field_in(:,:,k)))
3767 field_out(:,:,k) = field_in(:,:,k)
3768 ELSE WHERE
3769 field_out(:,:,k) = this%val1
3770 END WHERE
3771 ENDDO
3772
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
3778 ELSE WHERE
3779 field_out(:,:,:) = field_in(:,:,:)
3780 END WHERE
3781 ELSE IF (c_e(this%val1)) THEN
3782 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3783 field_out(:,:,:) = rmiss
3784 ELSE WHERE
3785 field_out(:,:,:) = field_in(:,:,:)
3786 END WHERE
3787 ELSE IF (c_e(this%val2)) THEN
3788 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3789 field_out(:,:,:) = rmiss
3790 ELSE WHERE
3791 field_out(:,:,:) = field_in(:,:,:)
3792 END WHERE
3793 ENDIF
3794 ENDIF
3795
3796ELSE IF (this%trans%trans_type == 'vertint') THEN
3797
3798 IF (this%trans%sub_type == 'linear') THEN
3799
3800 alloc_coord_3d_in_act = .false.
3801 IF (ASSOCIATED(this%inter_index_z)) THEN
3802
3803 DO k = 1, outnz
3804 IF (c_e(this%inter_index_z(k))) THEN
3805 z1 = real(this%inter_zp(k)) ! weight for k+1
3806 z2 = real(1.0d0 - this%inter_zp(k)) ! weight for k
3807 SELECT CASE(vartype)
3808
3809 CASE(var_dir360)
3810 DO j = 1, inny
3811 DO i = 1, innx
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)
3817 ENDIF
3818 ENDDO
3819 ENDDO
3820
3821 CASE(var_press)
3822 DO j = 1, inny
3823 DO i = 1, innx
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)
3831 ENDIF
3832 ENDDO
3833 ENDDO
3834
3835 CASE default
3836 DO j = 1, inny
3837 DO i = 1, innx
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
3843 ENDIF
3844 ENDDO
3845 ENDDO
3846
3847 END SELECT
3848
3849 ENDIF
3850 ENDDO
3851
3852 ELSE ! use coord_3d_in
3853
3854 IF (ALLOCATED(this%coord_3d_in)) THEN
3855 coord_3d_in_act => this%coord_3d_in
3856#ifdef DEBUG
3857 CALL l4f_category_log(this%category,l4f_debug, &
3858 "Using external vertical coordinate for vertint")
3859#endif
3860 ELSE
3861 IF (PRESENT(coord_3d_in)) THEN
3862#ifdef DEBUG
3863 CALL l4f_category_log(this%category,l4f_debug, &
3864 "Using internal vertical coordinate for vertint")
3865#endif
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)
3872 ELSEWHERE
3873 coord_3d_in_act = rmiss
3874 END WHERE
3875 ELSE
3876 coord_3d_in_act => coord_3d_in
3877 ENDIF
3878 ELSE
3879 CALL l4f_category_log(this%category,l4f_error, &
3880 "Missing internal and external vertical coordinate for vertint")
3881 CALL raise_error()
3882 RETURN
3883 ENDIF
3884 ENDIF
3885
3886 inused = SIZE(coord_3d_in_act,3)
3887 IF (inused < 2) RETURN ! to avoid algorithm failure
3888 kkcache = 1
3889
3890 DO k = 1, outnz
3891 IF (.NOT.c_e(this%vcoord_out(k))) cycle
3892 DO j = 1, inny
3893 DO i = 1, innx
3894 kfound = imiss
3895 DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
3896 kkup = kkcache + kk
3897 kkdown = kkcache - kk + 1
3898
3899 IF (kkdown >= 1) THEN ! search down
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
3904 kkcache = kkdown
3905 kfoundin = kkcache
3906 kfound = kkcache + this%levshift
3907 EXIT ! kk
3908 ENDIF
3909 ENDIF
3910 IF (kkup < inused) THEN ! search up
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
3915 kkcache = kkup
3916 kfoundin = kkcache
3917 kfound = kkcache + this%levshift
3918 EXIT ! kk
3919 ENDIF
3920 ENDIF
3921
3922 ENDDO
3923
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)))
3928 z2 = 1.0 - z1
3929 SELECT CASE(vartype)
3930
3931 CASE(var_dir360)
3932 field_out(i,j,k) = &
3933 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3934 CASE(var_press)
3935 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3936 log(field_in(i,j,kfound+1))*z1)
3937
3938 CASE default
3939 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3940 END SELECT
3941
3942 ENDIF
3943 ELSE
3944 ENDIF
3945 ENDDO ! i
3946 ENDDO ! j
3947 ENDDO ! k
3948 IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
3949 ENDIF
3950
3951 ELSE IF (this%trans%sub_type == 'linearsparse') THEN
3952
3953
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")
3957 RETURN
3958 ENDIF
3959
3960 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3961 DO j = 1, inny
3962 DO i = 1, innx
3963
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(:))
3967 ELSE
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,:))
3970 ENDIF
3971
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)
3975 ENDIF
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)
3980 ELSE
3981 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3982 ENDIF
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), &
3986 mask=mask_in))
3987 ELSE
3988 val_in(1:inused) = pack( &
3989 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3990 mask=mask_in)
3991 ENDIF
3992 kkcache = 1
3993 DO k = 1, outnz
3994
3995 kfound = imiss
3996 DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
3997 kkup = kkcache + kk
3998 kkdown = kkcache - kk + 1
3999
4000 IF (kkdown >= 1) THEN ! search down
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
4005 kkcache = kkdown
4006 kfoundin = kkcache
4007 kfound = kkcache
4008 EXIT ! kk
4009 ENDIF
4010 ENDIF
4011 IF (kkup < inused) THEN ! search up
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
4016 kkcache = kkup
4017 kfoundin = kkcache
4018 kfound = kkcache
4019 EXIT ! kk
4020 ENDIF
4021 ENDIF
4022
4023 ENDDO
4024
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)))
4028 z2 = 1.0 - z1
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)
4034 ELSE
4035 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
4036 ENDIF
4037 ENDIF
4038
4039 ENDDO
4040
4041 ENDIF
4042
4043 ENDDO
4044 ENDDO
4045 DEALLOCATE(coord_in,val_in)
4046
4047
4048 ENDIF
4049
4050ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN
4051
4052 field_out(:,:,:) = field_in(:,:,:)
4053
4054ENDIF
4055
4056
4057CONTAINS
4058
4059
4060! internal function for interpolating directions from 0 to 360 degree
4061! hope it is inlined by the compiler
4062FUNCTION interp_var_360(v1, v2, w1, w2)
4063REAL,INTENT(in) :: v1, v2, w1, w2
4064REAL :: interp_var_360
4065
4066REAL :: lv1, lv2
4067
4068IF (abs(v1 - v2) > 180.) THEN
4069 IF (v1 > v2) THEN
4070 lv1 = v1 - 360.
4071 lv2 = v2
4072 ELSE
4073 lv1 = v1
4074 lv2 = v2 - 360.
4075 ENDIF
4076 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4077ELSE
4078 interp_var_360 = v1*w2 + v2*w1
4079ENDIF
4080
4081END FUNCTION interp_var_360
4082
4083
4084RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
4085REAL,INTENT(in) :: v1(:,:)
4086REAL,INTENT(in) :: l, h, res
4087LOGICAL,INTENT(in),OPTIONAL :: mask(:,:)
4088REAL :: prev
4089
4090REAL :: m
4091INTEGER :: nh, nl
4092!REAL,PARAMETER :: res = 1.0
4093
4094m = (l + h)/2.
4095IF ((h - l) <= res) THEN
4096 prev = m
4097 RETURN
4098ENDIF
4099
4100IF (PRESENT(mask)) THEN
4101 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4102 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4103ELSE
4104 nl = count(v1 >= l .AND. v1 < m)
4105 nh = count(v1 >= m .AND. v1 < h)
4106ENDIF
4107IF (nh > nl) THEN
4108 prev = find_prevailing_direction(v1, m, h, res)
4109ELSE IF (nl > nh) THEN
4110 prev = find_prevailing_direction(v1, l, m, res)
4111ELSE IF (nl == 0 .AND. nh == 0) THEN
4112 prev = rmiss
4113ELSE
4114 prev = m
4115ENDIF
4116
4117END FUNCTION find_prevailing_direction
4118
4119
4120END SUBROUTINE grid_transform_compute
4121
4122
4128SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4129 coord_3d_in)
4130TYPE(grid_transform),INTENT(in) :: this
4131REAL, INTENT(in) :: field_in(:,:)
4132REAL, INTENT(out):: field_out(:,:,:)
4133TYPE(vol7d_var),INTENT(in),OPTIONAL :: var
4134REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:)
4135
4136real,allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4137INTEGER :: inn_p, ier, k
4138INTEGER :: innx, inny, innz, outnx, outny, outnz
4139
4140#ifdef DEBUG
4141call l4f_category_log(this%category,l4f_debug,"start v7d_grid_transform_compute")
4142#endif
4143
4144field_out(:,:,:) = rmiss
4145
4146IF (.NOT.this%valid) THEN
4147 call l4f_category_log(this%category,l4f_error, &
4148 "refusing to perform a non valid transformation")
4149 call raise_error()
4150 RETURN
4151ENDIF
4152
4153innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
4154outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
4155
4156! check size of field_in, field_out
4157IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation
4158
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))
4163 CALL raise_error()
4164 RETURN
4165 ENDIF
4166
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))
4171 CALL raise_error()
4172 RETURN
4173 ENDIF
4174
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))
4180 CALL raise_error()
4181 RETURN
4182 ENDIF
4183
4184ELSE ! horizontal interpolation
4185
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))
4191 CALL raise_error()
4192 RETURN
4193 ENDIF
4194
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))
4200 CALL raise_error()
4201 RETURN
4202 ENDIF
4203
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))
4208 CALL raise_error()
4209 RETURN
4210 ENDIF
4211
4212ENDIF
4213
4214#ifdef DEBUG
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))
4218#endif
4219
4220IF (this%trans%trans_type == 'inter') THEN
4221
4222 IF (this%trans%sub_type == 'linear') THEN
4223
4224#ifdef HAVE_LIBNGMATH
4225! optimization, allocate only once with a safe size
4226 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4227 DO k = 1, innz
4228 inn_p = count(c_e(field_in(:,k)))
4229
4230 CALL l4f_category_log(this%category,l4f_info, &
4231 "Number of sparse data points: "//t2c(inn_p)//','//t2c(SIZE(field_in(:,k))))
4232
4233 IF (inn_p > 2) THEN
4234
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)))
4238
4239 IF (.NOT.this%trans%extrap) THEN
4240 CALL nnseti('ext', 0) ! 0 = no extrapolation
4241 CALL nnsetr('nul', rmiss)
4242 ENDIF
4243
4244 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, & ! (1:inn_p) omitted
4245 this%outnx, this%outny, real(this%inter_x(:,1)), & ! no f90 interface
4246 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4247
4248 IF (ier /= 0) THEN
4249 CALL l4f_category_log(this%category,l4f_error, &
4250 "Error code from NCAR natgrids: "//t2c(ier))
4251 CALL raise_error()
4252 EXIT
4253 ENDIF ! exit loop to deallocate
4254 ELSE
4255
4256 CALL l4f_category_log(this%category,l4f_info, &
4257 "insufficient data in gridded region to triangulate")
4258
4259 ENDIF
4260 ENDDO
4261 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4262
4263#else
4264 CALL l4f_category_log(this%category,l4f_error, &
4265 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4266 CALL raise_error()
4267 RETURN
4268#endif
4269
4270 ENDIF
4271
4272ELSE 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 ! use the grid-to-grid method
4276
4277 CALL compute(this, &
4278 reshape(field_in, (/SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
4279 coord_3d_in)
4280
4281ENDIF
4282
4283END SUBROUTINE grid_transform_v7d_grid_compute
4284
4285
4286! Bilinear interpolation
4287! Effettua interpolazione bilineare dati i valori nei punti 1,2,3,4 e
4288! le coordinate dei punti 1 e 3 oltre a quelle del punto p dove viene
4289! valutato il campo.
4290!_____________________________________________________________
4291! disposizione punti
4292! 4 3
4293!
4294! p
4295!
4296! 1 2
4297! _____________________________________________________________
4298ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
4299REAL,INTENT(in) :: z1,z2,z3,z4 ! Z values on the four points
4300DOUBLE PRECISION,INTENT(in):: x1,y1 ! coordinate of the lower left point
|