pkg/exf/exf_interp.F

C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_interp.F,v 1.22 2008/01/23 16:41:01 mlosch 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

#ifdef TARGET_NEC_SX
!CDIR UNROLL=8
#endif /* TARGET_NEC_SX */
        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
#ifdef TARGET_NEC_SX
      integer  ic, ii, icnt
      integer  inx(snx*sny,2)
      _RL      ew_val1, ew_val2, ew_val3, ew_val4
#endif
      _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
         enddo
        enddo
#ifndef TARGET_NEC_SX
C     use the original and more readable variant of the algorithm that 
C     has unvectorizable while-loops for each (i,j)
        do i=1,snx
         do j=1,sny
          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
         enddo
        enddo
#else /* TARGET_NEC_SX defined */
C     this variant vectorizes more efficiently than the original one because
C     it moves the while loops out of the i,j loops (loop pushing) but
C     it is ugly and incomprehensible
        icnt = 0
        do j=1,sny
         do i=1,snx
          s_ind(i,j) = ny_in*.5
          icnt = icnt+1
          inx(icnt,1) = i
          inx(icnt,2) = j
         enddo
        enddo
        do while (icnt .gt. 0)
         ii = 0
!CDIR NODEP
         do ic=1,icnt
          i = inx(ic,1)
          j = inx(ic,2)
          if (yG(i,j,bi,bj) .lt. y_in(s_ind(i,j))) then
           s_ind(i,j) = (s_ind(i,j) - 1)*.5
           ii = ii+1
           inx(ii,1) = i
           inx(ii,2) = j
          endif
         enddo
         icnt = ii
        enddo
        icnt = 0
        do j=1,sny
         do i=1,snx
          icnt = icnt+1
          inx(icnt,1) = i
          inx(icnt,2) = j
         enddo
        enddo
        do while (icnt .gt. 0)
         ii = 0
!CDIR NODEP
         do ic=1,icnt
          i = inx(ic,1)
          j = inx(ic,2)
          if (yG(i,j,bi,bj) .ge. y_in(s_ind(i,j)+1)) then
           s_ind(i,j) = s_ind(i,j) + 1
           ii = ii+1
           inx(ii,1) = i
           inx(ii,2) = j
          endif
         enddo
         icnt = ii
        enddo
#endif /* TARGET_NEC_SX defined */
        do i=1,snx
         do j=1,sny
          n_ind(i,j) = s_ind(i,j) + 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
#ifndef TARGET_NEC_SX
           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
#else
           ew_val1 = arrayin(w_ind(i,j),s_ind(i,j)+1-1)
     &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j),s_ind(i,j)+1-1)
     &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
           ew_val2 = arrayin(w_ind(i,j),s_ind(i,j)+2-1)
     &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j),s_ind(i,j)+2-1)
     &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
           arrayout(i,j,bi,bj)=
     &             +ew_val1*lagran(1,yG(i,j,bi,bj),py_ind,sp)
     &             +ew_val2*lagran(2,yG(i,j,bi,bj),py_ind,sp)
#endif /* TARGET_NEC_SX defined */
          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
#ifndef TARGET_NEC_SX
           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
#else
           ew_val1 =
     &             arrayin(w_ind(i,j)-1,s_ind(i,j)+1-2)
     &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+1-2)
     &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+1-2)
     &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+1-2)
     &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)
            ew_val2 =
     &             arrayin(w_ind(i,j)-1,s_ind(i,j)+2-2)
     &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+2-2)
     &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+2-2)
     &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+2-2)
     &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)
            ew_val3 =
     &             arrayin(w_ind(i,j)-1,s_ind(i,j)+3-2)
     &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+3-2)
     &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+3-2)
     &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+3-2)
     &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)
            ew_val4 =
     &             arrayin(w_ind(i,j)-1,s_ind(i,j)+4-2)
     &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+4-2)
     &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+4-2)
     &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)
     &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+4-2)
     &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)
            arrayout(i,j,bi,bj)=
     &             +ew_val1*lagran(1,yG(i,j,bi,bj),py_ind,sp)
     &             +ew_val2*lagran(2,yG(i,j,bi,bj),py_ind,sp)
     &             +ew_val3*lagran(3,yG(i,j,bi,bj),py_ind,sp)
     &             +ew_val4*lagran(4,yG(i,j,bi,bj),py_ind,sp)
