/[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.4 by dimitri, Sun Dec 7 07:41:45 2003 UTC revision 1.18 by dimitri, Sun Jul 2 13:02:22 2006 UTC
# Line 36  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Line 36  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
36       I   infile,       I   infile,
37       I   filePrec,       I   filePrec,
38       O   arrayout,       O   arrayout,
39       I   irecord, xG, yG,       I   irecord, xG_in, yG,
40       I   lon_0, lon_inc,       I   lon_0, lon_inc,
41       I   lat_0, lat_inc,       I   lat_0, lat_inc,
42       I   nx_in, ny_in, method, mythid)       I   nx_in, ny_in, method, mythid)
# Line 52  C     lon_0, lat_0 = lon and lat of sw c Line 52  C     lon_0, lat_0 = lon and lat of sw c
52  C     lon_inc      = scalar x-grid increment  C     lon_inc      = scalar x-grid increment
53  C     lat_inc      = vector y-grid increments  C     lat_inc      = vector y-grid increments
54  C     nx_in, ny_in = input x-grid and y-grid size  C     nx_in, ny_in = input x-grid and y-grid size
55  C     method       = 1 for bilinear 2 for bicubic  C     method       = 1,11,21 for bilinear; 2,12,22 for bicubic
56    C                    1,2 for tracer; 11,12 for U; 21,22 for V
57  C     mythid       = thread id  C     mythid       = thread id
58  C  C
59    
60  #include "SIZE.h"  #include "SIZE.h"
61  #include "EEPARAMS.h"  #include "EEPARAMS.h"
62  #ifdef ALLOW_USE_MPI  #include "PARAMS.h"
 # include "EESUPPORT.h"  
 # include "PARAMS.h"  
 #endif /* ALLOW_USE_MPI */  
63    
64  C subroutine variables  C subroutine variables
65        character*(*) infile        character*(*) infile
66        integer       filePrec, irecord, nx_in, ny_in        integer       filePrec, irecord, nx_in, ny_in
67        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
68        _RS           xG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS           xG_in      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
69        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
70        _RL           lon_0, lon_inc        _RL           lon_0, lon_inc
71        _RL           lat_0, lat_inc(ny_in-1)        _RL           lat_0, lat_inc(ny_in-1)
72        integer       method, mythid        integer       method, mythid
73    
74  C local variables  C local variables
75        integer  ierr        integer  e_ind(snx,sny),w_ind(snx,sny)
76        real*8   ne_fac,nw_fac,se_fac,sw_fac        integer  n_ind(snx,sny),s_ind(snx,sny)
       integer  e_ind(snx),w_ind(snx)  
       integer  n_ind(sny),s_ind(sny)  
