/[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.1 by dimitri, Thu Aug 7 02:31:29 2003 UTC revision 1.2 by dimitri, Fri Aug 15 01:42:44 2003 UTC
# Line 9  C added Bicubic (bnc 1/2003) Line 9  C added Bicubic (bnc 1/2003)
9  C                                          C  C                                          C
10  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
11    
12          real*8 function lagran(i,x,a)          real*8 function lagran(i,x,a,sp)
13    
14          INTEGER i,k,sp          INTEGER i,k,sp
15          _RS x          _RS x
# Line 19  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Line 19  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
19          numer = 1.D0          numer = 1.D0
20          denom = 1.D0          denom = 1.D0
21    
   
 #ifdef BICUBIC  
         sp = 4  
 #else  
         sp = 2  
 #endif  
22          do k=1,sp          do k=1,sp
23          if ( k .ne. i) then          if ( k .ne. i) then
24            denom = denom*(a(i) - a(k))            denom = denom*(a(i) - a(k))
# Line 42  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, 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,mythid)       I   nx_in, ny_in, method, mythid)
43    
44  C  C
45  C *** infile       = name of the input file (global binary, before interp.)  C     infile       = name of the input file (direct access binary)
46  C     filePrec     = file precicision  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)
48  C     irecord      = record number in global file  C     irecord      = record number in global file
49  C     xG,yG        = coordinates for output grid  C     xG,yG        = coordinates for output grid
# Line 57  C     lon_0, lat_0 = lon and lat of sw c Line 51  C     lon_0, lat_0 = lon and lat of sw c
51  C     lon_inc      = scalar x-grid increment  C     lon_inc      = scalar x-grid increment
52  C     lat_inc      = vector y-grid increments  C     lat_inc      = vector y-grid increments
53  C     nx_in, ny_in = input x-grid and y-grid size  C     nx_in, ny_in = input x-grid and y-grid size
54    C     method       = 1 for bilinear 2 for bicubic
55    C     mythid       = thread id
56  C  C
57    
58  #include "SIZE.h"  #include "SIZE.h"
59  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "EESUPPORT.h"  
 #include "exf_param.h"  
   
 C local variables  
60    
61        integer ierr  C subroutine variables
62        integer nx_in, ny_in        character*(*) infile
63        integer irecord, filePrec        integer       filePrec, irecord, nx_in, ny_in
64        _RL arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nsx,nsy)        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
65          _RS           xG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
66        real*8 ne_fac,nw_fac,se_fac,sw_fac        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
67        integer e_ind(snx),w_ind(snx)        _RL           lon_0, lon_inc
68        integer n_ind(sny),s_ind(sny)        _RL           lat_0, lat_inc(ny_in-1)
69          integer       method, mythid
70    
71        real*8 px_ind(4), py_ind(4)  C local variables
72        real*8 ew_val(4)        integer  ierr
73          real*8   ne_fac,nw_fac,se_fac,sw_fac
74          integer  e_ind(snx),w_ind(snx)
75          integer  n_ind(sny),s_ind(sny)
76          real*8   px_ind(4), py_ind(4), ew_val(4)
77        external lagran        external lagran
78        real*8 lagran        real*8   lagran
79          real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)
80        _RL lon_0,lon_inc        real*8   x_in   (-1:nx_in+2), y_in(-1:ny_in+2)
81        _RL lat_0,lat_inc(ny_in-1)        integer  i, j, k, l, js, bi, bj, sp, interp_unit
82        real*4 arrayin(-1:nx_in+2,-1:ny_in+2)  
83        real*8 x_in(-1:nx_in+2),y_in(-1:ny_in+2)  C check input arguments
84          if ( .NOT. (filePrec .EQ. 32) )
85        _RS xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy)       & stop 'stop in exf_interp.F: value of filePrec not allowed'
       _RS yG(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
   
       character*(*) infile  
       integer i, j, k, l, js  
       integer bi,bj,mythid  
       integer interp_unit  
86    
87  C read in input data  C read in input data
   
88        call mdsfindunit( interp_unit, mythid)        call mdsfindunit( interp_unit, mythid)
89        open(interp_unit,file=infile,status='old',access='direct',        open(interp_unit,file=infile,status='old',access='direct',
90       & recl=nx_in*ny_in*4)       & recl=nx_in*ny_in*4)
# Line 105  C read in input data Line 95  C read in input data
95        close(interp_unit)        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    
 C check validity of input/output parameters  
   
       if ( xG(1,1)   .le. x_in(0)        .or.  
      &     xG(snx,1) .ge. x_in(nx_in+1)  .or.  
      &     yG(1,1)   .lt. y_in(1)        .or.  
      &     yG(1,sny) .gt. y_in(ny_in) )  then  
          print*,'ERROR in S/R EXF_INTERP:'  
          print*,'   input grid must encompass output grid.'  
          STOP   ' ABNORMAL END: S/R EXF_INTERP'  
       endif  
   
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
123    
124  C compute interpolation indices        do bj = mybylo(mythid), mybyhi(mythid)
125           do bi = mybxlo(mythid), mybxhi(mythid)
126    
127    C check validity of input/output coordinates
128            if ( xG(1,1  ,bi,bj) .le. x_in(0)       .or.
129         &       xG(snx,1,bi,bj) .ge. x_in(nx_in+1) .or.
130         &       yG(1,1  ,bi,bj) .lt. y_in(1)       .or.
131         &       yG(1,sny,bi,bj) .gt. y_in(ny_in) ) then
132               print*,'ERROR in S/R EXF_INTERP:'
133               print*,'   input grid must encompass output grid.'
134               STOP   ' ABNORMAL END: S/R EXF_INTERP'
135            endif
136    
137    C compute interpolation indices
138          do i=1,snx          do i=1,snx
139            if (xG(i,1)-x_in(1) .ge. 0.) then           if (xG(i,1,bi,bj)-x_in(1) .ge. 0.) then
140              w_ind(i) = int((xG(i,1)-x_in(1))/lon_inc) + 1                                  w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc) + 1
141            else           else
142              w_ind(i) = int((xG(i,1)-x_in(1))/lon_inc)              w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc)
143            endif           endif
144            e_ind(i) = w_ind(i) + 1               e_ind(i) = w_ind(i) + 1
145          enddo          enddo
   
