libsim Versione 7.1.11

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

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

Generated with Doxygen.