/[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.19 by jmc, Mon Apr 16 23:27:21 2007 UTC revision 1.20 by jmc, Thu May 10 22:21:08 2007 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  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 12  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    C-      local variables:
24            INTEGER k
25            _RL numer,denom
26    
27          numer = 1.D0          numer = 1. _d 0
28          denom = 1.D0          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(
# Line 46  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Line 51  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
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,11,21 for bilinear; 2,12,22 for bicubic  C     method             :: 1,11,21 for bilinear; 2,12,22 for bicubic
64  C                    1,2 for tracer; 11,12 for U; 21,22 for V  C                        :: 1,2 for tracer; 11,12 for U; 21,22 for V
65  C     mythid       = thread id  C  myThid      (integer) :: My Thread Id number
66  C  C
67    
68  #include "SIZE.h"  #include "SIZE.h"
# Line 68  C subroutine variables Line 73  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_in      (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
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)  
       real*8   ninety  
       PARAMETER ( ninety = 90. )  
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)        _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
95        _RS      threeSixtyRS, NorthValue        _RL      ninety
96          PARAMETER ( ninety = 90. )
97          _RS      threeSixtyRS
98        PARAMETER ( threeSixtyRS = 360. )        PARAMETER ( threeSixtyRS = 360. )
99    
100  C     put xG in interval [ lon_0 , lon_0+360 [  C     put xG in interval [ lon_0 , lon_0+360 [
# Line 106  C     put xG in interval [ lon_0 , lon_0 Line 114  C     put xG in interval [ lon_0 , lon_0
114       I   infile, filePrec,       I   infile, filePrec,
115       O   arrayin,       O   arrayin,
116       I   irecord, nx_in, ny_in, mythid)       I   irecord, nx_in, ny_in, mythid)
       _BARRIER  
   
 C     _BEGIN_MASTER( myThid )  
117    
118  C setup input longitude grid  C setup input longitude grid
119        do i=-1,nx_in+2        do i=-1,nx_in+2
# Line 141  C enlarge boundary Line 146  C enlarge boundary
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  C     For tracer (method=1,2) set to northernmost zonal-mean value  C     For tracer (method=1,2) set to northernmost zonal-mean value
# Line 153  C     as is already done above in order Line 158  C     as is already done above in order
158        do j=ny_in,ny_in+2        do j=ny_in,ny_in+2
159         if (y_in(j).eq.ninety) then         if (y_in(j).eq.ninety) then
160          if (method.eq.1 .or. method.eq.2) then          if (method.eq.1 .or. method.eq.2) then
161           NorthValue = 0           NorthValue = 0.
162           do i=1,nx_in           do i=1,nx_in
163            NorthValue = NorthValue + arrayin(i,j)            NorthValue = NorthValue + arrayin(i,j)
164           enddo           enddo
# Line 163  C     as is already done above in order Line 168  C     as is already done above in order
168           enddo           enddo
169          elseif (method.eq.11 .or. method.eq.12) then          elseif (method.eq.11 .or. method.eq.12) then
170           do i=-1,nx_in+2           do i=-1,nx_in+2
171            arrayin(i,j) = 0            arrayin(i,j) = 0.
172           enddo           enddo
173          endif          endif
174         endif         endif
175        enddo        enddo
176    
 C     _END_MASTER( myThid )  
         
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)
179    
# Line 197  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 257  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 271  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.19  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22