pkg/exf/exf_interp.F

C $Header:  $
C $Name:  $

#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_in, 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,11,21 for bilinear; 2,12,22 for bicubic
C                    1,2 for tracer; 11,12 for U; 21,22 for V
C     mythid       = thread id
C

#include "SIZE.h"
#include "EEPARAMS.h"
#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_in      (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 local variables
      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)
      real*8   ninety
      PARAMETER ( ninety = 90. )
      integer  i, j, k, l, js, bi, bj, sp, interp_unit
      _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS      threeSixtyRS, NorthValue
      PARAMETER ( threeSixtyRS = 360. )

C     put xG in interval [ lon_0 , lon_0+360 [
      do bj=myByLo(myThid),myByHi(myThid)
       do bi=myBxLo(myThid),myBxHi(myThid)
        do j=1-OLy,sNy+OLy
         do i=1-OLx,sNx+OLx
          xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
     &                  + threeSixtyRS*2.
          xG(i,j,bi,bj) = lon_0+mod(xG(i,j,bi,bj),threeSixtyRS)
         enddo
        enddo
       enddo
      enddo

       call exf_interp_read(
     I   infile, filePrec,
     O   arrayin,
     I   irecord, nx_in, ny_in, mythid)
      _BARRIER

C     _BEGIN_MASTER( myThid )

C setup input longitude grid
      do i=-1,nx_in+2
       x_in(i) = lon_0 + (i-1)*lon_inc
      enddo

C setup input latitude grid
      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
      do j=ny_in+1,ny_in+2
       if (y_in(j-1).eq.ninety) then
        y_in(j) = 2 * ninety - y_in(j-2)
       else
        y_in(j) = min( y_in(j-1)+lat_inc(ny_in-1), ninety )
       endif
      enddo

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

C     For tracer (method=1,2) set to northernmost zonal-mean value
C     at 90N to avoid sharp zonal gradients near the Pole.
C     For U (method=11,12) set to zero at 90N to minimize velocity
C     gradient at North Pole
C     For V (method=11,12) set to northernmost zonal value at 90N,
C     as is already done above in order to allow cross-PoleArctic flow
      do j=ny_in,ny_in+2
       if (y_in(j).eq.ninety) then
        if (method.eq.1 .or. method.eq.2) then
         NorthValue = 0
         do i=1,nx_in
          NorthValue = NorthValue + arrayin(i,j)
         enddo
         NorthValue = NorthValue / nx_in
         do i=-1,nx_in+2
          arrayin(i,j) = NorthValue
         enddo
        elseif (method.eq.11 .or. method.eq.12) then
         do i=-1,nx_in+2
          arrayin(i,j) = 0
         enddo
        endif
       endif
      enddo

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