146          js = ny_in/2          js = ny_in/2
147          do j=1,sny          do j=1,sny
148            do while (yG(1,j) .lt. y_in(js))           do while (yG(1,j,bi,bj) .lt. y_in(js))
149            js = (js + 1)/2              js = (js + 1)/2
150            enddo           enddo
151            do while (yG(1,j) .ge. y_in(js+1))           do while (yG(1,j,bi,bj) .ge. y_in(js+1))
152            js = js + 1              js = js + 1
153            enddo           enddo
154              s_ind(j) = js           s_ind(j) = js
155              n_ind(j) = js + 1           n_ind(j) = js + 1
156          enddo          enddo
157    
158  C interpolate          if (method .eq. 1) then
   
       do bj = mybylo(mythid), mybyhi(mythid)  
         do bi = mybxlo(mythid), mybxhi(mythid)  
   
 #ifdef BICUBIC  
         do j=1,sny  
          do i=1,snx  
159    
160    C bilinear interpolation
161             sp = 2
162             do j=1,sny
163              do i=1,snx
164             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
165               do l=0,1
166             do l=-1,2              px_ind(l+1) = x_in(w_ind(i)+l)
167               px_ind(l+2) = x_in(w_ind(i)+l)              py_ind(l+1) = y_in(s_ind(j)+l)
              py_ind(l+2) = y_in(s_ind(j)+l)  
168             enddo             enddo
169               do k=1,2
170             do k=1,4              ew_val(k) = arrayin(w_ind(i),s_ind(j)+k-1)
171               ew_val(k) =       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)
172       &   arrayin(w_ind(i)-1,s_ind(j)+k-2)*lagran(1,xG(i,1),px_ind)       &             +arrayin(e_ind(i),s_ind(j)+k-1)
173       & + arrayin(w_ind(i)  ,s_ind(j)+k-2)*lagran(2,xG(i,1),px_ind)       &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)
174       & + arrayin(e_ind(i)  ,s_ind(j)+k-2)*lagran(3,xG(i,1),px_ind)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
175       & + arrayin(e_ind(i)+1,s_ind(j)+k-2)*lagran(4,xG(i,1),px_ind)       &             +ew_val(k)*lagran(k,yG(1,j,bi,bj),py_ind,sp)
              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)  
      & + ew_val(k)*lagran(k,yG(1,j),py_ind)  
176             enddo             enddo
   
177            enddo            enddo
178           enddo           enddo
179  #else          elseif (method .eq. 2) then
         do j=1,sny  
          do i=1,snx  
180    
181    C bicubic interpolation
182             sp = 4
183             do j=1,sny
184              do i=1,snx
185             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
186               do l=-1,2
187             do l=0,1              px_ind(l+2) = x_in(w_ind(i)+l)
188               px_ind(l+1) = x_in(w_ind(i)+l)              py_ind(l+2) = y_in(s_ind(j)+l)
              py_ind(l+1) = y_in(s_ind(j)+l)  
189             enddo             enddo
190               do k=1,4
191             do k=1,2              ew_val(k) =
192               ew_val(k) =       &             arrayin(w_ind(i)-1,s_ind(j)+k-2)
193       &   arrayin(w_ind(i),s_ind(j)+k-1)*lagran(1,xG(i,1),px_ind)       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)
194       & + arrayin(e_ind(i),s_ind(j)+k-1)*lagran(2,xG(i,1),px_ind)       &             +arrayin(w_ind(i)  ,s_ind(j)+k-2)
195         &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)
196               arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)       &             +arrayin(e_ind(i)  ,s_ind(j)+k-2)
197       & + ew_val(k)*lagran(k,yG(1,j),py_ind)       &             *lagran(3,xG(i,1,bi,bj),px_ind,sp)
198         &             +arrayin(e_ind(i)+1,s_ind(j)+k-2)
199         &             *lagran(4,xG(i,1,bi,bj),px_ind,sp)
200                arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
201         &             +ew_val(k)*lagran(k,yG(1,j,bi,bj),py_ind,sp)
202             enddo             enddo
203              enddo
          enddo  
         enddo  
 #endif  
204           enddo           enddo
205          enddo          else
206             stop 'stop in exf_interp.F: interpolation method not supported'
207            endif
208           enddo
209          enddo
210    
211        END        END
 C ***  
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22