/[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.10 by heimbach, Mon May 2 15:31:46 2005 UTC revision 1.10.2.1 by dimitri, Thu Aug 25 16:06:10 2005 UTC
# Line 41  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Line 41  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
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)
43    
44        implicit none  C
   
45  C     infile       = name of the input file (direct access binary)  C     infile       = name of the input file (direct access binary)
46  C     filePrec     = file precicision (currently not used, assumes real*4)  C     filePrec     = file precicision (currently not used, assumes real*4)
47  C     arrout       = output arrays (different for each processor)  C     arrout       = output arrays (different for each processor)
# Line 58  C Line 57  C
57    
58  #include "SIZE.h"  #include "SIZE.h"
59  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "PARAMS.h"  
60    
61  C subroutine variables  C subroutine variables
62        character*(*) infile        character*(*) infile
# Line 71  C subroutine variables Line 69  C subroutine variables
69        integer       method, mythid        integer       method, mythid
70    
71  C local variables  C local variables
72          integer  ierr
73        real*8   ne_fac,nw_fac,se_fac,sw_fac        real*8   ne_fac,nw_fac,se_fac,sw_fac
74        integer  e_ind(snx,sny),w_ind(snx,sny)        integer  e_ind(snx),w_ind(snx)
75        integer  n_ind(snx,sny),s_ind(snx,sny)        integer  n_ind(sny),s_ind(sny)
76        real*8   px_ind(4), py_ind(4), ew_val(4)        real*8   px_ind(4), py_ind(4), ew_val(4)
77        external lagran        external lagran
78        real*8   lagran        real*8   lagran
# Line 81  C local variables Line 80  C local variables
80        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)
81        integer  i, j, k, l, js, bi, bj, sp, interp_unit        integer  i, j, k, l, js, bi, bj, sp, interp_unit
82    
83         call exf_interp_read(  C check input arguments
84       I   infile,        if ( .NOT. (filePrec .EQ. 32) )
85       I   filePrec,       & stop 'stop in exf_interp.F: value of filePrec not allowed'
86       O   arrayin,  
87       I   irecord, xG, yG,  C read in input data
88       I   lon_0, lon_inc,        call mdsfindunit( interp_unit, mythid)
89       I   lat_0, lat_inc,        open(interp_unit,file=infile,status='old',access='direct',
90       I   nx_in, ny_in, method, mythid)       & recl=nx_in*ny_in*4)
91          read(interp_unit,rec=irecord) ((arrayin(i,j),i=1,nx_in),j=1,ny_in)
92        _BEGIN_MASTER( myThid )  #ifdef _BYTESWAPIO
93          call MDS_BYTESWAPR4((nx_in+4)*(ny_in+4), arrayin )
94    #endif
95          close(interp_unit)
96    
97  C setup input grid  C setup input grid
98         do i=-1,nx_in+2        do i=-1,nx_in+2
99          x_in(i) = lon_0 + (i-1.)*lon_inc         x_in(i) = lon_0 + (i-1.)*lon_inc
100         enddo        enddo
101         y_in(0) = lat_0 - lat_inc(1)        y_in(0) = lat_0 - lat_inc(1)
102         y_in(-1)= lat_0 - 2.*lat_inc(1)        y_in(-1)= lat_0 - 2.*lat_inc(1)
103         y_in(1) = lat_0        y_in(1) = lat_0
104         do j=2,ny_in        do j=2,ny_in
105          y_in(j) = y_in(j-1) + lat_inc(j-1)        y_in(j) = y_in(j-1) + lat_inc(j-1)
106         enddo        enddo
107         y_in(ny_in+1) = y_in(ny_in) + lat_inc(ny_in-1)        y_in(ny_in+1) = y_in(ny_in) + lat_inc(ny_in-1)
108         y_in(ny_in+2) = y_in(ny_in) + 2.*lat_inc(ny_in-1)        y_in(ny_in+2) = y_in(ny_in) + 2.*lat_inc(ny_in-1)
109    
110  C enlarge boundary  C enlarge boundary
111         do j=1,ny_in        do j=1,ny_in
112          arrayin(0,j)       = arrayin(nx_in,j)           arrayin(0,j)       = arrayin(nx_in,j)
113          arrayin(-1,j)      = arrayin(nx_in-1,j)           arrayin(-1,j)      = arrayin(nx_in-1,j)
114          arrayin(nx_in+1,j) = arrayin(1,j)           arrayin(nx_in+1,j) = arrayin(1,j)
115          arrayin(nx_in+2,j) = arrayin(2,j)           arrayin(nx_in+2,j) = arrayin(2,j)
116         enddo        enddo
117         do i=-1,nx_in+2        do i=-1,nx_in+2
118          arrayin(i,0)       = arrayin(i,1)           arrayin(i,0)       = arrayin(i,1)
119          arrayin(i,-1)      = arrayin(i,1)           arrayin(i,-1)      = arrayin(i,1)
120          arrayin(i,ny_in+1) = arrayin(i,ny_in)           arrayin(i,ny_in+1) = arrayin(i,ny_in)
121          arrayin(i,ny_in+2) = arrayin(i,ny_in)           arrayin(i,ny_in+2) = arrayin(i,ny_in)
122         enddo        enddo
   
       _END_MASTER( myThid )  
