libsim  Versione6.3.0

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

2866  ENDDO
2867 
2868  ELSE IF (this%trans%sub_type == 'max') THEN
2869  DO k = 1, innz
2870  jj = 0
2871  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
2872  je = j+this%trans%box_info%npy-1
2873  jj = jj+1
2874  ii = 0
2875  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
2876  ie = i+this%trans%box_info%npx-1
2877  ii = ii+1
2878  navg = count(field_in(i:ie,j:je,k) /= rmiss)
2879  IF (navg > 0) THEN
2880  field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
2881  mask=(field_in(i:ie,j:je,k) /= rmiss))
2882  ENDIF
2883  ENDDO
2884  ENDDO
2885  ENDDO
2886 
2887  ELSE IF (this%trans%sub_type == 'min') THEN
2888  DO k = 1, innz
2889  jj = 0
2890  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
2891  je = j+this%trans%box_info%npy-1
2892  jj = jj+1
2893  ii = 0
2894  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
2895  ie = i+this%trans%box_info%npx-1
2896  ii = ii+1
2897  navg = count(field_in(i:ie,j:je,k) /= rmiss)
2898  IF (navg > 0) THEN
2899  field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
2900  mask=(field_in(i:ie,j:je,k) /= rmiss))
2901  ENDIF
2902  ENDDO
2903  ENDDO
2904  ENDDO
2905 
2906  ELSE IF (this%trans%sub_type == 'percentile') THEN
2907 
2908  navg = this%trans%box_info%npx*this%trans%box_info%npy
2909  DO k = 1, innz
2910  jj = 0
2911  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
2912  je = j+this%trans%box_info%npy-1
2913  jj = jj+1
2914  ii = 0
2915  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
2916  ie = i+this%trans%box_info%npx-1
2917  ii = ii+1
2918  field_out(ii:ii,jj,k) = stat_percentile( &
2919  reshape(field_in(i:ie,j:je,k),(/navg/)), &
2920  (/REAL(this%trans%stat_info%percentile)/))
2921  ENDDO
2922  ENDDO
2923  ENDDO
2924 
2925  ELSE IF (this%trans%sub_type == 'frequency') THEN
2926 
2927  DO k = 1, innz
2928  jj = 0
2929  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
2930  je = j+this%trans%box_info%npy-1
2931  jj = jj+1
2932  ii = 0
2933  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
2934  ie = i+this%trans%box_info%npx-1
2935  ii = ii+1
2936  navg = count(field_in(i:ie,j:je,k) /= rmiss)
2937  IF (navg > 0) THEN
2938  field_out(ii,jj,k) = sum(interval_info_valid( &
2939  this%trans%interval_info, field_in(i:ie,j:je,k)), &
2940  mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
2941  ENDIF
2942  ENDDO
2943  ENDDO
2944  ENDDO
2945 
2946  ENDIF
2947 
2948 ELSE IF (this%trans%trans_type == 'inter') THEN
2949 
2950  IF (this%trans%sub_type == 'near') THEN
2951 
2952  DO k = 1, innz
2953  DO j = 1, this%outny
2954  DO i = 1, this%outnx
2955 
2956  IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
2957  field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2958 
2959  ENDDO
2960  ENDDO
2961  ENDDO
2962 
2963  ELSE IF (this%trans%sub_type == 'bilin') THEN
2964 
2965  DO k = 1, innz
2966  DO j = 1, this%outny
2967  DO i = 1, this%outnx
2968 
2969  IF (c_e(this%inter_index_x(i,j))) THEN
2970 
2971  z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2972  z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
2973  z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
2974  z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
2975 
2976  IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN
2977 
2978  x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
2979  y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
2980  x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
2981  y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
2982 
2983  xp=this%inter_xp(i,j)
2984  yp=this%inter_yp(i,j)
2985 
2986  field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
2987 
2988  ENDIF
2989  ENDIF
2990 
2991  ENDDO
2992  ENDDO
2993  ENDDO
2994  ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
2995  DO k = 1, innz
2996  DO j = 1, this%outny
2997  DO i = 1, this%outnx
2998 
2999  IF (c_e(this%inter_index_x(i,j))) THEN
3000 
3001  IF(this%inter_index_x(i,j)-1>0)THEN
3002  z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3003  ELSE
3004  z(1) = rmiss
3005  END IF
3006  IF(this%inter_index_x(i,j)+1<=this%outnx)THEN
3007  z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3008  ELSE
3009  z(3) = rmiss
3010  END IF
3011  IF(this%inter_index_y(i,j)+1<=this%outny)THEN
3012  z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3013  ELSE
3014  z(2) = rmiss
3015  END IF
3016  IF(this%inter_index_y(i,j)-1>0)THEN
3017  z(4)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)-1,k)
3018  ELSE
3019  z(4) = rmiss
3020  END IF
3021  field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3022 
3023  ENDIF
3024 
3025  ENDDO
3026  ENDDO
3027  ENDDO
3028 
3029  ENDIF
3030 ELSE IF (this%trans%trans_type == 'boxinter' &
3031  .OR. this%trans%trans_type == 'polyinter' &
3032  .OR. this%trans%trans_type == 'maskinter') THEN
3033 
3034  IF (this%trans%sub_type == 'average') THEN
3035 
3036  IF (vartype == var_dir360) THEN
3037  DO k = 1, innz
3038  DO j = 1, this%outny
3039  DO i = 1, this%outnx
3040  field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3041  0.0, 360.0, 5.0, &
3042  mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3043  ENDDO
3044  ENDDO
3045  ENDDO
3046 
3047  ELSE
3048  ALLOCATE(nval(this%outnx, this%outny))
3049  field_out(:,:,:) = 0.0
3050  DO k = 1, innz
3051  nval(:,:) = 0
3052  DO j = 1, this%inny
3053  DO i = 1, this%innx
3054  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3055  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3056  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3057  field_in(i,j,k)
3058  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3059  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3060  ENDIF
3061  ENDDO
3062  ENDDO
3063  WHERE (nval(:,:) /= 0)
3064  field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3065  ELSEWHERE
3066  field_out(:,:,k) = rmiss
3067  END WHERE
3068  ENDDO
3069  DEALLOCATE(nval)
3070  ENDIF
3071 
3072  ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3073  this%trans%sub_type == 'stddevnm1') THEN
3074 
3075  IF (this%trans%sub_type == 'stddev') THEN
3076  nm1 = .false.
3077  ELSE
3078  nm1 = .true.
3079  ENDIF
3080  DO k = 1, innz
3081  DO j = 1, this%outny
3082  DO i = 1, this%outnx
3083 ! da paura
3084  field_out(i:i,j,k) = stat_stddev( &
3085  reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
3086  mask=reshape((this%inter_index_x == i .AND. &
3087  this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)), nm1=nm1)
3088  ENDDO
3089  ENDDO
3090  ENDDO
3091 
3092  ELSE IF (this%trans%sub_type == 'max') THEN
3093 
3094  DO k = 1, innz
3095  DO j = 1, this%inny
3096  DO i = 1, this%innx
3097  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3098  IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3099  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3100  max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3101  field_in(i,j,k))
3102  ELSE
3103  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3104  field_in(i,j,k)
3105  ENDIF
3106  ENDIF
3107  ENDDO
3108  ENDDO
3109  ENDDO
3110 
3111 
3112  ELSE IF (this%trans%sub_type == 'min') THEN
3113 
3114  DO k = 1, innz
3115  DO j = 1, this%inny
3116  DO i = 1, this%innx
3117  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3118  IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3119  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3120  min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3121  field_in(i,j,k))
3122  ELSE
3123  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3124  field_in(i,j,k)
3125  ENDIF
3126  ENDIF
3127  ENDDO
3128  ENDDO
3129  ENDDO
3130 
3131  ELSE IF (this%trans%sub_type == 'percentile') THEN
3132 
3133  DO k = 1, innz
3134  DO j = 1, this%outny
3135  DO i = 1, this%outnx
3136 ! da paura
3137  field_out(i:i,j,k) = stat_percentile( &
3138  reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
3139  (/REAL(this%trans%stat_info%percentile)/), &
3140  mask=reshape((this%inter_index_x == i .AND. &
3141  this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)))
3142  ENDDO
3143  ENDDO
3144  ENDDO
3145 
3146  ELSE IF (this%trans%sub_type == 'frequency') THEN
3147 
3148  ALLOCATE(nval(this%outnx, this%outny))
3149  field_out(:,:,:) = 0.0
3150  DO k = 1, innz
3151  nval(:,:) = 0
3152  DO j = 1, this%inny
3153  DO i = 1, this%innx
3154  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3155  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3156  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3157  interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3158  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3159  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3160  ENDIF
3161  ENDDO
3162  ENDDO
3163  WHERE (nval(:,:) /= 0)
3164  field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3165  ELSEWHERE
3166  field_out(:,:,k) = rmiss
3167  END WHERE
3168  ENDDO
3169  DEALLOCATE(nval)
3170 
3171  ENDIF
3172 
3173 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
3174  np = SIZE(this%stencil,1)/2
3175  ns = SIZE(this%stencil)
3176 
3177  IF (this%trans%sub_type == 'average') THEN
3178 
3179  IF (vartype == var_dir360) THEN
3180  DO k = 1, innz
3181  DO j = 1, this%outny
3182  DO i = 1, this%outnx
3183  IF (c_e(this%inter_index_x(i,j))) THEN
3184  i1 = this%inter_index_x(i,j) - np
3185  i2 = this%inter_index_x(i,j) + np
3186  j1 = this%inter_index_y(i,j) - np
3187  j2 = this%inter_index_y(i,j) + np
3188  field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3189  0.0, 360.0, 5.0, &
3190  mask=this%stencil(:,:)) ! simpler and equivalent form
3191 ! mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3192  ENDIF
3193  ENDDO
3194  ENDDO
3195  ENDDO
3196 
3197  ELSE
3198 !$OMP PARALLEL DEFAULT(SHARED)
3199 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3200  DO k = 1, innz
3201  DO j = 1, this%outny
3202  DO i = 1, this%outnx
3203  IF (c_e(this%inter_index_x(i,j))) THEN
3204  i1 = this%inter_index_x(i,j) - np
3205  i2 = this%inter_index_x(i,j) + np
3206  j1 = this%inter_index_y(i,j) - np
3207  j2 = this%inter_index_y(i,j) + np
3208  n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3209  IF (n > 0) THEN
3210  field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3211  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3212  ENDIF
3213  ENDIF
3214  ENDDO
3215  ENDDO
3216  ENDDO
3217 !$OMP END PARALLEL
3218  ENDIF
3219 
3220  ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3221  this%trans%sub_type == 'stddevnm1') THEN
3222 
3223  IF (this%trans%sub_type == 'stddev') THEN
3224  nm1 = .false.
3225  ELSE
3226  nm1 = .true.
3227  ENDIF
3228 
3229 !$OMP PARALLEL DEFAULT(SHARED)
3230 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3231  DO k = 1, innz
3232  DO j = 1, this%outny
3233  DO i = 1, this%outnx
3234  IF (c_e(this%inter_index_x(i,j))) THEN
3235  i1 = this%inter_index_x(i,j) - np
3236  i2 = this%inter_index_x(i,j) + np
3237  j1 = this%inter_index_y(i,j) - np
3238  j2 = this%inter_index_y(i,j) + np
3239 ! da paura
3240  field_out(i:i,j,k) = stat_stddev( &
3241  reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3242  mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3243  this%stencil(:,:), (/ns/)), nm1=nm1)
3244  ENDIF
3245  ENDDO
3246  ENDDO
3247  ENDDO
3248 !$OMP END PARALLEL
3249 
3250  ELSE IF (this%trans%sub_type == 'max') THEN
3251 
3252 !$OMP PARALLEL DEFAULT(SHARED)
3253 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3254  DO k = 1, innz
3255  DO j = 1, this%outny
3256  DO i = 1, this%outnx
3257  IF (c_e(this%inter_index_x(i,j))) THEN
3258  i1 = this%inter_index_x(i,j) - np
3259  i2 = this%inter_index_x(i,j) + np
3260  j1 = this%inter_index_y(i,j) - np
3261  j2 = this%inter_index_y(i,j) + np
3262  n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3263  IF (n > 0) THEN
3264  field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3265  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3266  ENDIF
3267  ENDIF
3268  ENDDO
3269  ENDDO
3270  ENDDO
3271 !$OMP END PARALLEL
3272 
3273  ELSE IF (this%trans%sub_type == 'min') THEN
3274 
3275 !$OMP PARALLEL DEFAULT(SHARED)
3276 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3277  DO k = 1, innz
3278  DO j = 1, this%outny
3279  DO i = 1, this%outnx
3280  IF (c_e(this%inter_index_x(i,j))) THEN
3281  i1 = this%inter_index_x(i,j) - np
3282  i2 = this%inter_index_x(i,j) + np
3283  j1 = this%inter_index_y(i,j) - np
3284  j2 = this%inter_index_y(i,j) + np
3285  n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3286  IF (n > 0) THEN
3287  field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3288  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3289  ENDIF
3290  ENDIF
3291  ENDDO
3292  ENDDO
3293  ENDDO
3294 !$OMP END PARALLEL
3295 
3296  ELSE IF (this%trans%sub_type == 'percentile') THEN
3297 
3298 !$OMP PARALLEL DEFAULT(SHARED)
3299 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3300  DO k = 1, innz
3301  DO j = 1, this%outny
3302  DO i = 1, this%outnx
3303  IF (c_e(this%inter_index_x(i,j))) THEN
3304  i1 = this%inter_index_x(i,j) - np
3305  i2 = this%inter_index_x(i,j) + np
3306  j1 = this%inter_index_y(i,j) - np
3307  j2 = this%inter_index_y(i,j) + np
3308 ! da paura
3309  field_out(i:i,j,k) = stat_percentile( &
3310  reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3311  (/REAL(this%trans%stat_info%percentile)/), &
3312  mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3313  this%stencil(:,:), (/ns/)))
3314  ENDIF
3315  ENDDO
3316  ENDDO
3317  ENDDO
3318 !$OMP END PARALLEL
3319 
3320  ELSE IF (this%trans%sub_type == 'frequency') THEN
3321 
3322 !$OMP PARALLEL DEFAULT(SHARED)
3323 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3324  DO k = 1, innz
3325  DO j = 1, this%outny
3326  DO i = 1, this%outnx
3327  IF (c_e(this%inter_index_x(i,j))) THEN
3328  i1 = this%inter_index_x(i,j) - np
3329  i2 = this%inter_index_x(i,j) + np
3330  j1 = this%inter_index_y(i,j) - np
3331  j2 = this%inter_index_y(i,j) + np
3332  n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3333  IF (n > 0) THEN
3334  field_out(i,j,k) = sum(interval_info_valid( &
3335  this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3336  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3337  ENDIF
3338  ENDIF
3339  ENDDO
3340  ENDDO
3341  ENDDO
3342 !$OMP END PARALLEL
3343 
3344  ENDIF
3345 
3346 ELSE IF (this%trans%trans_type == 'maskgen') THEN
3347 
3348  DO k = 1, innz
3349  WHERE(c_e(this%inter_index_x(:,:)))
3350  field_out(:,:,k) = REAL(this%inter_index_x(:,:))
3351  END WHERE
3352  ENDDO
3353 
3354 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
3355 
3356  IF (this%trans%sub_type == 'all') THEN
3357 
3358  field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3359 
3360  ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' &
3361  .OR. this%trans%sub_type == 'mask') THEN
3362 
3363  DO k = 1, innz
3364 ! this is to sparse-points only, so field_out(:,1,k) is acceptable
3365  field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)))
3366  ENDDO
3367 
3368  ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
3369  this%trans%sub_type == 'maskinvalid') THEN
3370 
3371  DO k = 1, innz
3372  WHERE (this%point_mask(:,:))
3373  field_out(:,:,k) = field_in(:,:,k)
3374  END WHERE
3375  ENDDO
3376 
3377  ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
3378 
3379  DO k = 1, innz
3380  WHERE (c_e(field_in(:,:,k)))
3381  field_out(:,:,k) = field_in(:,:,k)
3382  ELSE WHERE
3383  field_out(:,:,k) = this%val1
3384  END WHERE
3385  ENDDO
3386 
3387  ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
3388  IF (c_e(this%val1) .AND. c_e(this%val2)) THEN
3389  WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3390  .AND. field_in(:,:,:) <= this%val2)
3391  field_out(:,:,:) = rmiss
3392  ELSE WHERE
3393  field_out(:,:,:) = field_in(:,:,:)
3394  END WHERE
3395  ELSE IF (c_e(this%val1)) THEN
3396  WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3397  field_out(:,:,:) = rmiss
3398  ELSE WHERE
3399  field_out(:,:,:) = field_in(:,:,:)
3400  END WHERE
3401  ELSE IF (c_e(this%val2)) THEN
3402  WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3403  field_out(:,:,:) = rmiss
3404  ELSE WHERE
3405  field_out(:,:,:) = field_in(:,:,:)
3406  END WHERE
3407  ENDIF
3408  ENDIF
3409 
3410 ELSE IF (this%trans%trans_type == 'vertint') THEN
3411 
3412  IF (this%trans%sub_type == 'linear') THEN
3413 
3414  alloc_coord_3d_in_act = .false.
3415  IF (ASSOCIATED(this%inter_index_z)) THEN
3416 
3417  DO k = 1, outnz
3418  IF (c_e(this%inter_index_z(k))) THEN
3419  z1 = REAL(this%inter_zp(k)) ! weight for k+1
3420  z2 = REAL(1.0D0 - this%inter_zp(k)) ! weight for k
3421  SELECT CASE(vartype)
3422 
3423  CASE(var_dir360)
3424  DO j = 1, inny
3425  DO i = 1, innx
3426  IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3427  c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3428  field_out(i,j,k) = &
3429  interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3430  field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3431  ENDIF
3432  ENDDO
3433  ENDDO
3434 
3435  CASE(var_press)
3436  DO j = 1, inny
3437  DO i = 1, innx
3438  IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3439  c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3440  field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3441  field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN
3442  field_out(i,j,k) = exp( &
3443  log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3444  log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3445  ENDIF
3446  ENDDO
3447  ENDDO
3448 
3449  CASE default
3450  DO j = 1, inny
3451  DO i = 1, innx
3452  IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3453  c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3454  field_out(i,j,k) = &
3455  field_in(i,j,this%inter_index_z(k))*z2 + &
3456  field_in(i,j,this%inter_index_z(k)+1)*z1
3457  ENDIF
3458  ENDDO
3459  ENDDO
3460 
3461  END SELECT
3462 
3463  ENDIF
3464  ENDDO
3465 
3466  ELSE ! use coord_3d_in
3467 
3468  IF (ALLOCATED(this%coord_3d_in)) THEN
3469  coord_3d_in_act => this%coord_3d_in
3470 #ifdef DEBUG
3471  CALL l4f_category_log(this%category,l4f_debug, &
3472  "Using external vertical coordinate for vertint")
3473 #endif
3474  ELSE
3475  IF (PRESENT(coord_3d_in)) THEN
3476 #ifdef DEBUG
3477  CALL l4f_category_log(this%category,l4f_debug, &
3478  "Using internal vertical coordinate for vertint")
3479 #endif
3480  IF (this%dolog) THEN
3481  ALLOCATE(coord_3d_in_act(SIZE(coord_3d_in,1), &
3482  SIZE(coord_3d_in,2), SIZE(coord_3d_in,3)))
3483  alloc_coord_3d_in_act = .true.
3484  WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3485  coord_3d_in_act = log(coord_3d_in)
3486  ELSEWHERE
3487  coord_3d_in_act = rmiss
3488  END WHERE
3489  ELSE
3490  coord_3d_in_act => coord_3d_in
3491  ENDIF
3492  ELSE
3493  CALL l4f_category_log(this%category,l4f_error, &
3494  "Missing internal and external vertical coordinate for vertint")
3495  CALL raise_error()
3496  RETURN
3497  ENDIF
3498  ENDIF
3499 
3500  inused = SIZE(coord_3d_in_act,3)
3501  IF (inused < 2) RETURN ! to avoid algorithm failure
3502  kkcache = 1
3503 
3504  DO k = 1, outnz
3505  IF (.NOT.c_e(this%vcoord_out(k))) cycle
3506  DO j = 1, inny
3507  DO i = 1, innx
3508  kfound = imiss
3509  DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
3510  kkup = kkcache + kk
3511  kkdown = kkcache - kk + 1
3512 
3513  IF (kkdown >= 1) THEN ! search down
3514  IF (this%vcoord_out(k) >= &
3515  min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3516  this%vcoord_out(k) < &
3517  max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1))) THEN
3518  kkcache = kkdown
3519  kfoundin = kkcache
3520  kfound = kkcache + this%levshift
3521  EXIT ! kk
3522  ENDIF
3523  ENDIF
3524  IF (kkup < inused) THEN ! search up
3525  IF (this%vcoord_out(k) >= &
3526  min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3527  this%vcoord_out(k) < &
3528  max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1))) THEN
3529  kkcache = kkup
3530  kfoundin = kkcache
3531  kfound = kkcache + this%levshift
3532  EXIT ! kk
3533  ENDIF
3534  ENDIF
3535 
3536  ENDDO
3537 
3538  IF (c_e(kfound)) THEN
3539  IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound+1))) THEN
3540  z1 = REAL((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3541  (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3542  z2 = 1.0 - z1
3543  SELECT CASE(vartype)
3544 
3545  CASE(var_dir360)
3546  field_out(i,j,k) = &
3547  interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3548  CASE(var_press)
3549  field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3550  log(field_in(i,j,kfound+1))*z1)
3551 
3552  CASE default
3553  field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3554  END SELECT
3555 
3556  ENDIF
3557  ELSE
3558  ENDIF
3559  ENDDO ! i
3560  ENDDO ! j
3561  ENDDO ! k
3562  IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
3563  ENDIF
3564 
3565  ELSE IF (this%trans%sub_type == 'linearsparse') THEN
3566 
3567 
3568  IF (.NOT.ASSOCIATED(this%vcoord_in) .AND. .NOT.PRESENT(coord_3d_in)) THEN
3569  CALL l4f_category_log(this%category,l4f_error, &
3570  "linearsparse interpolation, no input vert coord available")
3571  RETURN
3572  ENDIF
3573 
3574  ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3575  DO j = 1, inny
3576  DO i = 1, innx
3577 
3578  IF (ASSOCIATED(this%vcoord_in)) THEN
3579  mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3580  .AND. c_e(this%vcoord_in(:))
3581  ELSE
3582  mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3583  .AND. c_e(coord_3d_in(i,j,:))
3584  ENDIF
3585 
3586  IF (vartype == var_press) THEN
3587  mask_in(:) = mask_in(:) .AND. &
3588  (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3589  ENDIF
3590  inused = count(mask_in)
3591  IF (inused > 1) THEN
3592  IF (ASSOCIATED(this%vcoord_in)) THEN
3593  coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3594  ELSE
3595  coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3596  ENDIF
3597  IF (vartype == var_press) THEN
3598  val_in(1:inused) = log(pack( &
3599  field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3600  mask=mask_in))
3601  ELSE
3602  val_in(1:inused) = pack( &
3603  field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3604  mask=mask_in)
3605  ENDIF
3606  kkcache = 1
3607  DO k = 1, outnz
3608 
3609  kfound = imiss
3610  DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
3611  kkup = kkcache + kk
3612  kkdown = kkcache - kk + 1
3613 
3614  IF (kkdown >= 1) THEN ! search down
3615  IF (this%vcoord_out(k) >= &
3616  min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
3617  this%vcoord_out(k) < &
3618  max(coord_in(kkdown), coord_in(kkdown+1))) THEN
3619  kkcache = kkdown
3620  kfoundin = kkcache
3621  kfound = kkcache
3622  EXIT ! kk
3623  ENDIF
3624  ENDIF
3625  IF (kkup < inused) THEN ! search up
3626  IF (this%vcoord_out(k) >= &
3627  min(coord_in(kkup), coord_in(kkup+1)) .AND. &
3628  this%vcoord_out(k) < &
3629  max(coord_in(kkup), coord_in(kkup+1))) THEN
3630  kkcache = kkup
3631  kfoundin = kkcache
3632  kfound = kkcache
3633  EXIT ! kk
3634  ENDIF
3635  ENDIF
3636 
3637  ENDDO
3638 
3639  IF (c_e(kfound)) THEN
3640  z1 = REAL((this%vcoord_out(k) - coord_in(kfound-1))/ &
3641  (coord_in(kfound) - coord_in(kfound-1)))
3642  z2 = 1.0 - z1
3643  IF (vartype == var_dir360) THEN
3644  field_out(i,j,k) = &
3645  interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
3646  ELSE IF (vartype == var_press) THEN
3647  field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
3648  ELSE
3649  field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
3650  ENDIF
3651  ENDIF
3652 
3653  ENDDO
3654 
3655  ENDIF
3656 
3657  ENDDO
3658  ENDDO
3659  DEALLOCATE(coord_in,val_in)
3660 
3661 
3662  ENDIF
3663 
3664 ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN
3665 
3666  field_out(:,:,:) = field_in(:,:,:)
3667 
3668 ENDIF
3669 
3670 
3671 CONTAINS
3672 
3673 
3674 ! internal function for interpolating directions from 0 to 360 degree
3675 ! hope it is inlined by the compiler
3676 FUNCTION interp_var_360(v1, v2, w1, w2)
3677 REAL,INTENT(in) :: v1, v2, w1, w2
3678 REAL :: interp_var_360
3679 
3680 REAL :: lv1, lv2
3681 
3682 IF (abs(v1 - v2) > 180.) THEN
3683  IF (v1 > v2) THEN
3684  lv1 = v1 - 360.
3685  lv2 = v2
3686  ELSE
3687  lv1 = v1
3688  lv2 = v2 - 360.
3689  ENDIF
3690  interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
3691 ELSE
3692  interp_var_360 = v1*w2 + v2*w1
3693 ENDIF
3694 
3695 END FUNCTION interp_var_360
3696 
3697 
3698 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
3699 REAL,INTENT(in) :: v1(:,:)
3700 REAL,INTENT(in) :: l, h, res
3701 LOGICAL,INTENT(in),OPTIONAL :: mask(:,:)
3702 REAL :: prev
3703 
3704 REAL :: m
3705 INTEGER :: nh, nl
3706 !REAL,PARAMETER :: res = 1.0
3707 
3708 m = (l + h)/2.
3709 IF ((h - l) <= res) THEN
3710  prev = m
3711  RETURN
3712 ENDIF
3713 
3714 IF (PRESENT(mask)) THEN
3715  nl = count(v1 >= l .AND. v1 < m .AND. mask)
3716  nh = count(v1 >= m .AND. v1 < h .AND. mask)
3717 ELSE
3718  nl = count(v1 >= l .AND. v1 < m)
3719  nh = count(v1 >= m .AND. v1 < h)
3720 ENDIF
3721 IF (nh > nl) THEN
3722  prev = find_prevailing_direction(v1, m, h, res)
3723 ELSE IF (nl > nh) THEN
3724  prev = find_prevailing_direction(v1, l, m, res)
3725 ELSE IF (nl == 0 .AND. nh == 0) THEN
3726  prev = rmiss
3727 ELSE
3728  prev = m
3729 ENDIF
3730 
3731 END FUNCTION find_prevailing_direction
3732 
3733 
3734 END SUBROUTINE grid_transform_compute
3735 
3736 
3742 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
3743  coord_3d_in)
3744 TYPE(grid_transform),INTENT(in) :: this
3745 REAL, INTENT(in) :: field_in(:,:)
3746 REAL, INTENT(out):: field_out(:,:,:)
3747 TYPE(vol7d_var),INTENT(in),OPTIONAL :: var
3748 REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:)
3749 
3750 real,allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
3751 INTEGER :: inn_p, ier, k
3752 INTEGER :: innx, inny, innz, outnx, outny, outnz
3753 
3754 #ifdef DEBUG
3755 call l4f_category_log(this%category,l4f_debug,"start v7d_grid_transform_compute")
3756 #endif
3757 
3758 field_out(:,:,:) = rmiss
3759 
3760 IF (.NOT.this%valid) THEN
3761  call l4f_category_log(this%category,l4f_error, &
3762  "refusing to perform a non valid transformation")
3763  call raise_error()
3764  RETURN
3765 ENDIF
3766 
3767 innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
3768 outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
3769 
3770 ! check size of field_in, field_out
3771 IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation
3772 
3773  IF (innz /= this%innz) THEN
3774  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
3775  CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
3776  t2c(this%innz)//" /= "//t2c(innz))
3777  CALL raise_error()
3778  RETURN
3779  ENDIF
3780 
3781  IF (outnz /= this%outnz) THEN
3782  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
3783  CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
3784  t2c(this%outnz)//" /= "//t2c(outnz))
3785  CALL raise_error()
3786  RETURN
3787  ENDIF
3788 
3789  IF (innx /= outnx .OR. inny /= outny) THEN
3790  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
3791  CALL l4f_category_log(this%category,l4f_error,"inconsistent hor. sizes: "//&
3792  t2c(innx)//","//t2c(inny)//" /= "//&
3793  t2c(outnx)//","//t2c(outny))
3794  CALL raise_error()
3795  RETURN
3796  ENDIF
3797 
3798 ELSE ! horizontal interpolation
3799 
3800  IF (innx /= this%innx .OR. inny /= this%inny) THEN
3801  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
3802  CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
3803  t2c(this%innx)//","//t2c(this%inny)//" /= "//&
3804  t2c(innx)//","//t2c(inny))
3805  CALL raise_error()
3806  RETURN
3807  ENDIF
3808 
3809  IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
3810  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
3811  CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
3812  t2c(this%outnx)//","//t2c(this%outny)//" /= "//&
3813  t2c(outnx)//","//t2c(outny))
3814  CALL raise_error()
3815  RETURN
3816  ENDIF
3817 
3818  IF (innz /= outnz) THEN
3819  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
3820  CALL l4f_category_log(this%category,l4f_error,"inconsistent vert. sizes: "//&
3821  t2c(innz)//" /= "//t2c(outnz))
3822  CALL raise_error()
3823  RETURN
3824  ENDIF
3825 
3826 ENDIF
3827 
3828 #ifdef DEBUG
3829  call l4f_category_log(this%category,l4f_debug, &
3830  "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//':'// &
3831  trim(this%trans%sub_type))
3832 #endif
3833 
3834 IF (this%trans%trans_type == 'inter') THEN
3835 
3836  IF (this%trans%sub_type == 'linear') THEN
3837 
3838 #ifdef HAVE_LIBNGMATH
3839 ! optimization, allocate only once with a safe size
3840  ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
3841  DO k = 1, innz
3842  inn_p = count(c_e(field_in(:,k)))
3843 
3844  CALL l4f_category_log(this%category,l4f_info, &
3845  "Number of sparse data points: "//t2c(inn_p)//','//t2c(SIZE(field_in(:,k))))
3846 
3847  IF (inn_p > 2) THEN
3848 
3849  field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k)))
3850  x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k)))
3851  y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k)))
3852 
3853  IF (.NOT.this%trans%extrap) THEN
3854  CALL nnseti('ext', 0) ! 0 = no extrapolation
3855  CALL nnsetr('nul', rmiss)
3856  ENDIF
3857 
3858  CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, & ! (1:inn_p) omitted
3859  this%outnx, this%outny, REAL(this%inter_x(:,1)), & ! no f90 interface
3860  REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
3861 
3862  IF (ier /= 0) THEN
3863  CALL l4f_category_log(this%category,l4f_error, &
3864  "Error code from NCAR natgrids: "//t2c(ier))
3865  CALL raise_error()
3866  EXIT
3867  ENDIF ! exit loop to deallocate
3868  ELSE
3869 
3870  CALL l4f_category_log(this%category,l4f_info, &
3871  "insufficient data in gridded region to triangulate")
3872 
3873  ENDIF
3874  ENDDO
3875  DEALLOCATE(field_in_p, x_in_p, y_in_p)
3876 
3877 #else
3878  CALL l4f_category_log(this%category,l4f_error, &
3879  "libsim compiled without NATGRIDD (ngmath ncarg library)")
3880  CALL raise_error()
3881  RETURN
3882 #endif
3883 
3884  ENDIF
3885 
3886 ELSE IF (this%trans%trans_type == 'boxinter' .OR. &
3887  this%trans%trans_type == 'polyinter' .OR. &
3888  this%trans%trans_type == 'vertint' .OR. &
3889  this%trans%trans_type == 'metamorphosis') THEN ! use the grid-to-grid method
3890 
3891  CALL compute(this, &
3892  reshape(field_in, (/SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
3893  coord_3d_in)
3894 
3895 ENDIF
3896 
3897 END SUBROUTINE grid_transform_v7d_grid_compute
3898 
3899 
3900 ! Bilinear interpolation
3901 ! Effettua interpolazione bilineare dati i valori nei punti 1,2,3,4 e
3902 ! le coordinate dei punti 1 e 3 oltre a quelle del punto p dove viene
3903 ! valutato il campo.
3904 !_____________________________________________________________
3905 ! disposizione punti
3906 ! 4 3
3907 !
3908 ! p
3909 !
3910 ! 1 2
3911 ! _____________________________________________________________
3912 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
3913 REAL,INTENT(in) :: z1,z2,z3,z4 ! Z values on the four points
3914 DOUBLE PRECISION,INTENT(in):: x1,y1 ! coordinate of the lower left point
3915 DOUBLE PRECISION,INTENT(in):: x3,y3 ! coordinate of the upper right point
3916 DOUBLE PRECISION,INTENT(in):: xp,yp ! coordinate of point where interpolate
3917 REAL :: zp
3918 
3919 REAL :: p1,p2
3920 REAL :: z5,z6
3921 
3922 
3923 p2 = REAL((yp-y1)/(y3-y1))
3924 p1 = REAL((xp-x1)/(x3-x1))
3925 
3926 z5 = (z4-z1)*p2+z1
3927 z6 = (z3-z2)*p2+z2
3928 
3929 zp = (z6-z5)*(p1)+z5
3930 
3931 END FUNCTION hbilin
3932 
3933 
3934 ! Shapiro filter of order 2
3935 FUNCTION shapiro (z,zp) RESULT(zs)
3936 !! z_smoothed(i,j) = z(i,j) * (1-S) + S * sum(z_vicini)/N
3937 !! = z(i,j) - S/N (N*z(i,j) - sum(z_vicini))
3938 REAL,INTENT(in) :: z(:) ! Z values on the four surrounding points
3939 REAL,INTENT(in) :: zp ! Z values on the central point
Functions that return a trimmed CHARACTER representation of the input variable.
Operatore di valore assoluto di un intervallo.

Generated with Doxygen.