libsim  Versione 7.1.9

◆ 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 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 
3186 ELSE 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
3268 ELSE 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 
3411 ELSE 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 
3584 ELSE 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 
3592 ELSE 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 
3688 ELSE 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 
3942 ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN
3943 
3944  field_out(:,:,:) = field_in(:,:,:)
3945 
3946 ENDIF
3947 
3948 
3949 CONTAINS
3950 
3951 
3952 ! internal function for interpolating directions from 0 to 360 degree
3953 ! hope it is inlined by the compiler
3954 FUNCTION interp_var_360(v1, v2, w1, w2)
3955 REAL,INTENT(in) :: v1, v2, w1, w2
3956 REAL :: interp_var_360
3957 
3958 REAL :: lv1, lv2
3959 
3960 IF (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.)
3969 ELSE
3970  interp_var_360 = v1*w2 + v2*w1
3971 ENDIF
3972 
3973 END FUNCTION interp_var_360
3974 
3975 
3976 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
3977 REAL,INTENT(in) :: v1(:,:)
3978 REAL,INTENT(in) :: l, h, res
3979 LOGICAL,INTENT(in),OPTIONAL :: mask(:,:)
3980 REAL :: prev
3981 
3982 REAL :: m
3983 INTEGER :: nh, nl
3984 !REAL,PARAMETER :: res = 1.0
3985 
3986 m = (l + h)/2.
3987 IF ((h - l) <= res) THEN
3988  prev = m
3989  RETURN
3990 ENDIF
3991 
3992 IF (PRESENT(mask)) THEN
3993  nl = count(v1 >= l .AND. v1 < m .AND. mask)
3994  nh = count(v1 >= m .AND. v1 < h .AND. mask)
3995 ELSE
3996  nl = count(v1 >= l .AND. v1 < m)
3997  nh = count(v1 >= m .AND. v1 < h)
3998 ENDIF
3999 IF (nh > nl) THEN
4000  prev = find_prevailing_direction(v1, m, h, res)
4001 ELSE IF (nl > nh) THEN
4002  prev = find_prevailing_direction(v1, l, m, res)
4003 ELSE IF (nl == 0 .AND. nh == 0) THEN
4004  prev = rmiss
4005 ELSE
4006  prev = m
4007 ENDIF
4008 
4009 END FUNCTION find_prevailing_direction
4010 
4011 
4012 END SUBROUTINE grid_transform_compute
4013 
4014 
4020 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4021  coord_3d_in)
4022 TYPE(grid_transform),INTENT(in) :: this
4023 REAL, INTENT(in) :: field_in(:,:)
4024 REAL, INTENT(out):: field_out(:,:,:)
4025 TYPE(vol7d_var),INTENT(in),OPTIONAL :: var
4026 REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:)
4027 
4028 real,allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4029 INTEGER :: inn_p, ier, k
4030 INTEGER :: innx, inny, innz, outnx, outny, outnz
4031 
4032 #ifdef DEBUG
4033 call l4f_category_log(this%category,L4F_DEBUG,"start v7d_grid_transform_compute")
4034 #endif
4035 
4036 field_out(:,:,:) = rmiss
4037 
4038 IF (.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
4043 ENDIF
4044 
4045 innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
4046 outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
4047 
4048 ! check size of field_in, field_out
4049 IF (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 
4076 ELSE ! 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 
4104 ENDIF
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 
4112 IF (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 
4173 ENDIF
4174 
4175 END 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 ! _____________________________________________________________
4190 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
4191 REAL,INTENT(in) :: z1,z2,z3,z4 ! Z values on the four points
4192 DOUBLE PRECISION,INTENT(in):: x1,y1 ! coordinate of the lower left point
4193 DOUBLE PRECISION,INTENT(in):: x3,y3 ! coordinate of the upper right point
4194 DOUBLE PRECISION,INTENT(in):: xp,yp ! coordinate of point where interpolate
4195 REAL :: zp
4196 
4197 REAL :: p1,p2
4198 REAL :: z5,z6

Generated with Doxygen.