#endif /* TARGET_NEC_SX defined */
          enddo
         enddo
        else
         stop 'stop in exf_interp.F: interpolation method not supported'
        endif
       enddo
      enddo

      RETURN
      END
1	mlosch	1.23	C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_interp.F,v 1.22 2008/01/23 16:41:01 mlosch 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	mlosch	1.22	#ifdef TARGET_NEC_SX
31			!CDIR UNROLL=8
32			#endif /* TARGET_NEC_SX */
33	dimitri	1.1	do k=1,sp
34	jmc	1.20	if ( k .ne. i) then
35	dimitri	1.1	denom = denom*(a(i) - a(k))
36			numer = numer*(x - a(k))
37	jmc	1.20	endif
38	dimitri	1.1	enddo
39
40			lagran = numer/denom
41
42	jmc	1.21	RETURN
43			END
44	dimitri	1.1
45
46			SUBROUTINE exf_interp(
47			I infile,
48			I filePrec,
49			O arrayout,
50	heimbach	1.13	I irecord, xG_in, yG,
51	dimitri	1.2	I lon_0, lon_inc,
52			I lat_0, lat_inc,
53			I nx_in, ny_in, method, mythid)
54	dimitri	1.1
55	dimitri	1.4	implicit none
56
57	jmc	1.20	C infile (string) :: name of the binary input file (direct access)
58			C filePrec (integer) :: number of bits per word in file (32 or 64)
59			C arrout ( _RL ) :: output array
60			C irecord (integer) :: record number to read
61			C xG,yG :: coordinates for output grid to interpolate to
62			C lon_0, lat_0 :: lon and lat of sw corner of global input grid
63			C lon_inc :: scalar x-grid increment
64			C lat_inc :: vector y-grid increments
65			C nx_in,ny_in (integer) :: size in x & y direction of input file to read
66			C method :: 1,11,21 for bilinear; 2,12,22 for bicubic
67			C :: 1,2 for tracer; 11,12 for U; 21,22 for V
68			C myThid (integer) :: My Thread Id number
69	dimitri	1.1	C
70
71			#include "SIZE.h"
72			#include "EEPARAMS.h"
73	adcroft	1.7	#include "PARAMS.h"
74	dimitri	1.2
75			C subroutine variables
76			character() infile
77			integer filePrec, irecord, nx_in, ny_in
78			_RL arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
79	jmc	1.20	_RS xG_in (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
80	dimitri	1.2	_RS yG (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
81			_RL lon_0, lon_inc
82			_RL lat_0, lat_inc(ny_in-1)
83			integer method, mythid
84	dimitri	1.1
85	jmc	1.20	C functions
86			external lagran
87			_RL lagran
88
89	dimitri	1.1	C local variables
90	dimitri	1.5	integer e_ind(snx,sny),w_ind(snx,sny)
91			integer n_ind(snx,sny),s_ind(snx,sny)
92	jmc	1.20	_RL px_ind(4), py_ind(4), ew_val(4)
93			_RL arrayin(-1:nx_in+2 , -1:ny_in+2)
94			_RL NorthValue
95			_RL x_in (-1:nx_in+2), y_in(-1:ny_in+2)
96	dimitri	1.5	integer i, j, k, l, js, bi, bj, sp, interp_unit
97	mlosch	1.22	#ifdef TARGET_NEC_SX
98			integer ic, ii, icnt
99			integer inx(snx*sny,2)
100	mlosch	1.23	_RL ew_val1, ew_val2, ew_val3, ew_val4
101	mlosch	1.22	#endif
102	heimbach	1.13	_RS xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
103	jmc	1.20	_RL ninety
104			PARAMETER ( ninety = 90. )
105			_RS threeSixtyRS
106	heimbach	1.13	PARAMETER ( threeSixtyRS = 360. )
107	heimbach	1.12
108	jmc	1.14	C put xG in interval [ lon_0 , lon_0+360 [
109	heimbach	1.12	do bj=myByLo(myThid),myByHi(myThid)
110			do bi=myBxLo(myThid),myBxHi(myThid)
111			do j=1-OLy,sNy+OLy
112			do i=1-OLx,sNx+OLx
113	jmc	1.14	xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
114			& + threeSixtyRS*2.
115			xG(i,j,bi,bj) = lon_0+mod(xG(i,j,bi,bj),threeSixtyRS)
116	heimbach	1.12	enddo
117			enddo
118			enddo
119			enddo
120	heimbach	1.9
121			call exf_interp_read(
122	dimitri	1.15	I infile, filePrec,
123	heimbach	1.9	O arrayin,
124	dimitri	1.15	I irecord, nx_in, ny_in, mythid)
125	dimitri	1.2
126	dimitri	1.18	C setup input longitude grid
127			do i=-1,nx_in+2
128			x_in(i) = lon_0 + (i-1)*lon_inc
129			enddo
130	heimbach	1.12
131	dimitri	1.18	C setup input latitude grid
132			y_in(0) = lat_0 - lat_inc(1)
133			y_in(-1)= lat_0 - 2.*lat_inc(1)
134			y_in(1) = lat_0
135			do j=2,ny_in
136			y_in(j) = y_in(j-1) + lat_inc(j-1)
137			enddo
138			do j=ny_in+1,ny_in+2
139			if (y_in(j-1).eq.ninety) then
140			y_in(j) = 2 * ninety - y_in(j-2)
141			else
142			y_in(j) = min( y_in(j-1)+lat_inc(ny_in-1), ninety )
143			endif
144			enddo
145	dimitri	1.1
146			C enlarge boundary
147	dimitri	1.18	do j=1,ny_in
148			arrayin(0,j) = arrayin(nx_in,j)
149			arrayin(-1,j) = arrayin(nx_in-1,j)
150			arrayin(nx_in+1,j) = arrayin(1,j)
151			arrayin(nx_in+2,j) = arrayin(2,j)
152			enddo
153			do i=-1,nx_in+2
154			arrayin(i,0) = arrayin(i,1)
155			arrayin(i,-1) = arrayin(i,1)
156			arrayin(i,ny_in+1) = arrayin(i,ny_in)
157	jmc	1.20	arrayin(i,ny_in+2) = arrayin(i,ny_in)
158	dimitri	1.18	enddo
159	dimitri	1.4
160	dimitri	1.15	C For tracer (method=1,2) set to northernmost zonal-mean value
161			C at 90N to avoid sharp zonal gradients near the Pole.
162			C For U (method=11,12) set to zero at 90N to minimize velocity
163			C gradient at North Pole
164			C For V (method=11,12) set to northernmost zonal value at 90N,
165			C as is already done above in order to allow cross-PoleArctic flow
166	dimitri	1.18	do j=ny_in,ny_in+2
167			if (y_in(j).eq.ninety) then
168	dimitri	1.15	if (method.eq.1 .or. method.eq.2) then
169	jmc	1.20	NorthValue = 0.
170	dimitri	1.15	do i=1,nx_in
171	dimitri	1.18	NorthValue = NorthValue + arrayin(i,j)
172	dimitri	1.15	enddo
173			NorthValue = NorthValue / nx_in
174			do i=-1,nx_in+2
175	dimitri	1.18	arrayin(i,j) = NorthValue
176	dimitri	1.15	enddo
177			elseif (method.eq.11 .or. method.eq.12) then
178			do i=-1,nx_in+2
179	jmc	1.20	arrayin(i,j) = 0.
180	dimitri	1.15	enddo
181			endif
182			endif
183	dimitri	1.18	enddo
184	dimitri	1.15
185	dimitri	1.2	do bj = mybylo(mythid), mybyhi(mythid)
186			do bi = mybxlo(mythid), mybxhi(mythid)
187
188			C check validity of input/output coordinates
189	dimitri	1.6	#ifdef ALLOW_DEBUG
190			if ( debugLevel .ge. debLevB ) then
191			do i=1,snx
192			do j=1,sny
193			if ( xG(i,j,bi,bj) .lt. x_in(0) .or.
194			& xG(i,j,bi,bj) .ge. x_in(nx_in+1) .or.
195			& yG(i,j,bi,bj) .lt. y_in(0) .or.
196			& yG(i,j,bi,bj) .ge. y_in(ny_in+1) ) then
197			print*,'ERROR in S/R EXF_INTERP:'
198			print*,' input grid must encompass output grid.'
199			print*,'i,j,bi,bj' ,i,j,bi,bj
200			print*,'xG,yG' ,xG(i,j,bi,bj),yG(i,j,bi,bj)
201			print*,'nx_in,ny_in' ,nx_in ,ny_in
202			print*,'x_in(0,nx_in+1)',x_in(0) ,x_in(nx_in+1)
203			print*,'y_in(0,ny_in+1)',y_in(0) ,y_in(ny_in+1)
204			STOP ' ABNORMAL END: S/R EXF_INTERP'
205			endif
206			enddo
207			enddo
208	dimitri	1.2	endif
209	dimitri	1.6	#endif /* ALLOW_DEBUG */
210	dimitri	1.1
211	jmc	1.20	C compute interpolation indices
212	dimitri	1.1	do i=1,snx
213	dimitri	1.5	do j=1,sny
214			if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then
215			w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1
216			else
217			w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)
218			endif
219			e_ind(i,j) = w_ind(i,j) + 1
220	mlosch	1.22	enddo
221			enddo
222			#ifndef TARGET_NEC_SX
223			C use the original and more readable variant of the algorithm that
224			C has unvectorizable while-loops for each (i,j)
225			do i=1,snx
226			do j=1,sny
227	dimitri	1.6	js = ny_in*.5
228	dimitri	1.5	do while (yG(i,j,bi,bj) .lt. y_in(js))
229	dimitri	1.6	js = (js - 1)*.5
230	dimitri	1.5	enddo
231			do while (yG(i,j,bi,bj) .ge. y_in(js+1))
232			js = js + 1
233			enddo
234			s_ind(i,j) = js
235	mlosch	1.22	enddo
236			enddo
237			#else /* TARGET_NEC_SX defined */
238			C this variant vectorizes more efficiently than the original one because
239			C it moves the while loops out of the i,j loops (loop pushing) but
240			C it is ugly and incomprehensible
241			icnt = 0
242			do j=1,sny
243			do i=1,snx
244			s_ind(i,j) = ny_in*.5
245			icnt = icnt+1
246			inx(icnt,1) = i
247			inx(icnt,2) = j
248			enddo
249			enddo
250			do while (icnt .gt. 0)
251			ii = 0
252			!CDIR NODEP
253			do ic=1,icnt
254			i = inx(ic,1)
255			j = inx(ic,2)
256			if (yG(i,j,bi,bj) .lt. y_in(s_ind(i,j))) then
257			s_ind(i,j) = (s_ind(i,j) - 1)*.5
258			ii = ii+1
259			inx(ii,1) = i
260			inx(ii,2) = j
261			endif
262			enddo
263			icnt = ii
264			enddo
265			icnt = 0
266			do j=1,sny
267			do i=1,snx
268			icnt = icnt+1
269			inx(icnt,1) = i
270			inx(icnt,2) = j
271			enddo
272			enddo
273			do while (icnt .gt. 0)
274			ii = 0
275			!CDIR NODEP
276			do ic=1,icnt
277			i = inx(ic,1)
278			j = inx(ic,2)
279			if (yG(i,j,bi,bj) .ge. y_in(s_ind(i,j)+1)) then
280			s_ind(i,j) = s_ind(i,j) + 1
281			ii = ii+1
282			inx(ii,1) = i
283			inx(ii,2) = j
284			endif
285			enddo
286			icnt = ii
287			enddo
288			#endif /* TARGET_NEC_SX defined */
289			do i=1,snx
290			do j=1,sny
291			n_ind(i,j) = s_ind(i,j) + 1
292	dimitri	1.2	enddo
293	dimitri	1.1	enddo
294
295	dimitri	1.15	if (method.eq.1 .or. method.eq.11 .or. method.eq.21) then
296	dimitri	1.1
297	dimitri	1.2	C bilinear interpolation
298			sp = 2
299			do j=1,sny
300			do i=1,snx
301	dimitri	1.1	arrayout(i,j,bi,bj) = 0.
302	dimitri	1.2	do l=0,1
303	dimitri	1.5	px_ind(l+1) = x_in(w_ind(i,j)+l)
304			py_ind(l+1) = y_in(s_ind(i,j)+l)
305	dimitri	1.1	enddo
306	mlosch	1.23	#ifndef TARGET_NEC_SX
307	dimitri	1.2	do k=1,2
308	dimitri	1.5	ew_val(k) = arrayin(w_ind(i,j),s_ind(i,j)+k-1)
309			& *lagran(1,xG(i,j,bi,bj),px_ind,sp)
310			& +arrayin(e_ind(i,j),s_ind(i,j)+k-1)
311			& *lagran(2,xG(i,j,bi,bj),px_ind,sp)
312	dimitri	1.2	arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
313	dimitri	1.5	& +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)
314	dimitri	1.1	enddo
315	mlosch	1.23	#else
316			ew_val1 = arrayin(w_ind(i,j),s_ind(i,j)+1-1)
317			& *lagran(1,xG(i,j,bi,bj),px_ind,sp)
318			& +arrayin(e_ind(i,j),s_ind(i,j)+1-1)
319			& *lagran(2,xG(i,j,bi,bj),px_ind,sp)
320			ew_val2 = arrayin(w_ind(i,j),s_ind(i,j)+2-1)
321			& *lagran(1,xG(i,j,bi,bj),px_ind,sp)
322			& +arrayin(e_ind(i,j),s_ind(i,j)+2-1)
323			& *lagran(2,xG(i,j,bi,bj),px_ind,sp)
324			arrayout(i,j,bi,bj)=
325			& +ew_val1*lagran(1,yG(i,j,bi,bj),py_ind,sp)
326			& +ew_val2*lagran(2,yG(i,j,bi,bj),py_ind,sp)
327			#endif /* TARGET_NEC_SX defined */
328	dimitri	1.1	enddo
329			enddo
330	dimitri	1.15	elseif (method .eq. 2 .or. method.eq.12 .or. method.eq.22) then
331	dimitri	1.1
332	dimitri	1.2	C bicubic interpolation
333			sp = 4
334			do j=1,sny
335			do i=1,snx
336	dimitri	1.1	arrayout(i,j,bi,bj) = 0.
337	dimitri	1.2	do l=-1,2
338	dimitri	1.5	px_ind(l+2) = x_in(w_ind(i,j)+l)
339			py_ind(l+2) = y_in(s_ind(i,j)+l)
340	dimitri	1.1	enddo
341	mlosch	1.23	#ifndef TARGET_NEC_SX
342	dimitri	1.2	do k=1,4
343			ew_val(k) =
344	dimitri	1.5	& arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)
345			& *lagran(1,xG(i,j,bi,bj),px_ind,sp)
346			& +arrayin(w_ind(i,j) ,s_ind(i,j)+k-2)
347			& *lagran(2,xG(i,j,bi,bj),px_ind,sp)
348			& +arrayin(e_ind(i,j) ,s_ind(i,j)+k-2)
349	jmc	1.20	& *lagran(3,xG(i,j,bi,bj),px_ind,sp)
350	dimitri	1.5	& +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)
351			& *lagran(4,xG(i,j,bi,bj),px_ind,sp)
352	jmc	1.20	arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
353	dimitri	1.5	& +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)
354	dimitri	1.1	enddo
355	mlosch	1.23	#else
356			ew_val1 =
357			& arrayin(w_ind(i,j)-1,s_ind(i,j)+1-2)
358			& *lagran(1,xG(i,j,bi,bj),px_ind,sp)
359			& +arrayin(w_ind(i,j) ,s_ind(i,j)+1-2)
360			& *lagran(2,xG(i,j,bi,bj),px_ind,sp)
361			& +arrayin(e_ind(i,j) ,s_ind(i,j)+1-2)
362			& *lagran(3,xG(i,j,bi,bj),px_ind,sp)
363			& +arrayin(e_ind(i,j)+1,s_ind(i,j)+1-2)
364			& *lagran(4,xG(i,j,bi,bj),px_ind,sp)
365			ew_val2 =
366			& arrayin(w_ind(i,j)-1,s_ind(i,j)+2-2)
367			& *lagran(1,xG(i,j,bi,bj),px_ind,sp)
368			& +arrayin(w_ind(i,j) ,s_ind(i,j)+2-2)
369			& *lagran(2,xG(i,j,bi,bj),px_ind,sp)
370			& +arrayin(e_ind(i,j) ,s_ind(i,j)+2-2)
371			& *lagran(3,xG(i,j,bi,bj),px_ind,sp)
372			& +arrayin(e_ind(i,j)+1,s_ind(i,j)+2-2)
373			& *lagran(4,xG(i,j,bi,bj),px_ind,sp)
374			ew_val3 =
375			& arrayin(w_ind(i,j)-1,s_ind(i,j)+3-2)
376			& *lagran(1,xG(i,j,bi,bj),px_ind,sp)
377			& +arrayin(w_ind(i,j) ,s_ind(i,j)+3-2)
378			& *lagran(2,xG(i,j,bi,bj),px_ind,sp)
379			& +arrayin(e_ind(i,j) ,s_ind(i,j)+3-2)
380			& *lagran(3,xG(i,j,bi,bj),px_ind,sp)
381			& +arrayin(e_ind(i,j)+1,s_ind(i,j)+3-2)
382			& *lagran(4,xG(i,j,bi,bj),px_ind,sp)
383			ew_val4 =
384			& arrayin(w_ind(i,j)-1,s_ind(i,j)+4-2)
385			& *lagran(1,xG(i,j,bi,bj),px_ind,sp)
386			& +arrayin(w_ind(i,j) ,s_ind(i,j)+4-2)
387			& *lagran(2,xG(i,j,bi,bj),px_ind,sp)
388			& +arrayin(e_ind(i,j) ,s_ind(i,j)+4-2)
389			& *lagran(3,xG(i,j,bi,bj),px_ind,sp)
390			& +arrayin(e_ind(i,j)+1,s_ind(i,j)+4-2)
391			& *lagran(4,xG(i,j,bi,bj),px_ind,sp)
392			arrayout(i,j,bi,bj)=
393			& +ew_val1*lagran(1,yG(i,j,bi,bj),py_ind,sp)
394			& +ew_val2*lagran(2,yG(i,j,bi,bj),py_ind,sp)
395			& +ew_val3*lagran(3,yG(i,j,bi,bj),py_ind,sp)
396			& +ew_val4*lagran(4,yG(i,j,bi,bj),py_ind,sp)
397			#endif /* TARGET_NEC_SX defined */
398	dimitri	1.2	enddo
399	dimitri	1.1	enddo
400	dimitri	1.2	else
401			stop 'stop in exf_interp.F: interpolation method not supported'
402			endif
403			enddo
404			enddo
405	dimitri	1.1
406	jmc	1.20	RETURN
407	dimitri	1.1	END