libsim  Versione 7.1.7

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

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

Generated with Doxygen.