/[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.9 by heimbach, Sat Apr 30 16:20:40 2005 UTC revision 1.20 by jmc, Thu May 10 22:21:08 2007 UTC
# Line 1  Line 1 
1    C $Header$
2    C $Name$
3    
4  #include "EXF_OPTIONS.h"  #include "EXF_OPTIONS.h"
5    
6  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
7  C Flux Coupler using                       C  C Flux Coupler using                       C
8  C Bilinear interpolation of forcing fields C  C Bilinear interpolation of forcing fields C
# Line 9  C added Bicubic (bnc 1/2003) Line 13  C added Bicubic (bnc 1/2003)
13  C                                          C  C                                          C
14  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
15    
16          real*8 function lagran(i,x,a,sp)         _RL FUNCTION LAGRAN(i,x,a,sp)
17    
18          INTEGER i,k,sp          INTEGER i
19          _RS x          _RS x
20          real*8 a(4)          _RL a(4)
21          real*8 numer,denom          INTEGER sp
22    
23          numer = 1.D0  C-      local variables:
24          denom = 1.D0          INTEGER k
25            _RL numer,denom
26    
27            numer = 1. _d 0
28            denom = 1. _d 0
29    
30          do k=1,sp          do k=1,sp
31          if ( k .ne. i) then           if ( k .ne. i) then
32            denom = denom*(a(i) - a(k))            denom = denom*(a(i) - a(k))
33            numer = numer*(x    - a(k))            numer = numer*(x    - a(k))
34          endif           endif
35          enddo          enddo
36    
37          lagran = numer/denom          lagran = numer/denom
38    
39          return         return
40          end         end
41    
42    
43         SUBROUTINE exf_interp(         SUBROUTINE exf_interp(
44       I   infile,       I   infile,
45       I   filePrec,       I   filePrec,
46       O   arrayout,       O   arrayout,
47       I   irecord, xG, yG,       I   irecord, xG_in, yG,
48       I   lon_0, lon_inc,       I   lon_0, lon_inc,
49       I   lat_0, lat_inc,       I   lat_0, lat_inc,
50       I   nx_in, ny_in, method, mythid)       I   nx_in, ny_in, method, mythid)
51    
52        implicit none        implicit none
53    
54  C     infile       = name of the input file (direct access binary)  C  infile      (string)  :: name of the binary input file (direct access)
55  C     filePrec     = file precicision (currently not used, assumes real*4)  C  filePrec    (integer) :: number of bits per word in file (32 or 64)
56  C     arrout       = output arrays (different for each processor)  C  arrout      ( _RL )   :: output array
57  C     irecord      = record number in global file  C  irecord     (integer) :: record number to read
58  C     xG,yG        = coordinates for output grid  C     xG,yG              :: coordinates for output grid to interpolate to
59  C     lon_0, lat_0 = lon and lat of sw corner of global input grid  C     lon_0, lat_0       :: lon and lat of sw corner of global input grid
60  C     lon_inc      = scalar x-grid increment  C     lon_inc            :: scalar x-grid increment
61  C     lat_inc      = vector y-grid increments  C     lat_inc            :: vector y-grid increments
62  C     nx_in, ny_in = input x-grid and y-grid size  C  nx_in,ny_in (integer) :: size in x & y direction of input file to read
63  C     method       = 1 for bilinear 2 for bicubic  C     method             :: 1,11,21 for bilinear; 2,12,22 for bicubic
64  C     mythid       = thread id  C                        :: 1,2 for tracer; 11,12 for U; 21,22 for V
65    C  myThid      (integer) :: My Thread Id number
66  C  C
67    
68  #include "SIZE.h"  #include "SIZE.h"
69  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #ifdef ALLOW_USE_MPI  
 # include "EESUPPORT.h"  
 #endif /* ALLOW_USE_MPI */  
70  #include "PARAMS.h"  #include "PARAMS.h"
71    
72  C subroutine variables  C subroutine variables
73        character*(*) infile        character*(*) infile
74        integer       filePrec, irecord, nx_in, ny_in        integer       filePrec, irecord, nx_in, ny_in
75        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
76        _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)
77        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
78        _RL           lon_0, lon_inc        _RL           lon_0, lon_inc
79        _RL           lat_0, lat_inc(ny_in-1)        _RL           lat_0, lat_inc(ny_in-1)
80        integer       method, mythid        integer       method, mythid
81    
82    C functions
83          external lagran
84          _RL      lagran
85    
86  C local variables  C local variables
       real*8   ne_fac,nw_fac,se_fac,sw_fac  
