pkg/exf/exf_interp.F

#include "EXF_OPTIONS.h"
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Flux Coupler using                       C
C Bilinear interpolation of forcing fields C
C                                          C
C B. Cheng (12/2002)                       C
C                                          C
C added Bicubic (bnc 1/2003)               C
C                                          C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        real*8 function lagran(i,x,a,sp)

        INTEGER i,k,sp
        _RS x
        real*8 a(4)
        real*8 numer,denom

        numer = 1.D0
        denom = 1.D0

        do k=1,sp
        if ( k .ne. i) then
          denom = denom*(a(i) - a(k))
          numer = numer*(x    - a(k))
        endif
        enddo

        lagran = numer/denom

        return
        end


       SUBROUTINE exf_interp(
     I   infile,
     I   filePrec,
     O   arrayout,
     I   irecord, xG, yG,
     I   lon_0, lon_inc,
     I   lat_0, lat_inc,
     I   nx_in, ny_in, method, mythid)

      implicit none

C     infile       = name of the input file (direct access binary)
C     filePrec     = file precicision (currently not used, assumes real*4)
C     arrout       = output arrays (different for each processor)
C     irecord      = record number in global file
C     xG,yG        = coordinates for output grid
C     lon_0, lat_0 = lon and lat of sw corner of global input grid
C     lon_inc      = scalar x-grid increment
C     lat_inc      = vector y-grid increments
C     nx_in, ny_in = input x-grid and y-grid size
C     method       = 1 for bilinear 2 for bicubic
C     mythid       = thread id
C

#include "SIZE.h"
#include "EEPARAMS.h"
#ifdef ALLOW_USE_MPI
# include "EESUPPORT.h"
#endif /* ALLOW_USE_MPI */
#include "PARAMS.h"

C subroutine variables
      character*(*) infile
      integer       filePrec, irecord, nx_in, ny_in
      _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS           xG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL           lon_0, lon_inc
      _RL           lat_0, lat_inc(ny_in-1)
      integer       method, mythid

C Functions
      integer MDS_RECLEN

C local variables
      integer  ierr, length_of_rec
      real*8   ne_fac,nw_fac,se_fac,sw_fac
      integer  e_ind(snx,sny),w_ind(snx,sny)
      integer  n_ind(snx,sny),s_ind(snx,sny)
      real*8   px_ind(4), py_ind(4), ew_val(4)
      external lagran
      real*8   lagran
      real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)
      real*8   x_in   (-1:nx_in+2), y_in(-1:ny_in+2)
      integer  i, j, k, l, js, bi, bj, sp, interp_unit
      real*4   global(nx_in,ny_in)

      _BEGIN_MASTER( myThid )

C check input arguments
       if ( .NOT. (filePrec .EQ. 32) )
     &     stop 'stop in exf_interp.F: value of filePrec not allowed'

C read in input data
#ifdef ALLOW_USE_MPI
       if (useSingleCPUIO) then

C master thread of process 0, only, opens a global file
        IF( mpiMyId .EQ. 0 ) THEN
         call mdsfindunit( interp_unit, mythid)
         length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, mythid )
         open(interp_unit,file=infile,status='old',access='direct',
     &        recl=length_of_rec)
         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)
         enddo
        enddo

       else
#endif /* ALLOW_USE_MPI */

        call mdsfindunit( interp_unit, mythid)
        length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, mythid )
        open(interp_unit,file=infile,status='old',access='direct',
     &       recl=length_of_rec)
        read(interp_unit,rec=irecord)
     &       ((arrayin(i,j),i=1,nx_in),j=1,ny_in)
        close(interp_unit)

#ifdef ALLOW_USE_MPI
       endif
#endif /* ALLOW_USE_MPI */

#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)

C enlarge boundary
       do j=1,ny_in
        arrayin(0,j)       = arrayin(nx_in,j)
        arrayin(-1,j)      = arrayin(nx_in-1,j)
        arrayin(nx_in+1,j) = arrayin(1,j)
        arrayin(nx_in+2,j) = arrayin(2,j)
       enddo
       do i=-1,nx_in+2
        arrayin(i,0)       = arrayin(i,1)
        arrayin(i,-1)      = arrayin(i,1)
        arrayin(i,ny_in+1) = arrayin(i,ny_in)
        arrayin(i,ny_in+2) = arrayin(i,ny_in) 
       enddo

      _END_MASTER( myThid )

      do bj = mybylo(mythid), mybyhi(mythid)
       do bi = mybxlo(mythid), mybxhi(mythid)

