/[MITgcm]/MITgcm/pkg/exf/exf_interp.F
ViewVC logotype

Diff of /MITgcm/pkg/exf/exf_interp.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.23 by mlosch, Thu Jan 24 08:29:51 2008 UTC revision 1.24 by jmc, Wed Apr 7 20:52:30 2010 UTC
# Line 42  C-      local variables: Line 42  C-      local variables:
42         RETURN         RETURN
43         END         END
44    
45    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
46    
47         SUBROUTINE exf_interp(         SUBROUTINE exf_interp(
48       I   infile,       I   infile,
# Line 79  C subroutine variables Line 80  C subroutine variables
80        _RS           xG_in   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS           xG_in   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
81        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
82        _RL           lon_0, lon_inc        _RL           lon_0, lon_inc
83        _RL           lat_0, lat_inc(ny_in-1)  c     _RL           lat_0, lat_inc(ny_in-1)
84          _RL           lat_0, lat_inc(*)
85        integer       method, mythid        integer       method, mythid
86    
87  C functions  C functions
# Line 139  C setup input latitude grid Line 141  C setup input latitude grid
141         if (y_in(j-1).eq.ninety) then         if (y_in(j-1).eq.ninety) then
142          y_in(j) = 2 * ninety - y_in(j-2)          y_in(j) = 2 * ninety - y_in(j-2)
143         else         else
144          y_in(j) = min( y_in(j-1)+lat_inc(ny_in-1), ninety )          i = max(1,ny_in-1)
145            y_in(j) = min( y_in(j-1)+lat_inc(i), ninety )
146         endif         endif
147        enddo        enddo
148    
# Line 188  C     as is already done above in order Line 191  C     as is already done above in order
191  C check validity of input/output coordinates  C check validity of input/output coordinates
192  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
193          if ( debugLevel .ge. debLevB ) then          if ( debugLevel .ge. debLevB ) then
194           do i=1,snx           do j=1,sny
195            do j=1,sny            do i=1,snx
196             if ( xG(i,j,bi,bj) .lt. x_in(0)         .or.             if ( xG(i,j,bi,bj) .lt. x_in(0)         .or.
197       &          xG(i,j,bi,bj) .ge. x_in(nx_in+1)   .or.       &          xG(i,j,bi,bj) .ge. x_in(nx_in+1)   .or.
198       &          yG(i,j,bi,bj) .lt. y_in(0)         .or.       &          yG(i,j,bi,bj) .lt. y_in(0)         .or.
# Line 202  C check validity of input/output coordin Line 205  C check validity of input/output coordin
205                print*,'x_in(0,nx_in+1)',x_in(0)      ,x_in(nx_in+1)                print*,'x_in(0,nx_in+1)',x_in(0)      ,x_in(nx_in+1)
206                print*,'y_in(0,ny_in+1)',y_in(0)      ,y_in(ny_in+1)                print*,'y_in(0,ny_in+1)',y_in(0)      ,y_in(ny_in+1)
207                STOP   ' ABNORMAL END: S/R EXF_INTERP'                STOP   ' ABNORMAL END: S/R EXF_INTERP'
208               endif             endif
209             enddo            enddo
210           enddo           enddo
211          endif          endif
212  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
213    
214  C compute interpolation indices  C compute interpolation indices
215          do i=1,snx          do j=1,sny
216           do j=1,sny           do i=1,snx
217            if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then            if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then
218             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1
219            else            else
# Line 220  C compute interpolation indices Line 223  C compute interpolation indices
223           enddo           enddo
224          enddo          enddo
225  #ifndef TARGET_NEC_SX  #ifndef TARGET_NEC_SX
226  C     use the original and more readable variant of the algorithm that  C     use the original and more readable variant of the algorithm that
227  C     has unvectorizable while-loops for each (i,j)  C     has unvectorizable while-loops for each (i,j)
228          do i=1,snx          do j=1,sny
229           do j=1,sny           do i=1,snx
230            js = ny_in*.5            js = ny_in*.5
231            do while (yG(i,j,bi,bj) .lt. y_in(js))            do while (yG(i,j,bi,bj) .lt. y_in(js))
232             js = (js - 1)*.5             js = (js - 1)*.5
# Line 238  C     has unvectorizable while-loops for Line 241  C     has unvectorizable while-loops for
241  C     this variant vectorizes more efficiently than the original one because  C     this variant vectorizes more efficiently than the original one because
242  C     it moves the while loops out of the i,j loops (loop pushing) but  C     it moves the while loops out of the i,j loops (loop pushing) but
243  C     it is ugly and incomprehensible  C     it is ugly and incomprehensible
244          icnt = 0          icnt = 0
245          do j=1,sny          do j=1,sny
246           do i=1,snx           do i=1,snx
247            s_ind(i,j) = ny_in*.5            s_ind(i,j) = ny_in*.5
248            icnt = icnt+1            icnt = icnt+1
249            inx(icnt,1) = i            inx(icnt,1) = i
250            inx(icnt,2) = j            inx(icnt,2) = j
251           enddo           enddo
252          enddo          enddo
253          do while (icnt .gt. 0)          do while (icnt .gt. 0)
254           ii = 0           ii = 0
255  !CDIR NODEP  !CDIR NODEP
256           do ic=1,icnt           do ic=1,icnt
257            i = inx(ic,1)            i = inx(ic,1)
258            j = inx(ic,2)            j = inx(ic,2)
259            if (yG(i,j,bi,bj) .lt. y_in(s_ind(i,j))) then            if (yG(i,j,bi,bj) .lt. y_in(s_ind(i,j))) then
260             s_ind(i,j) = (s_ind(i,j) - 1)*.5             s_ind(i,j) = (s_ind(i,j) - 1)*.5
261             ii = ii+1             ii = ii+1
# Line 260  C     it is ugly and incomprehensible Line 263  C     it is ugly and incomprehensible
263             inx(ii,2) = j             inx(ii,2) = j
264            endif            endif
265           enddo           enddo
266           icnt = ii           icnt = ii
267          enddo          enddo
268          icnt = 0          icnt = 0
269          do j=1,sny          do j=1,sny
270           do i=1,snx           do i=1,snx
271            icnt = icnt+1            icnt = icnt+1
272            inx(icnt,1) = i            inx(icnt,1) = i
273            inx(icnt,2) = j            inx(icnt,2) = j
274           enddo           enddo
275          enddo          enddo
276          do while (icnt .gt. 0)          do while (icnt .gt. 0)
277           ii = 0           ii = 0
278  !CDIR NODEP  !CDIR NODEP
279           do ic=1,icnt           do ic=1,icnt
280            i = inx(ic,1)            i = inx(ic,1)
281            j = inx(ic,2)            j = inx(ic,2)
282            if (yG(i,j,bi,bj) .ge. y_in(s_ind(i,j)+1)) then            if (yG(i,j,bi,bj) .ge. y_in(s_ind(i,j)+1)) then
283             s_ind(i,j) = s_ind(i,j) + 1             s_ind(i,j) = s_ind(i,j) + 1
284             ii = ii+1             ii = ii+1
# Line 283  C     it is ugly and incomprehensible Line 286  C     it is ugly and incomprehensible
286             inx(ii,2) = j             inx(ii,2) = j
287            endif            endif
288           enddo           enddo
289           icnt = ii           icnt = ii
290          enddo          enddo
291  #endif /* TARGET_NEC_SX defined */  #endif /* TARGET_NEC_SX defined */
292          do i=1,snx          do j=1,sny
293           do j=1,sny           do i=1,snx
294            n_ind(i,j) = s_ind(i,j) + 1            n_ind(i,j) = s_ind(i,j) + 1
295           enddo           enddo
296          enddo          enddo
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22