grid_transform_compute grid_transform_v7d_grid_compute

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, logical, dimension(:), intent(in), optional zlist ) 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 [in] zlist list of levels actually present in input volume, if provided, levels marked as .FALSE. will not be interpolated regardless of their content (optimisation to avoid interpolating missing data) Definizione alla linea 3145 del file grid_transform_class.F90. ii = ii+1 navg = count(field_in(i:ie,j:je,k) /= rmiss) IF (navg > 0) THEN field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), & mask=(field_in(i:ie,j:je,k) /= rmiss))/navg ENDIF ENDDO ENDDO ENDDO ENDIF ELSE IF (this%trans%sub_type == 'stddev' .OR. & this%trans%sub_type == 'stddevnm1') THEN IF (this%trans%sub_type == 'stddev') THEN nm1 = .false. ELSE nm1 = .true. ENDIF navg = this%trans%box_info%npx*this%trans%box_info%npy DO k = 1, innz IF (.NOT.lzlist(k)) cycle jj = 0 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy je = j+this%trans%box_info%npy-1 jj = jj+1 ii = 0 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx ie = i+this%trans%box_info%npx-1 ii = ii+1 field_out(ii,jj,k) = stat_stddev( & reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1) ENDDO ENDDO ENDDO ELSE IF (this%trans%sub_type == 'max') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle jj = 0 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy je = j+this%trans%box_info%npy-1 jj = jj+1 ii = 0 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx ie = i+this%trans%box_info%npx-1 ii = ii+1 navg = count(field_in(i:ie,j:je,k) /= rmiss) IF (navg > 0) THEN field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), & mask=(field_in(i:ie,j:je,k) /= rmiss)) ENDIF ENDDO ENDDO ENDDO ELSE IF (this%trans%sub_type == 'min') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle jj = 0 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy je = j+this%trans%box_info%npy-1 jj = jj+1 ii = 0 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx ie = i+this%trans%box_info%npx-1 ii = ii+1 navg = count(field_in(i:ie,j:je,k) /= rmiss) IF (navg > 0) THEN field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), & mask=(field_in(i:ie,j:je,k) /= rmiss)) ENDIF ENDDO ENDDO ENDDO ELSE IF (this%trans%sub_type == 'percentile') THEN navg = this%trans%box_info%npx*this%trans%box_info%npy DO k = 1, innz IF (.NOT.lzlist(k)) cycle jj = 0 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy je = j+this%trans%box_info%npy-1 jj = jj+1 ii = 0 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx ie = i+this%trans%box_info%npx-1 ii = ii+1 field_out(ii:ii,jj,k) = stat_percentile( & reshape(field_in(i:ie,j:je,k),(/navg/)), & (/real(this%trans%stat_info%percentile)/)) ENDDO ENDDO ENDDO ELSE IF (this%trans%sub_type == 'frequency') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle jj = 0 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy je = j+this%trans%box_info%npy-1 jj = jj+1 ii = 0 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx ie = i+this%trans%box_info%npx-1 ii = ii+1 navg = count(field_in(i:ie,j:je,k) /= rmiss) IF (navg > 0) THEN field_out(ii,jj,k) = sum(interval_info_valid( & this%trans%interval_info, field_in(i:ie,j:je,k)), & mask=(field_in(i:ie,j:je,k) /= rmiss))/navg ENDIF ENDDO ENDDO ENDDO ENDIF ELSE IF (this%trans%trans_type == 'inter') THEN IF (this%trans%sub_type == 'near') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = & field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k) ENDDO ENDDO ENDDO ELSE IF (this%trans%sub_type == 'bilin') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx IF (c_e(this%inter_index_x(i,j))) THEN z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k) z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k) z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k) z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k) IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j)) y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j)) x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1) y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1) xp=this%inter_xp(i,j) yp=this%inter_yp(i,j) field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) ENDIF ENDIF ENDDO ENDDO ENDDO ELSE IF (this%trans%sub_type == 'shapiro_near') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx IF (c_e(this%inter_index_x(i,j))) THEN IF(this%inter_index_x(i,j)-1>0)THEN z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k) ELSE z(1) = rmiss END IF IF(this%inter_index_x(i,j)+1<=this%outnx)THEN z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k) ELSE z(3) = rmiss END IF IF(this%inter_index_y(i,j)+1<=this%outny)THEN z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k) ELSE z(2) = rmiss END IF IF(this%inter_index_y(i,j)-1>0)THEN z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k) ELSE z(4) = rmiss END IF field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)) ENDIF ENDDO ENDDO ENDDO ENDIF ELSE IF (this%trans%trans_type == 'intersearch') THEN likethis = this likethis%trans%trans_type = 'inter' ! fake type and make a recursive call to compute base field CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in, zlist) CALL getenv('LIBSIM_DISABLEOPTSEARCH', env_var) optsearch = len_trim(env_var) == 0 DO k = 1, innz IF (.NOT.lzlist(k)) cycle IF ((.NOT.all(c_e(field_out(:,:,k)))) .AND. (any(c_e(field_in(:,:,k))))) THEN ! must fill some values DO j = 1, this%outny DO i = 1, this%outnx IF (.NOT.c_e(field_out(i,j,k))) THEN dist = huge(dist) nearcount = 0 IF (optsearch) THEN ! optimized, error-prone algorithm ix = this%inter_index_x(i,j) iy = this%inter_index_y(i,j) DO s = 0, max(this%innx, this%inny) farenough = .true. DO m = iy-s, iy+s, max(2*s, 1) ! y loop on upper and lower frames IF (m < 1 .OR. m > this%inny) cycle DO l = max(1, ix-s), min(this%innx, ix+s) ! x loop on upper and lower frames disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2 IF (c_e(field_in(l,m,k))) THEN IF (disttmp < dist) THEN dist = disttmp field_out(i,j,k) = field_in(l,m,k) nearcount = 1 ELSE IF (disttmp == dist) THEN field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k) nearcount = nearcount + 1 ENDIF ENDIF IF (disttmp < dist) farenough = .false. ENDDO ENDDO DO m = max(1, iy-s+1), min(this%inny, iy+s-1) ! y loop on left and right frames (avoid corners) DO l = ix-s, ix+s, 2*s ! x loop on left and right frames (exchange loops?) IF (l < 1 .OR. l > this%innx) cycle disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2 IF (c_e(field_in(l,m,k))) THEN IF (disttmp < dist) THEN dist = disttmp field_out(i,j,k) = field_in(l,m,k) nearcount = 1 ELSE IF (disttmp == dist) THEN field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k) nearcount = nearcount + 1 ENDIF ENDIF IF (disttmp < dist) farenough = .false. ENDDO ENDDO IF (s > 0 .AND. farenough) EXIT ! nearest point found, do not trust the same point, in case of bilin it could be not the nearest ENDDO ELSE ! linear, simple, slow algorithm DO m = 1, this%inny DO l = 1, this%innx IF (c_e(field_in(l,m,k))) THEN disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2 IF (disttmp < dist) THEN dist = disttmp field_out(i,j,k) = field_in(l,m,k) nearcount = 1 ELSE IF (disttmp == dist) THEN field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k) nearcount = nearcount + 1 ENDIF ENDIF ENDDO ENDDO ENDIF ! average points with same minimum distance IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount ENDIF ENDDO ENDDO ENDIF ENDDO ELSE IF (this%trans%trans_type == 'boxinter' & .OR. this%trans%trans_type == 'polyinter' & .OR. this%trans%trans_type == 'maskinter') THEN IF (this%trans%sub_type == 'average') THEN IF (vartype == var_dir360) THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), & 0.0, 360.0, 5.0, & mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j) ENDDO ENDDO ENDDO ELSE ALLOCATE(nval(this%outnx, this%outny)) field_out(:,:,:) = 0.0 DO k = 1, innz IF (.NOT.lzlist(k)) cycle nval(:,:) = 0 DO j = 1, this%inny DO i = 1, this%innx IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = & field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + & field_in(i,j,k) nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = & nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1 ENDIF ENDDO ENDDO WHERE (nval(:,:) /= 0) field_out(:,:,k) = field_out(:,:,k)/nval(:,:) ELSEWHERE field_out(:,:,k) = rmiss END WHERE ENDDO DEALLOCATE(nval) ENDIF ELSE IF (this%trans%sub_type == 'stddev' .OR. & this%trans%sub_type == 'stddevnm1') THEN IF (this%trans%sub_type == 'stddev') THEN nm1 = .false. ELSE nm1 = .true. ENDIF DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx ! da paura field_out(i:i,j,k) = stat_stddev( & reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), & mask=reshape((this%inter_index_x == i .AND. & this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)), nm1=nm1) ENDDO ENDDO ENDDO ELSE IF (this%trans%sub_type == 'max') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%inny DO i = 1, this%innx IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = & max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), & field_in(i,j,k)) ELSE field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = & field_in(i,j,k) ENDIF ENDIF ENDDO ENDDO ENDDO ELSE IF (this%trans%sub_type == 'min') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%inny DO i = 1, this%innx IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = & min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), & field_in(i,j,k)) ELSE field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = & field_in(i,j,k) ENDIF ENDIF ENDDO ENDDO ENDDO ELSE IF (this%trans%sub_type == 'percentile') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx ! da paura field_out(i:i,j,k) = stat_percentile( & reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), & (/real(this%trans%stat_info%percentile)/), & mask=reshape((this%inter_index_x == i .AND. & this%inter_index_y == j), (/SIZE(field_in(:,:,k))/))) ENDDO ENDDO ENDDO ELSE IF (this%trans%sub_type == 'frequency') THEN ALLOCATE(nval(this%outnx, this%outny)) DO k = 1, innz IF (.NOT.lzlist(k)) cycle field_out(:,:,k) = 0.0 nval(:,:) = 0 DO j = 1, this%inny DO i = 1, this%innx IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = & field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + & interval_info_valid(this%trans%interval_info, field_in(i,j,k)) nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = & nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1 ENDIF ENDDO ENDDO WHERE (nval(:,:) /= 0) field_out(:,:,k) = field_out(:,:,k)/nval(:,:) ELSEWHERE field_out(:,:,k) = rmiss END WHERE ENDDO DEALLOCATE(nval) ENDIF ELSE IF (this%trans%trans_type == 'stencilinter') THEN np = SIZE(this%stencil,1)/2 ns = SIZE(this%stencil) IF (this%trans%sub_type == 'average') THEN IF (vartype == var_dir360) THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx IF (c_e(this%inter_index_x(i,j))) THEN i1 = this%inter_index_x(i,j) - np i2 = this%inter_index_x(i,j) + np j1 = this%inter_index_y(i,j) - np j2 = this%inter_index_y(i,j) + np field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), & 0.0, 360.0, 5.0, & mask=this%stencil(:,:)) ! simpler and equivalent form ! mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:)) ENDIF ENDDO ENDDO ENDDO ELSE !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n) DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx IF (c_e(this%inter_index_x(i,j))) THEN i1 = this%inter_index_x(i,j) - np i2 = this%inter_index_x(i,j) + np j1 = this%inter_index_y(i,j) - np j2 = this%inter_index_y(i,j) + np n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:)) IF (n > 0) THEN field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), & mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n ENDIF ENDIF ENDDO ENDDO ENDDO !$OMP END PARALLEL ENDIF ELSE IF (this%trans%sub_type == 'stddev' .OR. & this%trans%sub_type == 'stddevnm1') THEN IF (this%trans%sub_type == 'stddev') THEN nm1 = .false. ELSE nm1 = .true. ENDIF !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2) DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx IF (c_e(this%inter_index_x(i,j))) THEN i1 = this%inter_index_x(i,j) - np i2 = this%inter_index_x(i,j) + np j1 = this%inter_index_y(i,j) - np j2 = this%inter_index_y(i,j) + np ! da paura field_out(i:i,j,k) = stat_stddev( & reshape(field_in(i1:i2,j1:j2,k), (/ns/)), & mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. & this%stencil(:,:), (/ns/)), nm1=nm1) ENDIF ENDDO ENDDO ENDDO !$OMP END PARALLEL ELSE IF (this%trans%sub_type == 'max') THEN !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n) DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx IF (c_e(this%inter_index_x(i,j))) THEN i1 = this%inter_index_x(i,j) - np i2 = this%inter_index_x(i,j) + np j1 = this%inter_index_y(i,j) - np j2 = this%inter_index_y(i,j) + np n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:)) IF (n > 0) THEN field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), & mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:)) ENDIF ENDIF ENDDO ENDDO ENDDO !$OMP END PARALLEL ELSE IF (this%trans%sub_type == 'min') THEN !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n) DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx IF (c_e(this%inter_index_x(i,j))) THEN i1 = this%inter_index_x(i,j) - np i2 = this%inter_index_x(i,j) + np j1 = this%inter_index_y(i,j) - np j2 = this%inter_index_y(i,j) + np n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:)) IF (n > 0) THEN field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), & mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:)) ENDIF ENDIF ENDDO ENDDO ENDDO !$OMP END PARALLEL ELSE IF (this%trans%sub_type == 'percentile') THEN !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2) DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx IF (c_e(this%inter_index_x(i,j))) THEN i1 = this%inter_index_x(i,j) - np i2 = this%inter_index_x(i,j) + np j1 = this%inter_index_y(i,j) - np j2 = this%inter_index_y(i,j) + np ! da paura field_out(i:i,j,k) = stat_percentile( & reshape(field_in(i1:i2,j1:j2,k), (/ns/)), & (/real(this%trans%stat_info%percentile)/), & mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. & this%stencil(:,:), (/ns/))) ENDIF ENDDO ENDDO ENDDO !$OMP END PARALLEL ELSE IF (this%trans%sub_type == 'frequency') THEN !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n) DO k = 1, innz IF (.NOT.lzlist(k)) cycle DO j = 1, this%outny DO i = 1, this%outnx IF (c_e(this%inter_index_x(i,j))) THEN i1 = this%inter_index_x(i,j) - np i2 = this%inter_index_x(i,j) + np j1 = this%inter_index_y(i,j) - np j2 = this%inter_index_y(i,j) + np n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:)) IF (n > 0) THEN field_out(i,j,k) = sum(interval_info_valid( & this%trans%interval_info, field_in(i1:i2,j1:j2,k)), & mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n ENDIF ENDIF ENDDO ENDDO ENDDO !$OMP END PARALLEL ENDIF ELSE IF (this%trans%trans_type == 'maskgen') THEN DO k = 1, innz WHERE(c_e(this%inter_index_x(:,:))) field_out(:,:,k) = real(this%inter_index_x(:,:)) END WHERE ENDDO ELSE IF (this%trans%trans_type == 'metamorphosis') THEN IF (this%trans%sub_type == 'all') THEN field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/)) ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' & .OR. this%trans%sub_type == 'mask' .OR. this%trans%sub_type == 'index') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle ! this is to sparse-points only, so field_out(:,1,k) is acceptable field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:))) ENDDO ELSE IF (this%trans%sub_type == 'maskvalid' .OR. & this%trans%sub_type == 'maskinvalid') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle WHERE (this%point_mask(:,:)) ! this is almost the same as the previous syntax field_out(:,:,k) = field_in(:,:,k) ! check which is more efficient END WHERE ENDDO ELSE IF (this%trans%sub_type == 'lemaskinvalid') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:)) field_out(:,:,k) = field_in(:,:,k) ELSEWHERE field_out(:,:,k) = rmiss END WHERE ENDDO ELSE IF (this%trans%sub_type == 'ltmaskinvalid') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:)) field_out(:,:,k) = field_in(:,:,k) ELSEWHERE field_out(:,:,k) = rmiss END WHERE ENDDO ELSE IF (this%trans%sub_type == 'gemaskinvalid') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:)) field_out(:,:,k) = field_in(:,:,k) ELSEWHERE field_out(:,:,k) = rmiss END WHERE ENDDO ELSE IF (this%trans%sub_type == 'gtmaskinvalid') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:)) field_out(:,:,k) = field_in(:,:,k) ELSEWHERE field_out(:,:,k) = rmiss END WHERE ENDDO ELSE IF (this%trans%sub_type == 'setinvalidto') THEN DO k = 1, innz IF (.NOT.lzlist(k)) cycle ! should i skip levels here? WHERE (c_e(field_in(:,:,k))) field_out(:,:,k) = field_in(:,:,k) ELSE WHERE field_out(:,:,k) = this%val1 END WHERE ENDDO ELSE IF (this%trans%sub_type == 'settoinvalid') THEN IF (c_e(this%val1) .AND. c_e(this%val2)) THEN WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 & .AND. field_in(:,:,:) <= this%val2) field_out(:,:,:) = rmiss ELSE WHERE field_out(:,:,:) = field_in(:,:,:) END WHERE ELSE IF (c_e(this%val1)) THEN WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1) field_out(:,:,:) = rmiss ELSE WHERE field_out(:,:,:) = field_in(:,:,:) END WHERE ELSE IF (c_e(this%val2)) THEN WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2) field_out(:,:,:) = rmiss ELSE WHERE field_out(:,:,:) = field_in(:,:,:) END WHERE ENDIF ENDIF ELSE IF (this%trans%trans_type == 'vertint') THEN IF (this%trans%sub_type == 'linear') THEN alloc_coord_3d_in_act = .false. IF (ASSOCIATED(this%inter_index_z)) THEN DO k = 1, outnz IF (c_e(this%inter_index_z(k))) THEN z1 = real(this%inter_zp(k)) ! weight for k+1 z2 = real(1.0d0 - this%inter_zp(k)) ! weight for k SELECT CASE(vartype) CASE(var_dir360) DO j = 1, inny DO i = 1, innx IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. & c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN field_out(i,j,k) = & interp_var_360(field_in(i,j,this%inter_index_z(k)), & field_in(i,j,this%inter_index_z(k)+1), z1, z2) ENDIF ENDDO ENDDO CASE(var_press) DO j = 1, inny DO i = 1, innx IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. & c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. & field_in(i,j,this%inter_index_z(k)) > 0. .AND. & field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN field_out(i,j,k) = exp( & log(field_in(i,j,this%inter_index_z(k)))*z2 + & log(field_in(i,j,this%inter_index_z(k)+1))*z1) ENDIF ENDDO ENDDO CASE default DO j = 1, inny DO i = 1, innx IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. & c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN field_out(i,j,k) = & field_in(i,j,this%inter_index_z(k))*z2 + & field_in(i,j,this%inter_index_z(k)+1)*z1 ENDIF ENDDO ENDDO END SELECT ENDIF ENDDO ELSE ! use coord_3d_in IF (ALLOCATED(this%coord_3d_in)) THEN coord_3d_in_act => this%coord_3d_in #ifdef DEBUG CALL l4f_category_log(this%category,l4f_debug, & "Using external vertical coordinate for vertint") #endif ELSE IF (PRESENT(coord_3d_in)) THEN #ifdef DEBUG CALL l4f_category_log(this%category,l4f_debug, & "Using internal vertical coordinate for vertint") #endif IF (this%dolog) THEN ALLOCATE(coord_3d_in_act(SIZE(coord_3d_in,1), & SIZE(coord_3d_in,2), SIZE(coord_3d_in,3))) alloc_coord_3d_in_act = .true. WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0) coord_3d_in_act = log(coord_3d_in) ELSEWHERE coord_3d_in_act = rmiss END WHERE ELSE coord_3d_in_act => coord_3d_in ENDIF ELSE CALL l4f_category_log(this%category,l4f_error, & "Missing internal and external vertical coordinate for vertint") CALL raise_error() RETURN ENDIF ENDIF inused = SIZE(coord_3d_in_act,3) IF (inused < 2) RETURN ! to avoid algorithm failure kkcache = 1 DO k = 1, outnz IF (.NOT.c_e(this%vcoord_out(k))) cycle DO j = 1, inny DO i = 1, innx kfound = imiss DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1 kkup = kkcache + kk kkdown = kkcache - kk + 1 IF (kkdown >= 1) THEN ! search down IF (this%vcoord_out(k) >= & min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. & this%vcoord_out(k) < & max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1))) THEN kkcache = kkdown kfoundin = kkcache kfound = kkcache + this%levshift EXIT ! kk ENDIF ENDIF IF (kkup < inused) THEN ! search up IF (this%vcoord_out(k) >= & min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. & this%vcoord_out(k) < & max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1))) THEN kkcache = kkup kfoundin = kkcache kfound = kkcache + this%levshift EXIT ! kk ENDIF ENDIF ENDDO IF (c_e(kfound)) THEN IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound+1))) THEN z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ & (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin))) z2 = 1.0 - z1 SELECT CASE(vartype) CASE(var_dir360) field_out(i,j,k) = & interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2) CASE(var_press) field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + & log(field_in(i,j,kfound+1))*z1) CASE default field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1 END SELECT ENDIF ELSE ENDIF ENDDO ! i ENDDO ! j ENDDO ! k IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act) ENDIF ELSE IF (this%trans%sub_type == 'linearsparse') THEN IF (.NOT.ASSOCIATED(this%vcoord_in) .AND. .NOT.PRESENT(coord_3d_in)) THEN CALL l4f_category_log(this%category,l4f_error, & "linearsparse interpolation, no input vert coord available") RETURN ENDIF ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz)) DO j = 1, inny DO i = 1, innx IF (ASSOCIATED(this%vcoord_in)) THEN mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) & .AND. c_e(this%vcoord_in(:)) ELSE mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) & .AND. c_e(coord_3d_in(i,j,:)) ENDIF IF (vartype == var_press) THEN mask_in(:) = mask_in(:) .AND. & (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0) ENDIF inused = count(mask_in) IF (inused > 1) THEN IF (ASSOCIATED(this%vcoord_in)) THEN coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in) ELSE coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in) ENDIF IF (vartype == var_press) THEN val_in(1:inused) = log(pack( & field_in(i,j,this%levshift+1:this%levshift+this%levused), & mask=mask_in)) ELSE val_in(1:inused) = pack( & field_in(i,j,this%levshift+1:this%levshift+this%levused), & mask=mask_in) ENDIF kkcache = 1 DO k = 1, outnz kfound = imiss DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1 kkup = kkcache + kk kkdown = kkcache - kk + 1 IF (kkdown >= 1) THEN ! search down IF (this%vcoord_out(k) >= & min(coord_in(kkdown), coord_in(kkdown+1)) .AND. & this%vcoord_out(k) < & max(coord_in(kkdown), coord_in(kkdown+1))) THEN kkcache = kkdown kfoundin = kkcache kfound = kkcache EXIT ! kk ENDIF ENDIF IF (kkup < inused) THEN ! search up IF (this%vcoord_out(k) >= & min(coord_in(kkup), coord_in(kkup+1)) .AND. & this%vcoord_out(k) < & max(coord_in(kkup), coord_in(kkup+1))) THEN kkcache = kkup kfoundin = kkcache kfound = kkcache EXIT ! kk ENDIF ENDIF ENDDO IF (c_e(kfound)) THEN z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ & (coord_in(kfound) - coord_in(kfound-1))) z2 = 1.0 - z1 IF (vartype == var_dir360) THEN field_out(i,j,k) = & interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2) ELSE IF (vartype == var_press) THEN field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1) ELSE field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1 ENDIF ENDIF ENDDO ENDIF ENDDO ENDDO DEALLOCATE(coord_in,val_in) ENDIF ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN field_out(:,:,:) = field_in(:,:,:) ENDIF CONTAINS ! internal function for interpolating directions from 0 to 360 degree ! hope it is inlined by the compiler FUNCTION interp_var_360(v1, v2, w1, w2) REAL,INTENT(in) :: v1, v2, w1, w2 REAL :: interp_var_360 REAL :: lv1, lv2 IF (abs(v1 - v2) > 180.) THEN IF (v1 > v2) THEN lv1 = v1 - 360. lv2 = v2 ELSE lv1 = v1 lv2 = v2 - 360. ENDIF interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.) ELSE interp_var_360 = v1*w2 + v2*w1 ENDIF END FUNCTION interp_var_360 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev) REAL,INTENT(in) :: v1(:,:) REAL,INTENT(in) :: l, h, res LOGICAL,INTENT(in),OPTIONAL :: mask(:,:) REAL :: prev REAL :: m INTEGER :: nh, nl !REAL,PARAMETER :: res = 1.0 m = (l + h)/2. IF ((h - l) <= res) THEN prev = m RETURN ENDIF IF (PRESENT(mask)) THEN nl = count(v1 >= l .AND. v1 < m .AND. mask) nh = count(v1 >= m .AND. v1 < h .AND. mask) ELSE nl = count(v1 >= l .AND. v1 < m) nh = count(v1 >= m .AND. v1 < h) ENDIF IF (nh > nl) THEN prev = find_prevailing_direction(v1, m, h, res) ELSE IF (nl > nh) THEN prev = find_prevailing_direction(v1, l, m, res) ELSE IF (nl == 0 .AND. nh == 0) THEN prev = rmiss ELSE prev = m ENDIF END FUNCTION find_prevailing_direction END SUBROUTINE grid_transform_compute !> Compute the output data array from input data array according to !! the defined transformation. The \a grid_transform object \a this !! must have been properly initialised, so that it contains all the !! information needed for computing the transformation. This is the !! sparse points-to-grid and sparse points-to-sparse points version. SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, & coord_3d_in) TYPE(grid_transform),INTENT(in) :: this !< grid_tranform object REAL, INTENT(in) :: field_in(:,:) !< input array REAL, INTENT(out):: field_out(:,:,:) !< output array TYPE(vol7d_var),INTENT(in),OPTIONAL :: var !< physical variable to be interpolated, if provided, some ad-hoc algorithms may be used where possible REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:) !< input vertical coordinate for vertical interpolation, if not provided by other means real,allocatable :: field_in_p(:),x_in_p(:),y_in_p(:) INTEGER :: inn_p, ier, k INTEGER :: innx, inny, innz, outnx, outny, outnz #ifdef DEBUG call l4f_category_log(this%category,L4F_DEBUG,"start v7d_grid_transform_compute") #endif field_out(:,:,:) = rmiss IF (.NOT.this%valid) THEN call l4f_category_log(this%category,l4f_error, & "refusing to perform a non valid transformation") call raise_error() RETURN ENDIF innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2) outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3) ! check size of field_in, field_out IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation IF (innz /= this%innz) THEN CALL l4f_category_log(this%category,l4f_error,"vertical interpolation") CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//& t2c(this%innz)//" /= "//t2c(innz)) CALL raise_error() RETURN ENDIF IF (outnz /= this%outnz) THEN CALL l4f_category_log(this%category,l4f_error,"vertical interpolation") CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//& t2c(this%outnz)//" /= "//t2c(outnz)) CALL raise_error() RETURN ENDIF IF (innx /= outnx .OR. inny /= outny) THEN CALL l4f_category_log(this%category,l4f_error,"vertical interpolation") CALL l4f_category_log(this%category,l4f_error,"inconsistent hor. sizes: "//& t2c(innx)//","//t2c(inny)//" /= "//& t2c(outnx)//","//t2c(outny)) CALL raise_error() RETURN ENDIF ELSE ! horizontal interpolation IF (innx /= this%innx .OR. inny /= this%inny) THEN CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation") CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//& t2c(this%innx)//","//t2c(this%inny)//" /= "//& t2c(innx)//","//t2c(inny)) CALL raise_error() RETURN ENDIF IF (outnx /= this%outnx .OR. outny /= this%outny) THEN CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation") CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//& t2c(this%outnx)//","//t2c(this%outny)//" /= "//& t2c(outnx)//","//t2c(outny)) CALL raise_error() RETURN ENDIF IF (innz /= outnz) THEN CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation") CALL l4f_category_log(this%category,l4f_error,"inconsistent vert. sizes: "//& t2c(innz)//" /= "//t2c(outnz)) CALL raise_error() RETURN ENDIF ENDIF #ifdef DEBUG call l4f_category_log(this%category,l4f_debug, & "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//':'// & trim(this%trans%sub_type)) #endif IF (this%trans%trans_type == 'inter') THEN IF (this%trans%sub_type == 'linear') THEN #ifdef HAVE_LIBNGMATH ! optimization, allocate only once with a safe size ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny)) DO k = 1, innz inn_p = count(c_e(field_in(:,k))) CALL l4f_category_log(this%category,l4f_info, & "Number of sparse data points: "//t2c(inn_p)//','//t2c(SIZE(field_in(:,k)))) IF (inn_p > 2) THEN field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k))) x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k))) y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k))) IF (.NOT.this%trans%extrap) THEN CALL nnseti('ext', 0) ! 0 = no extrapolation CALL nnsetr('nul', rmiss) ENDIF CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, & ! (1:inn_p) omitted this%outnx, this%outny, real(this%inter_x(:,1)), & ! no f90 interface REAL(this%inter_y(1,:)), field_out(1,1,k), ier) IF (ier /= 0) THEN CALL l4f_category_log(this%category,l4f_error, & "Error code from NCAR natgrids: "//t2c(ier)) CALL raise_error() EXIT ENDIF ! exit loop to deallocate ELSE CALL l4f_category_log(this%category,l4f_info, & "insufficient data in gridded region to triangulate") ENDIF ENDDO DEALLOCATE(field_in_p, x_in_p, y_in_p) #else CALL l4f_category_log(this%category,l4f_error, & "libsim compiled without natgridd(ngmath ncarg library)") CALL raise_error() RETURN #endif ENDIF .OR.ELSE IF (this%trans%trans_type == 'boxinter' & .OR. this%trans%trans_type == 'polyinter' & .OR. this%trans%trans_type == 'vertint' & this%trans%trans_type == 'metamorphosis') THEN ! use the grid-to-grid method CALL compute(this, & RESHAPE(field_in, (/SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, & coord_3d_in) ENDIF END SUBROUTINE grid_transform_v7d_grid_compute ! Bilinear interpolation ! Effettua interpolazione bilineare dati i valori nei punti 1,2,3,4 e ! le coordinate dei punti 1 e 3 oltre a quelle del punto p dove viene ! valutato il campo. !_____________________________________________________________ ! disposizione punti ! 4 3 ! ! p ! ! 1 2 ! _____________________________________________________________ ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp) REAL,INTENT(in) :: z1,z2,z3,z4 ! Z values on the four points DOUBLE PRECISION,INTENT(in):: x1,y1 ! coordinate of the lower left point