/[MITgcm]/MITgcm/pkg/exf/exf_interp.F
ViewVC logotype

Annotation of /MITgcm/pkg/exf/exf_interp.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.3 - (hide annotations) (download)
Thu Oct 9 04:19:19 2003 UTC (20 years, 8 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint51o_pre, checkpoint51n_pre, checkpoint51l_post, checkpoint51q_post, branch-netcdf, checkpoint52d_pre, checkpoint51r_post, checkpoint52b_pre, checkpoint51o_post, checkpoint51p_post, checkpoint52a_pre, checkpoint51i_post, checkpoint52, checkpoint52a_post, checkpoint52b_post, checkpoint52c_post, checkpoint51l_pre, ecco_c52_e35, checkpoint51t_post, checkpoint51n_post, checkpoint51u_post, checkpoint51m_post, checkpoint51s_post
Branch point for: branch-nonh, tg2-branch, checkpoint51n_branch
Changes since 1.2: +1 -1 lines 
 o first check-in for the "branch-genmake2" merge
 o verification suite as run on shelley (gcc 3.2.2):

Wed Oct  8 23:42:29 EDT 2003
                T           S           U           V
G D M    c        m  s        m  s        m  s        m  s
E p a R  g  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .
N n k u  2  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d
2 d e n  d  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .

OPTFILE=NONE

Y Y Y Y 13 16 16 16  0 16 16 16 16 16 16 16 16 13 12  0  0 pass  adjustment.128x64x1
Y Y Y Y 16 16 16 16  0 16 16 16 16 16 16  0  0 16 16  0  0 pass  adjustment.cs-32x32x1
Y Y Y Y 16 16 16 16  0 16 16 16 16 16 16 22  0 16 16 22  0 pass  adjust_nlfs.cs-32x32x1
Y Y Y Y -- 13 13 16 16 13 13 13 13 16 16 16 16 16 16 16 16 N/O   advect_cs
Y Y Y Y -- 22 16 16 16 16 16 16 13 16 16 16 16 16 16 16 16 N/O   advect_xy
Y Y Y Y -- 13 16 13 16 16 16 16 16 16 16 22 16 16 16 16 16 N/O   advect_xz
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  aim.5l_cs
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 16 16 16 16 16 13 16 pass  aim.5l_Equatorial_Channel
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 13 16 16 13 13 16 pass  aim.5l_LatLon
Y Y Y Y 13 16 16 16 16 16 16 16 16 16 13 12 13 13 16 13 16 pass  exp0
Y Y Y Y 14 16 16 16 16 16 16 16 22 16 16 16 13 16 16 22 16 pass  exp1
Y Y Y Y 13 13 16 13 16 16 16 16 16 13 13 16 16 13 13 13 13 pass  exp2
Y Y Y Y 16 16 16 16 16 16 16 16 22 16 16 16 16 16 16 16 16 pass  exp4
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 22 16 16 16 22 16 pass  exp5
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  front_relax
Y Y Y Y 14 16 16 13 13 16 16 13 13 16 13 13 16 12 13 13 16 pass  global_ocean.90x40x15
Y Y Y Y 10 16 16 13 13 16 13 16 16 13 13 13 13 16 16 13 16 FAIL  global_ocean.cs32x15
Y Y Y Y  6 11 12 13 13 12 13 16 13  9  9  9  9 10  9  9 11 FAIL  global_ocean_pressure
Y Y Y Y 14 16 16 13 16 16 16 13 13 13 13 13 16 12 16 13 16 pass  global_with_exf
Y Y Y Y 14 16 16 16 16 16 16 16 16 11 13 22 13 16 16  9 16 pass  hs94.128x64x5
Y Y Y Y 13 16 16 16 16 16 16 16 16 11 16 16 16 13 16 22 13 pass  hs94.1x64x5
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 16 13 13 16 16 22 13 pass  hs94.cs-32x32x5
Y Y Y Y 10 10 16 13 13 16 16 16 22 16 13 13 13 13 13 22 13 FAIL  ideal_2D_oce
Y Y Y Y  8 16 16 16 16 16 16 16 16 13 13  8 16 16 16 16 16 FAIL  internal_wave
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 13 22 13 13 13 22 16 pass  inverted_barometer
Y Y Y Y 12 16 16 16 16 16 16 16 16 16 13 12 13 13 13 13 13 FAIL  lab_sea
Y Y Y Y 11 16 16 16 16 16 16 16 13 13 13 12 13 16 13 12 13 FAIL  natl_box
Y Y Y Y 16 16 16 16 16 16 16 16 22 16 16 16 16 16 16 16 16 pass  plume_on_slope
Y Y Y Y 13 16 16 16 16 13 16 16 16 16 16 16 16 13 16 16 16 pass  solid-body.cs-32x32x1
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     C
45 dimitri 1.2 C infile = name of the input file (direct access binary)
46     C filePrec = file precicision (currently not used, assumes real*4)
47 dimitri 1.1 C arrout = output arrays (different for each processor)
48     C irecord = record number in global file
49     C xG,yG = coordinates for output grid
50     C lon_0, lat_0 = lon and lat of sw corner of global input grid
51     C lon_inc = scalar x-grid increment
52     C lat_inc = vector y-grid increments
53     C nx_in, ny_in = input x-grid and y-grid size
54 dimitri 1.2 C method = 1 for bilinear 2 for bicubic
55     C mythid = thread id
56 dimitri 1.1 C
57    
58     #include "SIZE.h"
59     #include "EEPARAMS.h"
60 dimitri 1.2
61     C subroutine variables
62     character*(*) infile
63     integer filePrec, irecord, nx_in, ny_in
64     _RL arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
65     _RS xG (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
66     _RS yG (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
67     _RL lon_0, lon_inc
68     _RL lat_0, lat_inc(ny_in-1)
69     integer method, mythid
70 dimitri 1.1
71     C local variables
72 dimitri 1.2 integer ierr
73     real*8 ne_fac,nw_fac,se_fac,sw_fac
74     integer e_ind(snx),w_ind(snx)
75     integer n_ind(sny),s_ind(sny)
76     real*8 px_ind(4), py_ind(4), ew_val(4)
77 dimitri 1.1 external lagran
78 dimitri 1.2 real*8 lagran
79     real*4 arrayin(-1:nx_in+2 , -1:ny_in+2)
80     real*8 x_in (-1:nx_in+2), y_in(-1:ny_in+2)
81     integer i, j, k, l, js, bi, bj, sp, interp_unit
82    
83     C check input arguments
84     if ( .NOT. (filePrec .EQ. 32) )
85     & stop 'stop in exf_interp.F: value of filePrec not allowed'
86 dimitri 1.1
87     C read in input data
88     call mdsfindunit( interp_unit, mythid)
89     open(interp_unit,file=infile,status='old',access='direct',
90     & recl=nx_in*ny_in*4)
91     read(interp_unit,rec=irecord) ((arrayin(i,j),i=1,nx_in),j=1,ny_in)
92     #ifdef _BYTESWAPIO
93     call MDS_BYTESWAPR4((nx_in+4)*(ny_in+4), arrayin )
94     #endif
95     close(interp_unit)
96    
97     C setup input grid
98     do i=-1,nx_in+2
99     x_in(i) = lon_0 + (i-1.)*lon_inc
100     enddo
101     y_in(0) = lat_0 - lat_inc(1)
102     y_in(-1)= lat_0 - 2.*lat_inc(1)
103     y_in(1) = lat_0
104     do j=2,ny_in
105     y_in(j) = y_in(j-1) + lat_inc(j-1)
106     enddo
107     y_in(ny_in+1) = y_in(ny_in) + lat_inc(ny_in-1)
108     y_in(ny_in+2) = y_in(ny_in) + 2.*lat_inc(ny_in-1)
109    
110     C enlarge boundary
111 dimitri 1.2 do j=1,ny_in
112 dimitri 1.1 arrayin(0,j) = arrayin(nx_in,j)
113     arrayin(-1,j) = arrayin(nx_in-1,j)
114     arrayin(nx_in+1,j) = arrayin(1,j)
115     arrayin(nx_in+2,j) = arrayin(2,j)
116 dimitri 1.2 enddo
117     do i=-1,nx_in+2
118 dimitri 1.1 arrayin(i,0) = arrayin(i,1)
119     arrayin(i,-1) = arrayin(i,1)
120     arrayin(i,ny_in+1) = arrayin(i,ny_in)
121     arrayin(i,ny_in+2) = arrayin(i,ny_in)
122 dimitri 1.2 enddo
123    
124     do bj = mybylo(mythid), mybyhi(mythid)
125     do bi = mybxlo(mythid), mybxhi(mythid)
126    
127     C check validity of input/output coordinates
128     if ( xG(1,1 ,bi,bj) .le. x_in(0) .or.
129     & xG(snx,1,bi,bj) .ge. x_in(nx_in+1) .or.
130     & yG(1,1 ,bi,bj) .lt. y_in(1) .or.
131     & yG(1,sny,bi,bj) .gt. y_in(ny_in) ) then
132     print*,'ERROR in S/R EXF_INTERP:'
133     print*,' input grid must encompass output grid.'
134     STOP ' ABNORMAL END: S/R EXF_INTERP'
135     endif
136 dimitri 1.1
137     C compute interpolation indices
138     do i=1,snx
139 dimitri 1.2 if (xG(i,1,bi,bj)-x_in(1) .ge. 0.) then
140     w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc) + 1
141     else
142     w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc)
143     endif
144     e_ind(i) = w_ind(i) + 1
145 dimitri 1.1 enddo
146     js = ny_in/2
147     do j=1,sny
148 dimitri 1.2 do while (yG(1,j,bi,bj) .lt. y_in(js))
149     js = (js + 1)/2
150     enddo
151     do while (yG(1,j,bi,bj) .ge. y_in(js+1))
152     js = js + 1
153     enddo
154     s_ind(j) = js
155     n_ind(j) = js + 1
156 dimitri 1.1 enddo
157    
158 dimitri 1.2 if (method .eq. 1) then
159 dimitri 1.1
160 dimitri 1.2 C bilinear interpolation
161     sp = 2
162     do j=1,sny
163     do i=1,snx
164 dimitri 1.1 arrayout(i,j,bi,bj) = 0.
165 dimitri 1.2 do l=0,1
166     px_ind(l+1) = x_in(w_ind(i)+l)
167     py_ind(l+1) = y_in(s_ind(j)+l)
168 dimitri 1.1 enddo
169 dimitri 1.2 do k=1,2
170     ew_val(k) = arrayin(w_ind(i),s_ind(j)+k-1)
171     & *lagran(1,xG(i,1,bi,bj),px_ind,sp)
172     & +arrayin(e_ind(i),s_ind(j)+k-1)
173     & *lagran(2,xG(i,1,bi,bj),px_ind,sp)
174     arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
175     & +ew_val(k)*lagran(k,yG(1,j,bi,bj),py_ind,sp)
176 dimitri 1.1 enddo
177     enddo
178     enddo
179 dimitri 1.2 elseif (method .eq. 2) then
180 dimitri 1.1
181 dimitri 1.2 C bicubic interpolation
182     sp = 4
183     do j=1,sny
184     do i=1,snx
185 dimitri 1.1 arrayout(i,j,bi,bj) = 0.
186 dimitri 1.2 do l=-1,2
187     px_ind(l+2) = x_in(w_ind(i)+l)
188     py_ind(l+2) = y_in(s_ind(j)+l)
189 dimitri 1.1 enddo
190 dimitri 1.2 do k=1,4
191     ew_val(k) =
192     & arrayin(w_ind(i)-1,s_ind(j)+k-2)
193     & *lagran(1,xG(i,1,bi,bj),px_ind,sp)
194     & +arrayin(w_ind(i) ,s_ind(j)+k-2)
195     & *lagran(2,xG(i,1,bi,bj),px_ind,sp)
196     & +arrayin(e_ind(i) ,s_ind(j)+k-2)
197     & *lagran(3,xG(i,1,bi,bj),px_ind,sp)
198     & +arrayin(e_ind(i)+1,s_ind(j)+k-2)
199     & *lagran(4,xG(i,1,bi,bj),px_ind,sp)
200     arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
201     & +ew_val(k)*lagran(k,yG(1,j,bi,bj),py_ind,sp)
202 dimitri 1.1 enddo
203 dimitri 1.2 enddo
204 dimitri 1.1 enddo
205 dimitri 1.2 else
206     stop 'stop in exf_interp.F: interpolation method not supported'
207     endif
208     enddo
209     enddo
210 dimitri 1.1
211     END
  ViewVC Help
Powered by ViewVC 1.1.22