87        integer  e_ind(snx,sny),w_ind(snx,sny)        integer  e_ind(snx,sny),w_ind(snx,sny)
88        integer  n_ind(snx,sny),s_ind(snx,sny)        integer  n_ind(snx,sny),s_ind(snx,sny)
89        real*8   px_ind(4), py_ind(4), ew_val(4)        _RL      px_ind(4), py_ind(4), ew_val(4)
90        external lagran        _RL      arrayin(-1:nx_in+2 ,      -1:ny_in+2)
91        real*8   lagran        _RL      NorthValue
92        real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)        _RL      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)  
93        integer  i, j, k, l, js, bi, bj, sp, interp_unit        integer  i, j, k, l, js, bi, bj, sp, interp_unit
94          _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
95          _RL      ninety
96          PARAMETER ( ninety = 90. )
97          _RS      threeSixtyRS
98          PARAMETER ( threeSixtyRS = 360. )
99    
100    C     put xG in interval [ lon_0 , lon_0+360 [
101          do bj=myByLo(myThid),myByHi(myThid)
102           do bi=myBxLo(myThid),myBxHi(myThid)
103            do j=1-OLy,sNy+OLy
104             do i=1-OLx,sNx+OLx
105              xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
106         &                  + threeSixtyRS*2.
107              xG(i,j,bi,bj) = lon_0+mod(xG(i,j,bi,bj),threeSixtyRS)
108             enddo
109            enddo
110           enddo
111          enddo
112    
113         call exf_interp_read(         call exf_interp_read(
114       I   infile,       I   infile, filePrec,
      I   filePrec,  
115       O   arrayin,       O   arrayin,
116       I   irecord, xG, yG,       I   irecord, nx_in, ny_in, mythid)
      I   lon_0, lon_inc,  
      I   lat_0, lat_inc,  
      I   nx_in, ny_in, method, mythid)  
117    
118        _BEGIN_MASTER( myThid )  C setup input longitude grid
119          do i=-1,nx_in+2
120           x_in(i) = lon_0 + (i-1)*lon_inc
121          enddo
122    
123  C setup input grid  C setup input latitude grid
124         do i=-1,nx_in+2        y_in(0) = lat_0 - lat_inc(1)
125          x_in(i) = lon_0 + (i-1.)*lon_inc        y_in(-1)= lat_0 - 2.*lat_inc(1)
126         enddo        y_in(1) = lat_0
127         y_in(0) = lat_0 - lat_inc(1)        do j=2,ny_in
128         y_in(-1)= lat_0 - 2.*lat_inc(1)         y_in(j) = y_in(j-1) + lat_inc(j-1)
129         y_in(1) = lat_0        enddo
130         do j=2,ny_in        do j=ny_in+1,ny_in+2
131          y_in(j) = y_in(j-1) + lat_inc(j-1)         if (y_in(j-1).eq.ninety) then
132         enddo          y_in(j) = 2 * ninety - y_in(j-2)
133         y_in(ny_in+1) = y_in(ny_in) + lat_inc(ny_in-1)         else
134         y_in(ny_in+2) = y_in(ny_in) + 2.*lat_inc(ny_in-1)          y_in(j) = min( y_in(j-1)+lat_inc(ny_in-1), ninety )
135           endif
136          enddo
137    
138  C enlarge boundary  C enlarge boundary
139         do j=1,ny_in        do j=1,ny_in
140          arrayin(0,j)       = arrayin(nx_in,j)         arrayin(0,j)       = arrayin(nx_in,j)
141          arrayin(-1,j)      = arrayin(nx_in-1,j)         arrayin(-1,j)      = arrayin(nx_in-1,j)
142          arrayin(nx_in+1,j) = arrayin(1,j)         arrayin(nx_in+1,j) = arrayin(1,j)
143          arrayin(nx_in+2,j) = arrayin(2,j)         arrayin(nx_in+2,j) = arrayin(2,j)
144         enddo        enddo
145         do i=-1,nx_in+2        do i=-1,nx_in+2
146          arrayin(i,0)       = arrayin(i,1)         arrayin(i,0)       = arrayin(i,1)
147          arrayin(i,-1)      = arrayin(i,1)         arrayin(i,-1)      = arrayin(i,1)
148          arrayin(i,ny_in+1) = arrayin(i,ny_in)         arrayin(i,ny_in+1) = arrayin(i,ny_in)
149          arrayin(i,ny_in+2) = arrayin(i,ny_in)         arrayin(i,ny_in+2) = arrayin(i,ny_in)
150         enddo        enddo
151    
152        _END_MASTER( myThid )  C     For tracer (method=1,2) set to northernmost zonal-mean value
153    C     at 90N to avoid sharp zonal gradients near the Pole.
154    C     For U (method=11,12) set to zero at 90N to minimize velocity
155    C     gradient at North Pole
156    C     For V (method=11,12) set to northernmost zonal value at 90N,
157    C     as is already done above in order to allow cross-PoleArctic flow
158          do j=ny_in,ny_in+2
159           if (y_in(j).eq.ninety) then
160            if (method.eq.1 .or. method.eq.2) then
161             NorthValue = 0.
162             do i=1,nx_in
163              NorthValue = NorthValue + arrayin(i,j)
164             enddo
165             NorthValue = NorthValue / nx_in
166             do i=-1,nx_in+2
167              arrayin(i,j) = NorthValue
168             enddo
169            elseif (method.eq.11 .or. method.eq.12) then
170             do i=-1,nx_in+2
171              arrayin(i,j) = 0.
172             enddo
173            endif
174           endif
175          enddo
176    
177        do bj = mybylo(mythid), mybyhi(mythid)        do bj = mybylo(mythid), mybyhi(mythid)
178         do bi = mybxlo(mythid), mybxhi(mythid)         do bi = mybxlo(mythid), mybxhi(mythid)
# Line 150  C check validity of input/output coordin Line 200  C check validity of input/output coordin
200          endif          endif
201  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
202    
203  C compute interpolation indices  C compute interpolation indices
204          do i=1,snx          do i=1,snx
205           do j=1,sny           do j=1,sny
206            if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then            if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then
# Line 171  C compute interpolation indices Line 221  C compute interpolation indices
221           enddo           enddo
222          enddo          enddo
223    
224          if (method .eq. 1) then          if (method.eq.1 .or. method.eq.11 .or. method.eq.21) then
225    
226  C bilinear interpolation  C bilinear interpolation
227           sp = 2           sp = 2
# Line 192  C bilinear interpolation Line 242  C bilinear interpolation
242             enddo             enddo
243            enddo            enddo
244           enddo           enddo
245          elseif (method .eq. 2) then          elseif (method .eq. 2 .or. method.eq.12 .or. method.eq.22) then
246    
247  C bicubic interpolation  C bicubic interpolation
248           sp = 4           sp = 4
# Line 210  C bicubic interpolation Line 260  C bicubic interpolation
260       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)
261       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
262       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)
263       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)
264       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)
265       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)
266              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
267       &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)       &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)
268             enddo             enddo
269            enddo            enddo
# Line 224  C bicubic interpolation Line 274  C bicubic interpolation
274         enddo         enddo
275        enddo        enddo
276    
277          RETURN
278        END        END
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22