123    
124        do bj = mybylo(mythid), mybyhi(mythid)        do bj = mybylo(mythid), mybyhi(mythid)
125         do bi = mybxlo(mythid), mybxhi(mythid)         do bi = mybxlo(mythid), mybxhi(mythid)
126    
127  C check validity of input/output coordinates  C check validity of input/output coordinates
128  #ifdef ALLOW_DEBUG          if ( xG(1,1  ,bi,bj) .le. x_in(0)       .or.
129          if ( debugLevel .ge. debLevB ) then       &       xG(snx,1,bi,bj) .ge. x_in(nx_in+1) .or.
130           do i=1,snx       &       yG(1,1  ,bi,bj) .lt. y_in(1)       .or.
131            do j=1,sny       &       yG(1,sny,bi,bj) .gt. y_in(ny_in) ) then
132             if ( xG(i,j,bi,bj) .lt. x_in(0)         .or.             print*,'ERROR in S/R EXF_INTERP:'
133       &          xG(i,j,bi,bj) .ge. x_in(nx_in+1)   .or.             print*,'   input grid must encompass output grid.'
134       &          yG(i,j,bi,bj) .lt. y_in(0)         .or.             STOP   ' ABNORMAL END: S/R EXF_INTERP'
      &          yG(i,j,bi,bj) .ge. y_in(ny_in+1) ) then  
               print*,'ERROR in S/R EXF_INTERP:'  
               print*,'   input grid must encompass output grid.'  
               print*,'i,j,bi,bj'      ,i,j,bi,bj  
               print*,'xG,yG'          ,xG(i,j,bi,bj),yG(i,j,bi,bj)  
               print*,'nx_in,ny_in'    ,nx_in        ,ny_in  
               print*,'x_in(0,nx_in+1)',x_in(0)      ,x_in(nx_in+1)  
               print*,'y_in(0,ny_in+1)',y_in(0)      ,y_in(ny_in+1)  
               STOP   ' ABNORMAL END: S/R EXF_INTERP'  
              endif  
            enddo  
          enddo  
135          endif          endif
 #endif /* ALLOW_DEBUG */  
136    
137  C compute interpolation indices  C compute interpolation indices
138          do i=1,snx          do i=1,snx
139           do j=1,sny           if (xG(i,1,bi,bj)-x_in(1) .ge. 0.) then
140            if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then              w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc) + 1
141             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1           else
142            else              w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc)
143             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)           endif
144            endif           e_ind(i) = w_ind(i) + 1
145            e_ind(i,j) = w_ind(i,j) + 1          enddo
146            js = ny_in*.5          js = ny_in/2
147            do while (yG(i,j,bi,bj) .lt. y_in(js))          do j=1,sny
148             js = (js - 1)*.5           do while (yG(1,j,bi,bj) .lt. y_in(js))
149            enddo              js = (js + 1)/2
150            do while (yG(i,j,bi,bj) .ge. y_in(js+1))           enddo
151             js = js + 1           do while (yG(1,j,bi,bj) .ge. y_in(js+1))
152            enddo              js = js + 1
           s_ind(i,j) = js  
           n_ind(i,j) = js + 1  
153           enddo           enddo
154             s_ind(j) = js
155             n_ind(j) = js + 1
156          enddo          enddo
157    
158          if (method .eq. 1) then          if (method .eq. 1) then
# Line 176  C bilinear interpolation Line 163  C bilinear interpolation
163            do i=1,snx            do i=1,snx
164             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
165             do l=0,1             do l=0,1
166              px_ind(l+1) = x_in(w_ind(i,j)+l)              px_ind(l+1) = x_in(w_ind(i)+l)
167              py_ind(l+1) = y_in(s_ind(i,j)+l)              py_ind(l+1) = y_in(s_ind(j)+l)
168             enddo             enddo
169             do k=1,2             do k=1,2
170              ew_val(k) = arrayin(w_ind(i,j),s_ind(i,j)+k-1)              ew_val(k) = arrayin(w_ind(i),s_ind(j)+k-1)
171       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)
172       &             +arrayin(e_ind(i,j),s_ind(i,j)+k-1)       &             +arrayin(e_ind(i),s_ind(j)+k-1)
173       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)
174              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
175       &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)       &             +ew_val(k)*lagran(k,yG(1,j,bi,bj),py_ind,sp)
176             enddo             enddo
177            enddo            enddo
178           enddo           enddo
# Line 197  C bicubic interpolation Line 184  C bicubic interpolation
184            do i=1,snx            do i=1,snx
185             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
186             do l=-1,2             do l=-1,2
187              px_ind(l+2) = x_in(w_ind(i,j)+l)              px_ind(l+2) = x_in(w_ind(i)+l)
188              py_ind(l+2) = y_in(s_ind(i,j)+l)              py_ind(l+2) = y_in(s_ind(j)+l)
189             enddo             enddo
190             do k=1,4             do k=1,4
191              ew_val(k) =              ew_val(k) =
192       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)       &             arrayin(w_ind(i)-1,s_ind(j)+k-2)
193       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)
194       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)       &             +arrayin(w_ind(i)  ,s_ind(j)+k-2)
195       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)
196       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)       &             +arrayin(e_ind(i)  ,s_ind(j)+k-2)
197       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)       &             *lagran(3,xG(i,1,bi,bj),px_ind,sp)
198       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)       &             +arrayin(e_ind(i)+1,s_ind(j)+k-2)
199       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)       &             *lagran(4,xG(i,1,bi,bj),px_ind,sp)
200              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
201       &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)       &             +ew_val(k)*lagran(k,yG(1,j,bi,bj),py_ind,sp)
202             enddo             enddo
203            enddo            enddo
204           enddo           enddo
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.10.2.1
  ViewVC Help
Powered by ViewVC 1.1.22