libsim  Versione 7.2.4

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

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

Generated with Doxygen.