pkg/exf/exf_interp.F

C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_interp.F,v 1.20 2007/05/10 22:21:08 jmc Exp $
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

       _RL FUNCTION LAGRAN(i,x,a,sp)

        INTEGER i
        _RS x
        _RL a(4)
        INTEGER sp

C-      local variables:
        INTEGER k
        _RL numer,denom

        numer = 1. _d 0
        denom = 1. _d 0

        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      (string)  :: name of the binary input file (direct access)
C  filePrec    (integer) :: number of bits per word in file (32 or 64)
C  arrout      ( _RL )   :: output array
C  irecord     (integer) :: record number to read
C     xG,yG              :: coordinates for output grid to interpolate to
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 (integer) :: size in x & y direction of input file to read
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      (integer) :: My Thread Id number
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 functions
      external lagran
      _RL      lagran

C local variables
      integer  e_ind(snx,sny),w_ind(snx,sny)
      integer  n_ind(snx,sny),s_ind(snx,sny)
      _RL      px_ind(4), py_ind(4), ew_val(4)
      _RL      arrayin(-1:nx_in+2 ,      -1:ny_in+2)
      _RL      NorthValue
      _RL      x_in   (-1:nx_in+2), y_in(-1:ny_in+2)
      integer  i, j, k, l, js, bi, bj, sp, interp_unit
      _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL      ninety
      PARAMETER ( ninety = 90. )
      _RS      threeSixtyRS
      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)

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

      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

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