libsim  Versione7.2.6

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

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

Generated with Doxygen.