77        real*8   px_ind(4), py_ind(4), ew_val(4)        real*8   px_ind(4), py_ind(4), ew_val(4)
78        external lagran        external lagran
79        real*8   lagran        real*8   lagran
80        real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)        real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)
81        real*8   x_in   (-1:nx_in+2), y_in(-1:ny_in+2)        real*8   x_in   (-1:nx_in+2), y_in(-1:ny_in+2)
82        integer  i, j, k, l, js, bi, bj, interp_unit, sp        real*8   ninety
83        real*4   global(nx_in,ny_in)        PARAMETER ( ninety = 90. )
84          integer  i, j, k, l, js, bi, bj, sp, interp_unit
85        _BEGIN_MASTER( myThid )        _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
86          _RS      threeSixtyRS, NorthValue
87  C check input arguments        PARAMETER ( threeSixtyRS = 360. )
88         if ( .NOT. (filePrec .EQ. 32) )  
89       &     stop 'stop in exf_interp.F: value of filePrec not allowed'  C     put xG in interval [ lon_0 , lon_0+360 [
90          do bj=myByLo(myThid),myByHi(myThid)
91  C read in input data         do bi=myBxLo(myThid),myBxHi(myThid)
92  #ifdef ALLOW_USE_MPI          do j=1-OLy,sNy+OLy
93         if (useSingleCPUIO) then           do i=1-OLx,sNx+OLx
94              xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
95  C master thread of process 0, only, opens a global file       &                  + threeSixtyRS*2.
96          IF( mpiMyId .EQ. 0 ) THEN            xG(i,j,bi,bj) = lon_0+mod(xG(i,j,bi,bj),threeSixtyRS)
          call mdsfindunit( interp_unit, mythid)  
          open(interp_unit,file=infile,status='old',access='direct',  
      &        recl=nx_in*ny_in*4)  
          read(interp_unit,rec=irecord)  
      &        ((global(i,j),i=1,nx_in),j=1,ny_in)  
          close(interp_unit)  
         ENDIF  
   
 C broadcast to all processes  
         call MPI_BCAST(global,nx_in*ny_in,MPI_REAL,  
      &       0,MPI_COMM_MODEL,ierr)  
         do j=1,ny_in  
          do i=1,nx_in  
           arrayin(i,j)=global(i,j)  
97           enddo           enddo
98          enddo          enddo
99           enddo
100          enddo
101    
102           call exf_interp_read(
103         I   infile, filePrec,
104         O   arrayin,
105         I   irecord, nx_in, ny_in, mythid)
106          _BARRIER
107    
108    C     _BEGIN_MASTER( myThid )
109    
110    C setup input longitude grid
111          do i=-1,nx_in+2
112           x_in(i) = lon_0 + (i-1)*lon_inc
113          enddo
114    
115    C setup input latitude grid
116          y_in(0) = lat_0 - lat_inc(1)
117          y_in(-1)= lat_0 - 2.*lat_inc(1)
118          y_in(1) = lat_0
119          do j=2,ny_in
120           y_in(j) = y_in(j-1) + lat_inc(j-1)
121          enddo
122          do j=ny_in+1,ny_in+2
123           if (y_in(j-1).eq.ninety) then
124            y_in(j) = 2 * ninety - y_in(j-2)
125         else         else
126  #endif /* ALLOW_USE_MPI */          y_in(j) = min( y_in(j-1)+lat_inc(ny_in-1), ninety )
   
         call mdsfindunit( interp_unit, mythid)  
         open(interp_unit,file=infile,status='old',access='direct',  
      &       recl=nx_in*ny_in*4)  
         read(interp_unit,rec=irecord)  
      &       ((arrayin(i,j),i=1,nx_in),j=1,ny_in)  
         close(interp_unit)  
   
 #ifdef ALLOW_USE_MPI  
127         endif         endif
128  #endif /* ALLOW_USE_MPI */        enddo
   
 #ifdef _BYTESWAPIO  
        call MDS_BYTESWAPR4((nx_in+4)*(ny_in+4), arrayin )  
 #endif /* _BYTESWAPIO */  
   
 C setup input grid  
        do i=-1,nx_in+2  
         x_in(i) = lon_0 + (i-1.)*lon_inc  
        enddo  
        y_in(0) = lat_0 - lat_inc(1)  
        y_in(-1)= lat_0 - 2.*lat_inc(1)  
        y_in(1) = lat_0  
        do j=2,ny_in  
         y_in(j) = y_in(j-1) + lat_inc(j-1)  
        enddo  
        y_in(ny_in+1) = y_in(ny_in) + lat_inc(ny_in-1)  
        y_in(ny_in+2) = y_in(ny_in) + 2.*lat_inc(ny_in-1)  
129    
130  C enlarge boundary  C enlarge boundary
131         do j=1,ny_in        do j=1,ny_in
132          arrayin(0,j)       = arrayin(nx_in,j)         arrayin(0,j)       = arrayin(nx_in,j)
133          arrayin(-1,j)      = arrayin(nx_in-1,j)         arrayin(-1,j)      = arrayin(nx_in-1,j)
134          arrayin(nx_in+1,j) = arrayin(1,j)         arrayin(nx_in+1,j) = arrayin(1,j)
135          arrayin(nx_in+2,j) = arrayin(2,j)         arrayin(nx_in+2,j) = arrayin(2,j)
136         enddo        enddo
137         do i=-1,nx_in+2        do i=-1,nx_in+2
138          arrayin(i,0)       = arrayin(i,1)         arrayin(i,0)       = arrayin(i,1)
139          arrayin(i,-1)      = arrayin(i,1)         arrayin(i,-1)      = arrayin(i,1)
140          arrayin(i,ny_in+1) = arrayin(i,ny_in)         arrayin(i,ny_in+1) = arrayin(i,ny_in)
141          arrayin(i,ny_in+2) = arrayin(i,ny_in)         arrayin(i,ny_in+2) = arrayin(i,ny_in)
142         enddo        enddo
143    
144        _END_MASTER( myThid )  C     For tracer (method=1,2) set to northernmost zonal-mean value
145    C     at 90N to avoid sharp zonal gradients near the Pole.
146    C     For U (method=11,12) set to zero at 90N to minimize velocity
147    C     gradient at North Pole
148    C     For V (method=11,12) set to northernmost zonal value at 90N,
149    C     as is already done above in order to allow cross-PoleArctic flow
150          do j=ny_in,ny_in+2
151           if (y_in(j).eq.ninety) then
152            if (method.eq.1 .or. method.eq.2) then
153             NorthValue = 0
154             do i=1,nx_in
155              NorthValue = NorthValue + arrayin(i,j)
156             enddo
157             NorthValue = NorthValue / nx_in
158             do i=-1,nx_in+2
159              arrayin(i,j) = NorthValue
160             enddo
161            elseif (method.eq.11 .or. method.eq.12) then
162             do i=-1,nx_in+2
163              arrayin(i,j) = 0
164             enddo
165            endif
166           endif
167          enddo
168    
169    C     _END_MASTER( myThid )
170          
171        do bj = mybylo(mythid), mybyhi(mythid)        do bj = mybylo(mythid), mybyhi(mythid)
172         do bi = mybxlo(mythid), mybxhi(mythid)         do bi = mybxlo(mythid), mybxhi(mythid)
173    
174  C check validity of input/output coordinates  C check validity of input/output coordinates
175          if ( xG(1,1  ,bi,bj) .le. x_in(0)       .or.  #ifdef ALLOW_DEBUG
176       &       xG(snx,1,bi,bj) .ge. x_in(nx_in+1) .or.          if ( debugLevel .ge. debLevB ) then
177       &       yG(1,1  ,bi,bj) .lt. y_in(1)       .or.           do i=1,snx
178       &       yG(1,sny,bi,bj) .gt. y_in(ny_in) ) then            do j=1,sny
179             print*,'ERROR in S/R EXF_INTERP:'             if ( xG(i,j,bi,bj) .lt. x_in(0)         .or.
180             print*,'   input grid must encompass output grid.'       &          xG(i,j,bi,bj) .ge. x_in(nx_in+1)   .or.
181             STOP   ' ABNORMAL END: S/R EXF_INTERP'       &          yG(i,j,bi,bj) .lt. y_in(0)         .or.
182         &          yG(i,j,bi,bj) .ge. y_in(ny_in+1) ) then
183                  print*,'ERROR in S/R EXF_INTERP:'
184                  print*,'   input grid must encompass output grid.'
185                  print*,'i,j,bi,bj'      ,i,j,bi,bj
186                  print*,'xG,yG'          ,xG(i,j,bi,bj),yG(i,j,bi,bj)
187                  print*,'nx_in,ny_in'    ,nx_in        ,ny_in
188                  print*,'x_in(0,nx_in+1)',x_in(0)      ,x_in(nx_in+1)
189                  print*,'y_in(0,ny_in+1)',y_in(0)      ,y_in(ny_in+1)
190                  STOP   ' ABNORMAL END: S/R EXF_INTERP'
191                 endif
192               enddo
193             enddo
194          endif          endif
195    #endif /* ALLOW_DEBUG */
196    
197  C compute interpolation indices  C compute interpolation indices
198          do i=1,snx          do i=1,snx
199           if (xG(i,1,bi,bj)-x_in(1) .ge. 0.) then           do j=1,sny
200              w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc) + 1            if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then
201           else             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1
202              w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc)            else
203           endif             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)
204           e_ind(i) = w_ind(i) + 1            endif
205          enddo            e_ind(i,j) = w_ind(i,j) + 1
206          js = ny_in/2            js = ny_in*.5
207          do j=1,sny            do while (yG(i,j,bi,bj) .lt. y_in(js))
208           do while (yG(1,j,bi,bj) .lt. y_in(js))             js = (js - 1)*.5
209              js = (js + 1)/2            enddo
210           enddo            do while (yG(i,j,bi,bj) .ge. y_in(js+1))
211           do while (yG(1,j,bi,bj) .ge. y_in(js+1))             js = js + 1
212              js = js + 1            enddo
213              s_ind(i,j) = js
214              n_ind(i,j) = js + 1
215           enddo           enddo
          s_ind(j) = js  
          n_ind(j) = js + 1  
