|
◆ 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] | this | grid_transformation object |
[in] | field_in | input array |
[out] | field_out | output array |
[in] | var | physical variable to be interpolated, if provided, some ad-hoc algorithms may be used where possible |
[in] | coord_3d_in | input vertical coordinate for vertical interpolation, if not provided by other means |
Definizione alla linea 2866 del file grid_transform_class.F90.
2868 ELSE IF (this%trans%sub_type == 'max') THEN 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 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 2878 navg = count(field_in(i:ie,j:je,k) /= rmiss) 2880 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), & 2881 mask=(field_in(i:ie,j:je,k) /= rmiss)) 2887 ELSE IF (this%trans%sub_type == 'min') THEN 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 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 2897 navg = count(field_in(i:ie,j:je,k) /= rmiss) 2899 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), & 2900 mask=(field_in(i:ie,j:je,k) /= rmiss)) 2906 ELSE IF (this%trans%sub_type == 'percentile') THEN 2908 navg = this%trans%box_info%npx*this%trans%box_info%npy 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 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 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)/)) 2925 ELSE IF (this%trans%sub_type == 'frequency') THEN 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 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 2936 navg = count(field_in(i:ie,j:je,k) /= rmiss) 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 2948 ELSE IF (this%trans%trans_type == 'inter') THEN 2950 IF (this%trans%sub_type == 'near') THEN 2953 DO j = 1, this%outny 2954 DO i = 1, this%outnx 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) 2963 ELSE IF (this%trans%sub_type == 'bilin') THEN 2966 DO j = 1, this%outny 2967 DO i = 1, this%outnx 2969 IF ( c_e(this%inter_index_x(i,j))) THEN 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) 2976 IF ( c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN 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) 2983 xp=this%inter_xp(i,j) 2984 yp=this%inter_yp(i,j) 2986 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) 2994 ELSE IF (this%trans%sub_type == 'shapiro_near') THEN 2996 DO j = 1, this%outny 2997 DO i = 1, this%outnx 2999 IF ( c_e(this%inter_index_x(i,j))) THEN 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) 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) 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) 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) 3021 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)) 3030 ELSE IF (this%trans%trans_type == 'boxinter' & 3031 .OR. this%trans%trans_type == 'polyinter' & 3032 .OR. this%trans%trans_type == 'maskinter') THEN 3034 IF (this%trans%sub_type == 'average') THEN 3036 IF (vartype == var_dir360) THEN 3038 DO j = 1, this%outny 3039 DO i = 1, this%outnx 3040 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), & 3042 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j) 3048 ALLOCATE(nval(this%outnx, this%outny)) 3049 field_out(:,:,:) = 0.0 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) + & 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 3063 WHERE (nval(:,:) /= 0) 3064 field_out(:,:,k) = field_out(:,:,k)/nval(:,:) 3066 field_out(:,:,k) = rmiss 3072 ELSE IF (this%trans%sub_type == 'stddev' .OR. & 3073 this%trans%sub_type == 'stddevnm1') THEN 3075 IF (this%trans%sub_type == 'stddev') THEN 3081 DO j = 1, this%outny 3082 DO i = 1, this%outnx 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) 3092 ELSE IF (this%trans%sub_type == 'max') THEN 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), & 3103 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = & 3112 ELSE IF (this%trans%sub_type == 'min') THEN 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), & 3123 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = & 3131 ELSE IF (this%trans%sub_type == 'percentile') THEN 3134 DO j = 1, this%outny 3135 DO i = 1, this%outnx 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))/))) 3146 ELSE IF (this%trans%sub_type == 'frequency') THEN 3148 ALLOCATE(nval(this%outnx, this%outny)) 3149 field_out(:,:,:) = 0.0 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 3163 WHERE (nval(:,:) /= 0) 3164 field_out(:,:,k) = field_out(:,:,k)/nval(:,:) 3166 field_out(:,:,k) = rmiss 3173 ELSE IF (this%trans%trans_type == 'stencilinter') THEN 3174 np = SIZE(this%stencil,1)/2 3175 ns = SIZE(this%stencil) 3177 IF (this%trans%sub_type == 'average') THEN 3179 IF (vartype == var_dir360) THEN 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), & 3190 mask=this%stencil(:,:)) 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(:,:)) 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 3220 ELSE IF (this%trans%sub_type == 'stddev' .OR. & 3221 this%trans%sub_type == 'stddevnm1') THEN 3223 IF (this%trans%sub_type == 'stddev') THEN 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 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) 3250 ELSE IF (this%trans%sub_type == 'max') THEN 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(:,:)) 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(:,:)) 3273 ELSE IF (this%trans%sub_type == 'min') THEN 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(:,:)) 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(:,:)) 3296 ELSE IF (this%trans%sub_type == 'percentile') THEN 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 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/))) 3320 ELSE IF (this%trans%sub_type == 'frequency') THEN 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(:,:)) 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 3346 ELSE IF (this%trans%trans_type == 'maskgen') THEN 3349 WHERE( c_e(this%inter_index_x(:,:))) 3350 field_out(:,:,k) = REAL(this%inter_index_x(:,:)) 3354 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN 3356 IF (this%trans%sub_type == 'all') THEN 3358 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/)) 3360 ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' & 3361 .OR. this%trans%sub_type == 'mask') THEN 3365 field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:))) 3368 ELSE IF (this%trans%sub_type == 'maskvalid' .OR. & 3369 this%trans%sub_type == 'maskinvalid') THEN 3372 WHERE (this%point_mask(:,:)) 3373 field_out(:,:,k) = field_in(:,:,k) 3377 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN 3380 WHERE ( c_e(field_in(:,:,k))) 3381 field_out(:,:,k) = field_in(:,:,k) 3383 field_out(:,:,k) = this%val1 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 3393 field_out(:,:,:) = field_in(:,:,:) 3395 ELSE IF ( c_e(this%val1)) THEN 3396 WHERE ( c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1) 3397 field_out(:,:,:) = rmiss 3399 field_out(:,:,:) = field_in(:,:,:) 3401 ELSE IF ( c_e(this%val2)) THEN 3402 WHERE ( c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2) 3403 field_out(:,:,:) = rmiss 3405 field_out(:,:,:) = field_in(:,:,:) 3410 ELSE IF (this%trans%trans_type == 'vertint') THEN 3412 IF (this%trans%sub_type == 'linear') THEN 3414 alloc_coord_3d_in_act = .false. 3415 IF ( ASSOCIATED(this%inter_index_z)) THEN 3418 IF ( c_e(this%inter_index_z(k))) THEN 3419 z1 = REAL(this%inter_zp(k)) 3420 z2 = REAL(1.0D0 - this%inter_zp(k)) 3421 SELECT CASE(vartype) 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) 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) 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 3468 IF ( ALLOCATED(this%coord_3d_in)) THEN 3469 coord_3d_in_act => this%coord_3d_in 3471 CALL l4f_category_log(this%category,l4f_debug, & 3472 "Using external vertical coordinate for vertint") 3475 IF ( PRESENT(coord_3d_in)) THEN 3477 CALL l4f_category_log(this%category,l4f_debug, & 3478 "Using internal vertical coordinate for vertint") 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) 3487 coord_3d_in_act = rmiss 3490 coord_3d_in_act => coord_3d_in 3493 CALL l4f_category_log(this%category,l4f_error, & 3494 "Missing internal and external vertical coordinate for vertint") 3500 inused = SIZE(coord_3d_in_act,3) 3501 IF (inused < 2) RETURN 3505 IF (.NOT. c_e(this%vcoord_out(k))) cycle 3509 DO kk = 1, max(inused-kkcache-1, kkcache) 3511 kkdown = kkcache - kk + 1 3513 IF (kkdown >= 1) THEN 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 3520 kfound = kkcache + this%levshift 3524 IF (kkup < inused) THEN 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 3531 kfound = kkcache + this%levshift 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))) 3543 SELECT CASE(vartype) 3546 field_out(i,j,k) = & 3547 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2) 3549 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + & 3550 log(field_in(i,j,kfound+1))*z1) 3553 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1 3562 IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act) 3565 ELSE IF (this%trans%sub_type == 'linearsparse') THEN 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") 3574 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz)) 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(:)) 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,:)) 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) 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) 3595 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in) 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), & 3602 val_in(1:inused) = pack( & 3603 field_in(i,j,this%levshift+1:this%levshift+this%levused), & 3610 DO kk = 1, max(inused-kkcache-1, kkcache) 3612 kkdown = kkcache - kk + 1 3614 IF (kkdown >= 1) THEN 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 3625 IF (kkup < inused) THEN 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 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))) 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) 3649 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1 3659 DEALLOCATE(coord_in,val_in) 3664 ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN 3666 field_out(:,:,:) = field_in(:,:,:) 3676 FUNCTION interp_var_360(v1, v2, w1, w2) 3677 REAL, INTENT(in) :: v1, v2, w1, w2 3678 REAL :: interp_var_360 3682 IF ( abs(v1 - v2) > 180.) THEN 3690 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.) 3692 interp_var_360 = v1*w2 + v2*w1 3695 END FUNCTION interp_var_360 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(:,:) 3709 IF ((h - l) <= res) THEN 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) 3718 nl = count(v1 >= l .AND. v1 < m) 3719 nh = count(v1 >= m .AND. v1 < h) 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 3731 END FUNCTION find_prevailing_direction 3734 END SUBROUTINE grid_transform_compute 3742 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, & 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(:,:,:) 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 3755 call l4f_category_log(this%category,l4f_debug, "start v7d_grid_transform_compute") 3758 field_out(:,:,:) = rmiss 3760 IF (.NOT.this%valid) THEN 3761 call l4f_category_log(this%category,l4f_error, & 3762 "refusing to perform a non valid transformation") 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) 3771 IF (this%trans%trans_type == 'vertint') THEN 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)) 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)) 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)) 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)) 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)) 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)) 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)) 3834 IF (this%trans%trans_type == 'inter') THEN 3836 IF (this%trans%sub_type == 'linear') THEN 3838 #ifdef HAVE_LIBNGMATH 3840 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny)) 3842 inn_p = count( c_e(field_in(:,k))) 3844 CALL l4f_category_log(this%category,l4f_info, & 3845 "Number of sparse data points: "// t2c(inn_p)// ','// t2c( SIZE(field_in(:,k)))) 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))) 3853 IF (.NOT.this%trans%extrap) THEN 3854 CALL nnseti( 'ext', 0) 3855 CALL nnsetr( 'nul', rmiss) 3858 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, & 3859 this%outnx, this%outny, REAL(this%inter_x(:,1)), & 3860 REAL(this%inter_y(1,:)), field_out(1,1,k), ier) 3863 CALL l4f_category_log(this%category,l4f_error, & 3864 "Error code from NCAR natgrids: "// t2c(ier)) 3870 CALL l4f_category_log(this%category,l4f_info, & 3871 "insufficient data in gridded region to triangulate") 3875 DEALLOCATE(field_in_p, x_in_p, y_in_p) 3878 CALL l4f_category_log(this%category,l4f_error, & 3879 "libsim compiled without NATGRIDD (ngmath ncarg library)") 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 3891 CALL compute(this, & 3892 reshape(field_in, (/ SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, & 3897 END SUBROUTINE grid_transform_v7d_grid_compute 3912 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp) 3913 REAL, INTENT(in) :: z1,z2,z3,z4 3914 DOUBLE PRECISION, INTENT(in):: x1,y1 3915 DOUBLE PRECISION, INTENT(in):: x3,y3 3916 DOUBLE PRECISION, INTENT(in):: xp,yp 3923 p2 = REAL((yp-y1)/(y3-y1)) 3924 p1 = REAL((xp-x1)/(x3-x1)) 3929 zp = (z6-z5)*(p1)+z5 3935 FUNCTION shapiro (z,zp) RESULT(zs) 3938 REAL, INTENT(in) :: z(:) 3939 REAL, INTENT(in) :: zp Functions that return a trimmed CHARACTER representation of the input variable.
Operatore di valore assoluto di un intervallo.
|