C check validity of input/output coordinates
#ifdef ALLOW_DEBUG
        if ( debugLevel .ge. debLevB ) then
         do i=1,snx
          do j=1,sny
           if ( xG(i,j,bi,bj) .lt. x_in(0)         .or.
     &          xG(i,j,bi,bj) .ge. x_in(nx_in+1)   .or.
     &          yG(i,j,bi,bj) .lt. y_in(0)         .or.
     &          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
        endif
#endif /* ALLOW_DEBUG */

C compute interpolation indices 
        do i=1,snx
         do j=1,sny
          if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then
           w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1
          else
           w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)
          endif
          e_ind(i,j) = w_ind(i,j) + 1
          js = ny_in*.5
          do while (yG(i,j,bi,bj) .lt. y_in(js))
           js = (js - 1)*.5
          enddo
          do while (yG(i,j,bi,bj) .ge. y_in(js+1))
           js = js + 1
          enddo
          s_ind(i,j) = js
          n_ind(i,j) = js + 1
         enddo
        enddo

        if (method .eq. 1) then

C bilinear interpolation
         sp = 2
         do j=1,sny
          do i=1,snx
           arrayout(i,j,bi,bj) = 0.
           do l=0,1
            px_ind(l+1) = x_in(w_ind(i,j)+l)
            py_ind(l+1) = y_in(s_ind(i,j)+l)
           enddo
           do k=1,2
            ew_val(k) = arrayin(w_ind(i,j),s_ind(i,j)+k-1)
     &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j),s_ind(i,j)+k-1)
     &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
            arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
     &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)
           enddo
          enddo
         enddo
        elseif (method .eq. 2) then

C bicubic interpolation
         sp = 4
         do j=1,sny
          do i=1,snx
           arrayout(i,j,bi,bj) = 0.
           do l=-1,2
            px_ind(l+2) = x_in(w_ind(i,j)+l)
            py_ind(l+2) = y_in(s_ind(i,j)+l)
           enddo
           do k=1,4
            ew_val(k) =
     &             arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)
     &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)
     &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)
     &             *lagran(3,xG(i,j,bi,bj),px_ind,sp) 
     &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)
     &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)
            arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj) 
     &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)
           enddo
          enddo
         enddo
        else
         stop 'stop in exf_interp.F: interpolation method not supported'
        endif
       enddo
      enddo

      END
