/[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.30 - (show annotations) (download)
Mon Apr 2 21:24:54 2012 UTC (12 years, 2 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64x, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint64
Changes since 1.29: +3 -1 lines
avoid multiple print to the same file (STDOUT) if multi-threads

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

  ViewVC Help
Powered by ViewVC 1.1.22