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

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

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


Revision 1.32 - (show annotations) (download)
Fri Mar 13 14:16:36 2015 UTC (9 years, 8 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint65k
Changes since 1.31: +76 -28 lines
add code to allow input grids with latitude starting in the north
(when j=1 corresponds to northern edge of field)

1 C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_interp.F,v 1.31 2014/06/05 15:37:46 jmc Exp $
2 C $Name: $
3
4 #include "EXF_OPTIONS.h"
5
6 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7
8 CBOP
9 C !ROUTINE: EXF_INTERP
10 C !INTERFACE:
11 SUBROUTINE EXF_INTERP(
12 I inFile, filePrec,
13 O arrayout,
14 I irecord, xG_in, yG,
15 I lon_0, lon_inc, lat_0, lat_inc,
16 I nxIn, nyIn, method, myIter, myThid )
17
18 C !DESCRIPTION: \bv
19 C *==========================================================*
20 C | SUBROUTINE EXF_INTERP
21 C | o Load from file a regular lat-lon input field
22 C | and interpolate on to the model grid location
23 C *==========================================================*
24 C \ev
25
26 C !USES:
27 IMPLICIT NONE
28 C === Global variables ===
29 #include "SIZE.h"
30 #include "EEPARAMS.h"
31 #include "PARAMS.h"
32 #ifdef ALLOW_DEBUG
33 # include "EXF_PARAM.h"
34 #endif
35
36 C !INPUT/OUTPUT PARAMETERS:
37 C inFile (string) :: name of the binary input file (direct access)
38 C filePrec (integer) :: number of bits per word in file (32 or 64)
39 C arrayout ( _RL ) :: output array
40 C irecord (integer) :: record number to read
41 C xG_in,yG :: coordinates for output grid to interpolate to
42 C lon_0, lat_0 :: lon and lat of sw corner of global input grid
43 C lon_inc :: scalar x-grid increment
44 C lat_inc :: vector y-grid increments
45 C nxIn,nyIn (integer) :: size in x & y direction of input file to read
46 C method :: 1,11,21 for bilinear; 2,12,22 for bicubic
47 C :: 1,2 for tracer; 11,12 for U; 21,22 for V
48 C myIter (integer) :: current iteration number
49 C myThid (integer) :: My Thread Id number
50
51 CHARACTER*(*) inFile
52 INTEGER filePrec, irecord, nxIn, nyIn
53 _RL arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54 _RS xG_in (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55 _RS yG (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
56 _RL lon_0, lon_inc
57 c _RL lat_0, lat_inc(nyIn-1)
58 _RL lat_0, lat_inc(*)
59 INTEGER method, myIter, myThid
60
61 C !FUNCTIONS:
62 #ifdef ALLOW_DEBUG
63 INTEGER ILNBLNK
64 EXTERNAL ILNBLNK
65 #endif
66
67 C !LOCAL VARIABLES:
68 C arrayin :: input field array (loaded from file)
69 C x_in :: longitude vector defining input field grid
70 C y_in :: latitude vector defining input field grid
71 C w_ind :: input field longitudinal index, on western side of model grid pt
72 C s_ind :: input field latitudinal index, on southern side of model grid pt
73 C bi, bj :: tile indices
74 C i, j, k, l :: loop indices
75 C msgBuf :: Informational/error message buffer
76 _RL arrayin( -1:nxIn+2, -1:nyIn+2 )
77 _RL x_in(-1:nxIn+2), y_in(-1:nyIn+2)
78 INTEGER w_ind(sNx,sNy), s_ind(sNx,sNy)
79 INTEGER bi, bj
80 INTEGER i, j, k, l
81 INTEGER nLoop
82 _RL tmpVar
83 _RS xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
84 _RS threeSixtyRS
85 _RL yPole, symSign, poleValue
86 PARAMETER ( threeSixtyRS = 360. )
87 PARAMETER ( yPole = 90. )
88 INTEGER nxd2
89 LOGICAL xIsPeriodic, poleSymmetry, latFromS2N
90 #ifdef ALLOW_DEBUG
91 LOGICAL debugFlag, gridNotRight
92 CHARACTER*(MAX_LEN_MBUF) msgBuf
93 _RS prtPole(-1:4)
94 #endif
95 CEOP
96
97 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
98 C--- Load inut field
99
100 CALL EXF_INTERP_READ(
101 I inFile, filePrec,
102 O arrayin,
103 I irecord, nxIn, nyIn, myThid )
104
105 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
106 C--- Prepare input grid and input field
107
108 C-- setup input longitude grid
109 DO i=-1,nxIn+2
110 x_in(i) = lon_0 + (i-1)*lon_inc
111 ENDDO
112 xIsPeriodic = nxIn.EQ.NINT( threeSixtyRS / lon_inc )
113 nxd2 = NINT( nxIn*0.5 )
114 poleSymmetry = xIsPeriodic .AND. ( nxIn.EQ.2*nxd2 )
115 #ifdef EXF_USE_OLD_INTERP_POLE
116 poleSymmetry = .FALSE.
117 #endif
118
119 C-- setup input latitude grid
120 y_in(1) = lat_0
121 DO j=1,nyIn+1
122 i = MIN(j,nyIn-1)
123 y_in(j+1) = y_in(j) + lat_inc(i)
124 ENDDO
125 y_in(0) = y_in(1) - lat_inc(1)
126 y_in(-1)= y_in(0) - lat_inc(1)
127 #ifdef ALLOW_DEBUG
128 DO l=-1,4
129 prtPole(l) = 0.
130 ENDDO
131 #endif
132 C-- For tracer (method=1,2) add 1 row @ the pole (if last row is not)
133 C and will fill it with the polarmost-row zonally averaged value.
134 C For vector component, cannot do much without the other component;
135 C averaging properly done if uvInterp=T in S/R EXF_INTERP_UV
136 #ifdef EXF_USE_OLD_INTERP_POLE
137 IF ( .TRUE. ) THEN
138 #else
139 IF ( method.LT.10 ) THEN
140 C- Add 2 row @ southern end; if one is beyond S.pole, put one @ S.pole
141 j = 0
142 IF ( ABS(y_in(j+1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
143 y_in(j) = -yPole
144 y_in(j-1) = -2.*yPole - y_in(1)
145 #ifdef ALLOW_DEBUG
146 prtPole(j) = 1.
147 prtPole(j-1) = 2.
148 #endif
149 ENDIF
150 j = -1
151 IF ( ABS(y_in(j+1)).GT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
152 y_in(j) = -yPole
153 #ifdef ALLOW_DEBUG
154 prtPole(j) = 1.
155 #endif
156 ENDIF
157 #endif /* EXF_USE_OLD_INTERP_POLE */
158 C- Add 2 row @ northern end; if one is beyond N.pole, put one @ N.pole
159 j = nyIn+1
160 IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
161 y_in(j) = yPole
162 y_in(j+1) = 2.*yPole - y_in(j-1)
163 #ifdef ALLOW_DEBUG
164 prtPole(3) = 1.
165 prtPole(3+1) = 2.
166 #endif
167 ENDIF
168 j = nyIn+2
169 IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
170 y_in(j) = yPole
171 #ifdef ALLOW_DEBUG
172 prtPole(4) = 1.
173 #endif
174 ENDIF
175 ENDIF
176 C store orientation of latitude
177 latFromS2N = y_in(0) .LT. y_in(nyIn+1)
178
179 C-- Enlarge boundary
180 IF ( xIsPeriodic ) THEN
181 C- fill-in added column assuming periodicity
182 DO j=1,nyIn
183 arrayin( 0,j) = arrayin(nxIn ,j)
184 arrayin(-1,j) = arrayin(nxIn-1,j)
185 arrayin(nxIn+1,j) = arrayin(1,j)
186 arrayin(nxIn+2,j) = arrayin(2,j)
187 ENDDO
188 ELSE
189 C- fill-in added column from nearest column
190 DO j=1,nyIn
191 arrayin( 0,j) = arrayin(1,j)
192 arrayin(-1,j) = arrayin(1,j)
193 arrayin(nxIn+1,j) = arrayin(nxIn,j)
194 arrayin(nxIn+2,j) = arrayin(nxIn,j)
195 ENDDO
196 ENDIF
197 symSign = 1. _d 0
198 IF ( method.GE.10 ) symSign = -1. _d 0
199 DO l=-1,2
200 j = l
201 IF ( l.GE.1 ) j = nyIn+l
202 k = MAX(1,MIN(j,nyIn))
203 IF ( poleSymmetry .AND. ABS(y_in(j)).GT.yPole ) THEN
204 C- fill-in added row assuming pole-symmetry
205 DO i=-1,nxd2
206 arrayin(i,j) = symSign*arrayin(i+nxd2,k)
207 ENDDO
208 DO i=1,nxd2+2
209 arrayin(i+nxd2,j) = symSign*arrayin(i,k)
210 ENDDO
211 #ifdef ALLOW_DEBUG
212 i = l + 2*( (l+1)/2 )
213 prtPole(i) = prtPole(i) + 0.2
214 #endif
215 ELSE
216 C- fill-in added row from nearest column values
217 DO i=-1,nxIn+2
218 arrayin(i,j) = arrayin(i,k)
219 ENDDO
220 ENDIF
221 ENDDO
222
223 C-- For tracer (method=1,2) set to northernmost zonal-mean value
224 C at 90N to avoid sharp zonal gradients near the Pole.
225 C For vector component, cannot do much without the other component;
226 C averaging properly done if uvInterp=T in S/R EXF_INTERP_UV
227 #ifdef EXF_USE_OLD_INTERP_POLE
228 IF ( .TRUE. ) THEN
229 DO l= 3,4
230 #else
231 IF ( method.LT.10 ) THEN
232 DO l=-1,4
233 #endif
234 j = l
235 IF ( l.GE.2 ) j = nyIn+l-2
236 IF ( ABS(y_in(j)).EQ.yPole ) THEN
237 IF (method.EQ.1 .OR. method.EQ.2) THEN
238 poleValue = 0.
239 DO i=1,nxIn
240 poleValue = poleValue + arrayin(i,j)
241 ENDDO
242 poleValue = poleValue / nxIn
243 DO i=-1,nxIn+2
244 arrayin(i,j) = poleValue
245 ENDDO
246 #ifdef ALLOW_DEBUG
247 prtPole(l) = prtPole(l) + 0.1
248 #endif
249 #ifdef EXF_USE_OLD_INTERP_POLE
250 ELSEIF (method.EQ.11 .OR. method.EQ.12) THEN
251 DO i=-1,nxIn+2
252 arrayin(i,j) = 0.
253 ENDDO
254 #ifdef ALLOW_DEBUG
255 prtPole(l) = prtPole(l) + 0.9
256 #endif
257 #endif /* EXF_USE_OLD_INTERP_POLE */
258 ENDIF
259 ENDIF
260 ENDDO
261 ENDIF
262
263 #ifdef ALLOW_DEBUG
264 debugFlag = ( exf_debugLev.GE.debLevC )
265 & .OR. ( exf_debugLev.GE.debLevB .AND. myIter.LE.nIter0 )
266 C prtPole(l)=0 : extended, =1 : changed to pole, =2 : changed to symetric
267 IF ( debugFlag ) THEN
268 l = ILNBLNK(inFile)
269 _BEGIN_MASTER(myThid)
270 WRITE(msgBuf,'(3A,I6,A,2L5)')
271 & 'EXF_INTERP: file="',inFile(1:l),'", rec=', irecord,
272 & ' , x-Per,P.Sym=', xIsPeriodic, poleSymmetry
273 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
274 & SQUEEZE_RIGHT, myThid )
275 WRITE(msgBuf,'(2A,3F4.1,A,3F12.6)') 'S.edge (j=-1,0,1) :',
276 & ' proc=', (prtPole(j),j=-1,1), ', yIn=', (y_in(j),j=-1,1)
277 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
278 & SQUEEZE_RIGHT, myThid )
279 WRITE(msgBuf,'(2A,3F4.1,A,3F12.6)') 'N.edge (j=+0,+1,+2)',
280 & ' proc=', (prtPole(j),j=2,4), ', yIn=',(y_in(j),j=nyIn,nyIn+2)
281 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
282 & SQUEEZE_RIGHT, myThid )
283 _END_MASTER(myThid)
284 ENDIF
285 #endif /* ALLOW_DEBUG */
286
287 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
288 C--- Prepare output grid and interpolate for each tile
289
290 C-- put xG in interval [ lon_0 , lon_0+360 [
291 DO bj=myByLo(myThid),myByHi(myThid)
292 DO bi=myBxLo(myThid),myBxHi(myThid)
293 DO j=1-OLy,sNy+OLy
294 DO i=1-OLx,sNx+OLx
295 xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
296 & + threeSixtyRS*2.
297 xG(i,j,bi,bj) = lon_0+MOD(xG(i,j,bi,bj),threeSixtyRS)
298 ENDDO
299 ENDDO
300 #ifdef ALLOW_DEBUG
301 C-- Check validity of input/output coordinates
302 IF ( debugFlag ) THEN
303 DO j=1,sNy
304 DO i=1,sNx
305 IF ( latFromS2N ) THEN
306 C y_in(0) is the southern edge
307 gridNotRight = (
308 & xG(i,j,bi,bj) .LT. x_in(0) .OR.
309 & xG(i,j,bi,bj) .GE. x_in(nxIn+1) .OR.
310 & yG(i,j,bi,bj) .LT. y_in(0) .OR.
311 & yG(i,j,bi,bj) .GE. y_in(nyIn+1) )
312 ELSE
313 C y_in(0) is the northern edge
314 gridNotRight = (
315 & xG(i,j,bi,bj) .LT. x_in(0) .OR.
316 & xG(i,j,bi,bj) .GE. x_in(nxIn+1) .OR.
317 & yG(i,j,bi,bj) .GE. y_in(0) .OR.
318 & yG(i,j,bi,bj) .LT. y_in(nyIn+1) )
319 ENDIF
320 IF ( gridNotRight ) THEN
321 l = ILNBLNK(inFile)
322 WRITE(msgBuf,'(3A,I6)')
323 & 'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord
324 CALL PRINT_ERROR( msgBuf, myThid )
325 WRITE(msgBuf,'(A)')
326 & 'EXF_INTERP: input grid must encompass output grid.'
327 CALL PRINT_ERROR( msgBuf, myThid )
328 WRITE(msgBuf,'(A,2I8,2I6,A,1P2E14.6)') 'i,j,bi,bj=',
329 & i,j,bi,bj, ' , xG,yG=', xG(i,j,bi,bj), yG(i,j,bi,bj)
330 CALL PRINT_ERROR( msgBuf, myThid )
331 WRITE(msgBuf,'(A,I9,A,1P2E14.6)') 'nxIn=', nxIn,
332 & ' , x_in(0,nxIn+1)=', x_in(0) ,x_in(nxIn+1)
333 CALL PRINT_ERROR( msgBuf, myThid )
334 WRITE(msgBuf,'(A,I9,A,1P2E14.6)') 'nyIn=', nyIn,
335 & ' , y_in(0,nyIn+1)=', y_in(0) ,y_in(nyIn+1)
336 CALL PRINT_ERROR( msgBuf, myThid )
337 STOP 'ABNORMAL END: S/R EXF_INTERP'
338 ENDIF
339 ENDDO
340 ENDDO
341 ENDIF
342 #endif /* ALLOW_DEBUG */
343 ENDDO
344 ENDDO
345
346 C-- Compute interpolation lon & lat index mapping
347 C # of pts = nyIn+2 ; # of interval = nyIn+1 ; evaluate nLoop as
348 C 1 + truncated log2(# interval -1); add epsil=1.e-3 for safey
349 tmpVar = nyIn + 1. _d -3
350 nLoop = 1 + INT( LOG(tmpVar)/LOG(2. _d 0) )
351 C
352 DO bj = myByLo(myThid), myByHi(myThid)
353 DO bi = myBxLo(myThid), myBxHi(myThid)
354
355 IF ( latFromS2N ) THEN
356 C-- Case 1: y_in(0) is the southern edge
357 C-- latitude index
358 DO j=1,sNy
359 DO i=1,sNx
360 s_ind(i,j) = 0
361 w_ind(i,j) = nyIn+1
362 ENDDO
363 ENDDO
364 DO l=1,nLoop
365 DO j=1,sNy
366 DO i=1,sNx
367 IF ( w_ind(i,j).GT.s_ind(i,j)+1 ) THEN
368 k = NINT( (s_ind(i,j)+w_ind(i,j))*0.5 )
369 IF ( yG(i,j,bi,bj) .LT. y_in(k) ) THEN
370 w_ind(i,j) = k
371 ELSE
372 s_ind(i,j) = k
373 ENDIF
374 ENDIF
375 ENDDO
376 ENDDO
377 ENDDO
378 ELSE
379 C-- Case 2: y_in(0) is the northern edge
380 C-- In this case the roles of s_ind and w_ind are exchanged: now
381 C-- s_ind is the index on the northern side grid box
382 C-- latitude index
383 DO j=1,sNy
384 DO i=1,sNx
385 s_ind(i,j) = nyIn+1
386 w_ind(i,j) = 0
387 ENDDO
388 ENDDO
389 DO l=1,nLoop
390 DO j=1,sNy
391 DO i=1,sNx
392 IF ( s_ind(i,j).GT.w_ind(i,j)+1 ) THEN
393 k = NINT( (s_ind(i,j)+w_ind(i,j))*0.5 )
394 IF ( yG(i,j+1,bi,bj) .GT. y_in(k) ) THEN
395 s_ind(i,j) = k
396 ELSE
397 w_ind(i,j) = k
398 ENDIF
399 ENDIF
400 ENDDO
401 ENDDO
402 ENDDO
403 ENDIF
404 #ifdef ALLOW_DEBUG
405 IF ( debugFlag ) THEN
406 C- Check that we found the right lat. index
407 DO j=1,sNy
408 DO i=1,sNx
409 IF ( latFromS2N ) THEN
410 gridNotRight = w_ind(i,j).NE.s_ind(i,j)+1
411 ELSE
412 gridNotRight = s_ind(i,j).NE.w_ind(i,j)+1
413 ENDIF
414 IF ( gridNotRight ) THEN
415 l = ILNBLNK(inFile)
416 WRITE(msgBuf,'(3A,I4,A,I4)')
417 & 'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord,
418 & ', nLoop=', nLoop
419 CALL PRINT_ERROR( msgBuf, myThid )
420 WRITE(msgBuf,'(A)')
421 & 'EXF_INTERP: did not find Latitude index for grid-pt:'
422 CALL PRINT_ERROR( msgBuf, myThid )
423 WRITE(msgBuf,'(A,2I8,2I6,A,1PE16.8)')
424 & 'EXF_INTERP: i,j,bi,bj=',i,j,bi,bj,' , yG=',yG(i,j,bi,bj)
425 CALL PRINT_ERROR( msgBuf, myThid )
426 WRITE(msgBuf,'(A,I8,A,1PE16.8)')
427 & 'EXF_INTERP: s_ind=',s_ind(i,j),', lat=',y_in(s_ind(i,j))
428 CALL PRINT_ERROR( msgBuf, myThid )
429 WRITE(msgBuf,'(A,I8,A,1PE16.8)')
430 & 'EXF_INTERP: n_ind=',w_ind(i,j),', lat=',y_in(w_ind(i,j))
431 CALL PRINT_ERROR( msgBuf, myThid )
432 STOP 'ABNORMAL END: S/R EXF_INTERP'
433 ENDIF
434 ENDDO
435 ENDDO
436 ENDIF
437 #endif /* ALLOW_DEBUG */
438 C-- longitude index
439 DO j=1,sNy
440 DO i=1,sNx
441 w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(-1))/lon_inc) - 1
442 ENDDO
443 ENDDO
444
445 C-- Do interpolation using lon & lat index mapping
446 CALL EXF_INTERPOLATE(
447 I inFile, irecord, method,
448 I nxIn, nyIn, x_in, y_in,
449 O arrayin,
450 O arrayout,
451 I xG, yG,
452 I w_ind, s_ind,
453 I bi, bj, myThid )
454
455 ENDDO
456 ENDDO
457
458 RETURN
459 END

  ViewVC Help
Powered by ViewVC 1.1.22