1	edhill	1.3	#include "EXF_OPTIONS.h"
2	dimitri	1.1	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3			C Flux Coupler using C
4			C Bilinear interpolation of forcing fields C
5			C C
6			C B. Cheng (12/2002) C
7			C C
8			C added Bicubic (bnc 1/2003) C
9			C C
10			CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
11
12	dimitri	1.2	real*8 function lagran(i,x,a,sp)
13	dimitri	1.1
14			INTEGER i,k,sp
15			_RS x
16			real*8 a(4)
17			real*8 numer,denom
18
19			numer = 1.D0
20			denom = 1.D0
21
22			do k=1,sp
23			if ( k .ne. i) then
24			denom = denom*(a(i) - a(k))
25			numer = numer*(x - a(k))
26			endif
27			enddo
28
29			lagran = numer/denom
30
31			return
32			end
33
34
35			SUBROUTINE exf_interp(
36			I infile,
37			I filePrec,
38			O arrayout,
39	dimitri	1.2	I irecord, xG, yG,
40			I lon_0, lon_inc,
41			I lat_0, lat_inc,
42			I nx_in, ny_in, method, mythid)
43	dimitri	1.1
44	dimitri	1.4	implicit none
45
46	dimitri	1.2	C infile = name of the input file (direct access binary)
47			C filePrec = file precicision (currently not used, assumes real*4)
48	dimitri	1.1	C arrout = output arrays (different for each processor)
49			C irecord = record number in global file
50			C xG,yG = coordinates for output grid
51			C lon_0, lat_0 = lon and lat of sw corner of global input grid
52			C lon_inc = scalar x-grid increment
53			C lat_inc = vector y-grid increments
54			C nx_in, ny_in = input x-grid and y-grid size
55	dimitri	1.2	C method = 1 for bilinear 2 for bicubic
56			C mythid = thread id
57	dimitri	1.1	C
58
59			#include "SIZE.h"
60			#include "EEPARAMS.h"
61	dimitri	1.4	#ifdef ALLOW_USE_MPI
62			# include "EESUPPORT.h"
63			#endif /* ALLOW_USE_MPI */
64	adcroft	1.7	#include "PARAMS.h"
65	dimitri	1.2
66			C subroutine variables
67			character() infile
68			integer filePrec, irecord, nx_in, ny_in
69			_RL arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
70			_RS xG (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
71			_RS yG (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
72			_RL lon_0, lon_inc
73			_RL lat_0, lat_inc(ny_in-1)
74			integer method, mythid
75	dimitri	1.1
76	dimitri	1.8	C Functions
77			integer MDS_RECLEN
78
79	dimitri	1.1	C local variables
80	dimitri	1.8	integer ierr, length_of_rec
81	dimitri	1.2	real*8 ne_fac,nw_fac,se_fac,sw_fac
82	dimitri	1.5	integer e_ind(snx,sny),w_ind(snx,sny)
83			integer n_ind(snx,sny),s_ind(snx,sny)
84	dimitri	1.2	real*8 px_ind(4), py_ind(4), ew_val(4)
85	dimitri	1.1	external lagran
86	dimitri	1.2	real*8 lagran
87			real*4 arrayin(-1:nx_in+2 , -1:ny_in+2)
88			real*8 x_in (-1:nx_in+2), y_in(-1:ny_in+2)
89	dimitri	1.5	integer i, j, k, l, js, bi, bj, sp, interp_unit
90	dimitri	1.4	real*4 global(nx_in,ny_in)
91
92			_BEGIN_MASTER( myThid )
93	dimitri	1.2
94			C check input arguments
95	dimitri	1.4	if ( .NOT. (filePrec .EQ. 32) )
96			& stop 'stop in exf_interp.F: value of filePrec not allowed'
97	dimitri	1.1
98			C read in input data
99	dimitri	1.4	#ifdef ALLOW_USE_MPI
100			if (useSingleCPUIO) then
101
102			C master thread of process 0, only, opens a global file
103			IF( mpiMyId .EQ. 0 ) THEN
104			call mdsfindunit( interp_unit, mythid)
105	dimitri	1.8	length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, mythid )
106	dimitri	1.4	open(interp_unit,file=infile,status='old',access='direct',
107	dimitri	1.8	& recl=length_of_rec)
108	dimitri	1.4	read(interp_unit,rec=irecord)
109			& ((global(i,j),i=1,nx_in),j=1,ny_in)
110			close(interp_unit)
111			ENDIF
112
113			C broadcast to all processes
114			call MPI_BCAST(global,nx_in*ny_in,MPI_REAL,
115			& 0,MPI_COMM_MODEL,ierr)
116			do j=1,ny_in
117			do i=1,nx_in
118			arrayin(i,j)=global(i,j)
119			enddo
120			enddo
121
122			else
123			#endif /* ALLOW_USE_MPI */
124
125			call mdsfindunit( interp_unit, mythid)
126	dimitri	1.8	length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, mythid )
127	dimitri	1.4	open(interp_unit,file=infile,status='old',access='direct',
128	dimitri	1.8	& recl=length_of_rec)
129	dimitri	1.4	read(interp_unit,rec=irecord)
130			& ((arrayin(i,j),i=1,nx_in),j=1,ny_in)
131			close(interp_unit)
132
133			#ifdef ALLOW_USE_MPI
134			endif
135			#endif /* ALLOW_USE_MPI */
136
137	dimitri	1.1	#ifdef _BYTESWAPIO
138	dimitri	1.4	call MDS_BYTESWAPR4((nx_in+4)*(ny_in+4), arrayin )
139			#endif /* _BYTESWAPIO */
140	dimitri	1.1
141			C setup input grid
142	dimitri	1.4	do i=-1,nx_in+2
143			x_in(i) = lon_0 + (i-1.)*lon_inc
144			enddo
145			y_in(0) = lat_0 - lat_inc(1)
146			y_in(-1)= lat_0 - 2.*lat_inc(1)
147			y_in(1) = lat_0
148			do j=2,ny_in
149			y_in(j) = y_in(j-1) + lat_inc(j-1)
150			enddo
151			y_in(ny_in+1) = y_in(ny_in) + lat_inc(ny_in-1)
152			y_in(ny_in+2) = y_in(ny_in) + 2.*lat_inc(ny_in-1)
153	dimitri	1.1
154			C enlarge boundary
155	dimitri	1.4	do j=1,ny_in
156			arrayin(0,j) = arrayin(nx_in,j)
157			arrayin(-1,j) = arrayin(nx_in-1,j)
158			arrayin(nx_in+1,j) = arrayin(1,j)
159			arrayin(nx_in+2,j) = arrayin(2,j)
160			enddo
161			do i=-1,nx_in+2
162			arrayin(i,0) = arrayin(i,1)
163			arrayin(i,-1) = arrayin(i,1)
164			arrayin(i,ny_in+1) = arrayin(i,ny_in)
165			arrayin(i,ny_in+2) = arrayin(i,ny_in)
166			enddo
167
168			_END_MASTER( myThid )
169	dimitri	1.2
170			do bj = mybylo(mythid), mybyhi(mythid)
171			do bi = mybxlo(mythid), mybxhi(mythid)
172
173			C check validity of input/output coordinates
174	dimitri	1.6	#ifdef ALLOW_DEBUG
175			if ( debugLevel .ge. debLevB ) then
176			do i=1,snx
177			do j=1,sny
178			if ( xG(i,j,bi,bj) .lt. x_in(0) .or.
179			& xG(i,j,bi,bj) .ge. x_in(nx_in+1) .or.
180			& yG(i,j,bi,bj) .lt. y_in(0) .or.
181			& yG(i,j,bi,bj) .ge. y_in(ny_in+1) ) then
182			print*,'ERROR in S/R EXF_INTERP:'
183			print*,' input grid must encompass output grid.'
184			print*,'i,j,bi,bj' ,i,j,bi,bj
185			print*,'xG,yG' ,xG(i,j,bi,bj),yG(i,j,bi,bj)
186			print*,'nx_in,ny_in' ,nx_in ,ny_in
187			print*,'x_in(0,nx_in+1)',x_in(0) ,x_in(nx_in+1)
188			print*,'y_in(0,ny_in+1)',y_in(0) ,y_in(ny_in+1)
189			STOP ' ABNORMAL END: S/R EXF_INTERP'
190			endif
191			enddo
192			enddo
193	dimitri	1.2	endif
194	dimitri	1.6	#endif /* ALLOW_DEBUG */
195	dimitri	1.1
196			C compute interpolation indices
197			do i=1,snx
198	dimitri	1.5	do j=1,sny
199			if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then
200			w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1
201			else
202			w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)
203			endif
204			e_ind(i,j) = w_ind(i,j) + 1
205	dimitri	1.6	js = ny_in*.5
206	dimitri	1.5	do while (yG(i,j,bi,bj) .lt. y_in(js))
207	dimitri	1.6	js = (js - 1)*.5
208	dimitri	1.5	enddo
209			do while (yG(i,j,bi,bj) .ge. y_in(js+1))
210			js = js + 1
211			enddo
212			s_ind(i,j) = js
213			n_ind(i,j) = js + 1
214	dimitri	1.2	enddo
215	dimitri	1.1	enddo
216
217	dimitri	1.2	if (method .eq. 1) then
218	dimitri	1.1
219	dimitri	1.2	C bilinear interpolation
220			sp = 2
221			do j=1,sny
222			do i=1,snx
223	dimitri	1.1	arrayout(i,j,bi,bj) = 0.
224	dimitri	1.2	do l=0,1
225	dimitri	1.5	px_ind(l+1) = x_in(w_ind(i,j)+l)
226			py_ind(l+1) = y_in(s_ind(i,j)+l)
227	dimitri	1.1	enddo
228	dimitri	1.2	do k=1,2
229	dimitri	1.5	ew_val(k) = arrayin(w_ind(i,j),s_ind(i,j)+k-1)
230			& *lagran(1,xG(i,j,bi,bj),px_ind,sp)
231			& +arrayin(e_ind(i,j),s_ind(i,j)+k-1)
232			& *lagran(2,xG(i,j,bi,bj),px_ind,sp)
233	dimitri	1.2	arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
234	dimitri	1.5	& +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)
235	dimitri	1.1	enddo
236			enddo
237			enddo
238	dimitri	1.2	elseif (method .eq. 2) then
239	dimitri	1.1
240	dimitri	1.2	C bicubic interpolation
241			sp = 4
242			do j=1,sny
243			do i=1,snx
244	dimitri	1.1	arrayout(i,j,bi,bj) = 0.
245	dimitri	1.2	do l=-1,2
246	dimitri	1.5	px_ind(l+2) = x_in(w_ind(i,j)+l)
247			py_ind(l+2) = y_in(s_ind(i,j)+l)
248	dimitri	1.1	enddo
249	dimitri	1.2	do k=1,4
250			ew_val(k) =
251	dimitri	1.5	& arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)
252			& *lagran(1,xG(i,j,bi,bj),px_ind,sp)
253			& +arrayin(w_ind(i,j) ,s_ind(i,j)+k-2)
254			& *lagran(2,xG(i,j,bi,bj),px_ind,sp)
255			& +arrayin(e_ind(i,j) ,s_ind(i,j)+k-2)
256			& *lagran(3,xG(i,j,bi,bj),px_ind,sp)
257			& +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)
258			& *lagran(4,xG(i,j,bi,bj),px_ind,sp)
259	dimitri	1.2	arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
260	dimitri	1.5	& +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)
261	dimitri	1.1	enddo
262	dimitri	1.2	enddo
263	dimitri	1.1	enddo
264	dimitri	1.2	else
265			stop 'stop in exf_interp.F: interpolation method not supported'
266			endif
267			enddo
268			enddo
269	dimitri	1.1
270			END