libsim  Versione7.2.1

◆ grid_transform_compute()

recursive subroutine grid_transform_class::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 3107 del file grid_transform_class.F90.

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

Generated with Doxygen.