216          enddo          enddo
217    
218          if (method .eq. 1) then          if (method.eq.1 .or. method.eq.11 .or. method.eq.21) then
219    
220  C bilinear interpolation  C bilinear interpolation
221           sp = 2           sp = 2
# Line 204  C bilinear interpolation Line 223  C bilinear interpolation
223            do i=1,snx            do i=1,snx
224             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
225             do l=0,1             do l=0,1
226              px_ind(l+1) = x_in(w_ind(i)+l)              px_ind(l+1) = x_in(w_ind(i,j)+l)
227              py_ind(l+1) = y_in(s_ind(j)+l)              py_ind(l+1) = y_in(s_ind(i,j)+l)
228             enddo             enddo
229             do k=1,2             do k=1,2
230              ew_val(k) = arrayin(w_ind(i),s_ind(j)+k-1)              ew_val(k) = arrayin(w_ind(i,j),s_ind(i,j)+k-1)
231       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
232       &             +arrayin(e_ind(i),s_ind(j)+k-1)       &             +arrayin(e_ind(i,j),s_ind(i,j)+k-1)
233       &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
234              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
235       &             +ew_val(k)*lagran(k,yG(1,j,bi,bj),py_ind,sp)       &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)
236             enddo             enddo
237            enddo            enddo
238           enddo           enddo
239          elseif (method .eq. 2) then          elseif (method .eq. 2 .or. method.eq.12 .or. method.eq.22) then
240    
241  C bicubic interpolation  C bicubic interpolation
242           sp = 4           sp = 4
# Line 225  C bicubic interpolation Line 244  C bicubic interpolation
244            do i=1,snx            do i=1,snx
245             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
246             do l=-1,2             do l=-1,2
247              px_ind(l+2) = x_in(w_ind(i)+l)              px_ind(l+2) = x_in(w_ind(i,j)+l)
248              py_ind(l+2) = y_in(s_ind(j)+l)              py_ind(l+2) = y_in(s_ind(i,j)+l)
249             enddo             enddo
250             do k=1,4             do k=1,4
251              ew_val(k) =              ew_val(k) =
252       &             arrayin(w_ind(i)-1,s_ind(j)+k-2)       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)
253       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
254       &             +arrayin(w_ind(i)  ,s_ind(j)+k-2)       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)
255       &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
256       &             +arrayin(e_ind(i)  ,s_ind(j)+k-2)       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)
257       &             *lagran(3,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)
258       &             +arrayin(e_ind(i)+1,s_ind(j)+k-2)       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)
259       &             *lagran(4,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)
260              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
261       &             +ew_val(k)*lagran(k,yG(1,j,bi,bj),py_ind,sp)       &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)
262             enddo             enddo
263            enddo            enddo
264           enddo           enddo
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.22