Line data Source code
1 : module rgrd_mod
2 : implicit none
3 : private
4 :
5 : public :: rgrd2
6 :
7 : contains
8 :
9 : !
10 : !
11 : ! ... file rgrd2.f
12 :
13 : ! this file contains documentation for subroutine rgrd2 followed by
14 : ! fortran code for rgrd2 and additional subroutines.
15 :
16 : ! ... author
17 :
18 : ! John C. Adams (NCAR 1997)
19 :
20 : ! ... subroutine rgrd2(nx,ny,x,y,p,mx,my,xx,yy,q,intpol,w,lw,iw,liw,ier)
21 :
22 : ! ... purpose
23 :
24 : ! subroutine rgrd2 interpolates the values p(i,j) on the orthogonal
25 : ! grid (x(i),y(j)) for i=1,...,nx and j=1,...,ny onto q(ii,jj) on the
26 : ! orthogonal grid (xx(ii),yy(jj)) for ii=1,...,mx and jj=1,...,my.
27 :
28 : ! ... language
29 :
30 : ! coded in portable FORTRAN77 and FORTRAN90
31 :
32 : ! ... test program
33 :
34 : ! file trgrd2.f on regridpack includes a test program for subroutine rgrd2
35 :
36 : ! ... method
37 :
38 : ! linear or cubic interpolation is used (independently) in
39 : ! each direction (see argument intpol).
40 :
41 : ! ... required files
42 :
43 : ! file rgrd1.f must be loaded with rgrd2.f. it includes
44 : ! subroutines called by the routines in rgrd2.f
45 :
46 : ! ... requirements
47 :
48 : ! each of the x,y grids must be strictly montonically increasing
49 : ! and each of the xx,yy grids must be montonically increasing (see
50 : ! ier = 4). in addition the (X,Y) region
51 :
52 : ! [xx(1),xx(mx)] X [yy(1),yy(my)]
53 :
54 : ! must lie within the (X,Y) region
55 :
56 : ! [x(1),x(nx)] X [y(1),y(ny)].
57 :
58 : ! extrapolation is not allowed (see ier=3). if these (X,Y)
59 : ! regions are identical and the orthogonal grids are UNIFORM
60 : ! in each direction then subroutine rgrd2u (see file rgrd2u.f)
61 : ! should be used instead of rgrd2.
62 :
63 : ! ... efficiency
64 :
65 : ! inner most loops in regridpack software vectorize.
66 : ! If the arguments mx,my (see below) have different values, optimal
67 : ! vectorization will be achieved if mx > my.
68 :
69 :
70 : ! *** input argument
71 :
72 :
73 : ! ... nx
74 :
75 : ! the integer dimension of the grid vector x and the first dimension
76 : ! of p. nx > 1 if intpol(1) = 1 or nx > 3 if intpol(1) = 3 is required.
77 :
78 : ! ... ny
79 :
80 : ! the integer dimension of the grid vector y and the second dimension
81 : ! of p. ny > 1 if intpol(2) = 1 or ny > 3 if intpol(2) = 3 is required.
82 :
83 : ! ... x
84 :
85 : ! a real nx vector of strictly increasing values which defines the x
86 : ! portion of the orthogonal grid on which p is given
87 :
88 : ! ... y
89 :
90 : ! a real ny vector of strictly increasing values which defines the y
91 : ! portion of the orthogonal grid on which p is given
92 :
93 : ! ... p
94 :
95 : ! a real nx by ny array of values given on the orthogonal (x,y) grid
96 :
97 : ! ... mx
98 :
99 : ! the integer dimension of the grid vector xx and the first dimension
100 : ! of q. mx > 0 is required.
101 :
102 : ! ... my
103 :
104 : ! the integer dimension of the grid vector yy and the second dimension
105 : ! of q. my > 0 is required.
106 :
107 : ! ... xx
108 :
109 : ! a real mx vector of increasing values which defines the x portion of the
110 : ! orthogonal grid on which q is defined. xx(1) < x(1) or xx(mx) > x(nx)
111 : ! is not allowed (see ier = 3)
112 :
113 : ! ... yy
114 :
115 : ! a real my vector of increasing values which defines the y portion of the
116 : ! orthogonal grid on which q is defined. yy(1) < y(1) or yy(my) > y(ny)
117 : ! is not allowed (see ier = 3)
118 :
119 : ! ... intpol
120 :
121 : ! an integer vector of dimension 2 which sets linear or cubic
122 : ! interpolation in the x,y directions as follows:
123 :
124 : ! intpol(1) = 1 sets linear interpolation in the x direction
125 : ! intpol(1) = 3 sets cubic interpolation in the x direction.
126 :
127 : ! intpol(2) = 1 sets linear interpolation in the y direction
128 : ! intpol(2) = 3 sets cubic interpolation in the y direction.
129 :
130 : ! values other than 1 or 3 in intpol are not allowed (ier = 5).
131 :
132 : ! ... w
133 :
134 : ! a real work space of length at least lw which must be provided in the
135 : ! routine calling rgrd2
136 :
137 : ! ... lw
138 :
139 : ! the integer length of the real work space w. let
140 :
141 : ! lwx = mx if intpol(1) = 1
142 : ! lwx = 4*mx if intpol(1) = 3
143 :
144 : ! lwy = my+2*mx if intpol(2) = 1
145 : ! lwy = 4*(mx+my) if intpol(2) = 3
146 :
147 : ! then lw must be greater than or equal to lwx+lwy
148 :
149 : ! ... iw
150 :
151 : ! an integer work space of length at least liw which must be provided in the
152 : ! routine calling rgrd2
153 :
154 : ! ... liw
155 :
156 : ! the integer length of the integer work space iw. liw must be at least mx+my
157 :
158 : ! *** output arguments
159 :
160 :
161 : ! ... q
162 :
163 : ! a real mx by my array of values on the (xx,yy) grid which are
164 : ! interpolated from p on the (x,y) grid
165 :
166 : ! ... ier
167 :
168 : ! an integer error flag set as follows:
169 :
170 : ! ier = 0 if no errors in input arguments are detected
171 :
172 : ! ier = 1 if min0(mx,my) < 1
173 :
174 : ! ier = 2 if nx < 2 when intpol(1)=1 or nx < 4 when intpol(1)=3 (or)
175 : ! ny < 2 when intpol(2)=1 or ny < 4 when intpol(2)=3
176 :
177 : ! ier = 3 if xx(1) < x(1) or x(nx) < xx(mx) (or)
178 : ! yy(1) < y(1) or y(ny) < yy(my) (or)
179 :
180 : ! *** to avoid this flag when end points are intended to be the
181 : ! same but may differ slightly due to roundoff error, they
182 : ! should be set exactly in the calling routine (e.g., if both
183 : ! grids have the same y boundaries then yy(1)=y(1) and yy(my)=y(ny)
184 : ! should be set before calling rgrd2)
185 :
186 : ! ier = 4 if one of the grids x,y is not strictly monotonically
187 : ! increasing or if one of the grids xx,yy is not
188 : ! montonically increasing. more precisely if:
189 :
190 : ! x(i+1) <= x(i) for some i such that 1 <= i < nx (or)
191 :
192 : ! y(j+1) <= y(j) for some j such that 1 <= j < ny (or)
193 :
194 : ! xx(ii+1) < xx(ii) for some ii such that 1 <= ii < mx (or)
195 :
196 : ! yy(jj+1) < yy(jj) for some jj such that 1 <= jj < my
197 :
198 : ! ier = 5 if lw or liw is to small (insufficient work space)
199 :
200 : ! ier = 6 if intpol(1) or intpol(2) is not equal to 1 or 3
201 :
202 : ! ************************************************************************
203 :
204 : ! end of rgrd2 documentation, fortran code follows:
205 :
206 : ! ************************************************************************
207 :
208 0 : subroutine rgrd2(nx,ny,x,y,p,mx,my,xx,yy,q,intpol,w,lw,iw,liw,ier)
209 :
210 : use shr_kind_mod,only: r8 => shr_kind_r8
211 :
212 : implicit none
213 :
214 : integer, intent(in) :: nx
215 : integer, intent(in) :: ny
216 : real(r8), intent(in) :: x(nx)
217 : real(r8), intent(in) :: y(ny)
218 : real(r8), intent(inout) :: p(nx,ny)
219 : integer, intent(in) :: mx
220 : integer, intent(in) :: my
221 : real(r8), intent(in) :: xx(mx)
222 : real(r8), intent(in) :: yy(my)
223 : real(r8), intent(out) :: q(mx,my)
224 : integer, intent(in) :: intpol(2)
225 : integer, intent(inout) :: lw
226 : real(r8), intent(inout) :: w(lw)
227 : integer, intent(inout) :: liw
228 : integer, intent(inout) :: iw(liw)
229 : integer, intent(out) :: ier
230 :
231 :
232 : integer :: i,ii,j,jj,j2,j3,j4,j5,j6,j7,j8,j9,i2,i3,i4,i5
233 : integer :: jy,lwx,lwy
234 :
235 :
236 : ! check input arguments
237 :
238 0 : ier = 1
239 :
240 : ! check (xx,yy) grid resolution
241 :
242 0 : IF (MIN0(mx,my) < 1) RETURN
243 :
244 : ! check intpol
245 :
246 0 : ier = 6
247 0 : IF (intpol(1) /= 1 .AND. intpol(1) /= 3) RETURN
248 0 : IF (intpol(2) /= 1 .AND. intpol(2) /= 3) RETURN
249 : ! write(6,"('rgrd2: nx,ny,mx,my = ',6i6)") nx,ny,mx,my,intpol
250 : ! write(6,"('rgrd2: alatm = ',/,(6e12.4))") y
251 : ! write(6,"('rgrd2: ylatm = ',/,(6e12.4))") yy
252 : ! write(6,"('rgrd2: ylonm = ',/,(6e12.4))") xx
253 : ! write(6,"('rgrd2: potm(1,:) = ',/,(8e12.4))") p(1,:)
254 :
255 : ! check (x,y) grid resolution
256 :
257 0 : ier = 2
258 0 : IF (intpol(1) == 1 .AND. nx < 2) RETURN
259 0 : IF (intpol(1) == 3 .AND. nx < 4) RETURN
260 0 : IF (intpol(2) == 1 .AND. ny < 2) RETURN
261 0 : IF (intpol(2) == 3 .AND. ny < 4) RETURN
262 :
263 : ! check work space lengths
264 :
265 0 : ier = 5
266 0 : IF (intpol(1) == 1) THEN
267 : lwx = mx
268 : ELSE
269 0 : lwx = 4*mx
270 : END IF
271 0 : IF (intpol(2) == 1) THEN
272 0 : lwy = my+2*mx
273 : ELSE
274 0 : lwy = 4*(mx+my)
275 : END IF
276 0 : IF (lw < lwx+lwy) RETURN
277 0 : IF (liw < mx+my) RETURN
278 :
279 : ! check (xx,yy) grid contained in (x,y) grid
280 :
281 0 : ier = 3
282 0 : IF (xx(1) < x(1) .OR. xx(mx) > x(nx)) RETURN
283 0 : IF (yy(1) < y(1) .OR. yy(my) > y(ny)) RETURN
284 :
285 : ! check montonicity of grids
286 :
287 0 : ier = 4
288 0 : DO i=2,nx
289 0 : IF (x(i-1) >= x(i)) RETURN
290 : END DO
291 0 : DO j=2,ny
292 0 : IF (y(j-1) >= y(j)) RETURN
293 : END DO
294 0 : DO ii=2,mx
295 0 : IF (xx(ii-1) > xx(ii)) RETURN
296 : END DO
297 0 : DO jj=2,my
298 0 : IF (yy(jj-1) > yy(jj)) RETURN
299 : END DO
300 :
301 : ! arguments o.k.
302 :
303 0 : ier = 0
304 :
305 : ! set pointer in integer work space
306 :
307 0 : jy = mx+1
308 0 : IF (intpol(2) == 1) THEN
309 :
310 : ! linearly interpolate in y
311 :
312 0 : j2 = 1
313 0 : j3 = j2
314 0 : j4 = j3+my
315 0 : j5 = j4
316 0 : j6 = j5
317 0 : j7 = j6
318 0 : j8 = j7+mx
319 0 : j9 = j8+mx
320 :
321 : ! set y interpolation indices and scales and linearly interpolate
322 :
323 0 : CALL linmx(ny,y,my,yy,iw(jy),w(j3))
324 0 : i2 = j9
325 :
326 : ! set work space portion and indices which depend on x interpolation
327 :
328 0 : IF (intpol(1) == 1) THEN
329 0 : i3 = i2
330 0 : i4 = i3
331 0 : i5 = i4
332 0 : CALL linmx(nx,x,mx,xx,iw,w(i3))
333 : ELSE
334 0 : i3 = i2+mx
335 0 : i4 = i3+mx
336 0 : i5 = i4+mx
337 0 : CALL cubnmx(nx,x,mx,xx,iw,w(i2),w(i3),w(i4),w(i5))
338 : END IF
339 : ! write(6,"('rgrd2: w(i2),w(i3),w(i4),w(i5) = ',4e12.4)")
340 : ! + w(i2),w(i3),w(i4),w(i5)
341 : CALL lint2(nx,ny,p,mx,my,q,intpol,iw(jy),w(j3), &
342 0 : w(j7),w(j8),iw,w(i2),w(i3),w(i4),w(i5))
343 :
344 : ! write(6,"('rgrd2: potout1 = ',/,(6e12.4))") q
345 :
346 0 : RETURN
347 :
348 : ELSE
349 :
350 : ! cubically interpolate in y, set indice pointers
351 :
352 0 : j2 = 1
353 0 : j3 = j2+my
354 0 : j4 = j3+my
355 0 : j5 = j4+my
356 0 : j6 = j5+my
357 0 : j7 = j6+mx
358 0 : j8 = j7+mx
359 0 : j9 = j8+mx
360 0 : CALL cubnmx(ny,y,my,yy,iw(jy),w(j2),w(j3),w(j4),w(j5))
361 0 : i2 = j9+mx
362 :
363 : ! set work space portion and indices which depend on x interpolation
364 :
365 0 : IF (intpol(1) == 1) THEN
366 0 : i3 = i2
367 0 : i4 = i3
368 0 : i5 = i4
369 0 : CALL linmx(nx,x,mx,xx,iw,w(i3))
370 : ELSE
371 0 : i3 = i2+mx
372 0 : i4 = i3+mx
373 0 : i5 = i4+mx
374 0 : CALL cubnmx(nx,x,mx,xx,iw,w(i2),w(i3),w(i4),w(i5))
375 : END IF
376 : CALL cubt2(nx,ny,p,mx,my,q,intpol,iw(jy),w(j2),w(j3), &
377 0 : w(j4),w(j5),w(j6),w(j7),w(j8),w(j9),iw,w(i2),w(i3),w(i4),w(i5))
378 :
379 : ! write(6,"('rgrd2: potout2 = ',/,(6e12.4))") q
380 :
381 0 : RETURN
382 : END IF
383 : END SUBROUTINE rgrd2
384 : !----------------------------------------------------------------------------
385 :
386 0 : subroutine lint2(nx,ny,p,mx,my,q,intpol,jy,dy,pj,pjp, ix,dxm,dx,dxp,dxpp)
387 :
388 : use shr_kind_mod,only: r8 => shr_kind_r8
389 :
390 : implicit none
391 :
392 : integer, intent(in) :: nx
393 : integer, intent(in) :: ny
394 : real(r8), intent(inout) :: p(nx,ny)
395 : integer, intent(in) :: mx
396 : integer, intent(in) :: my
397 : real(r8), intent(out) :: q(mx,my)
398 : integer, intent(in) :: intpol(2)
399 : integer, intent(in) :: jy(my)
400 : real(r8), intent(in) :: dy(my)
401 : real(r8), intent(out) :: pj(mx)
402 : real(r8), intent(inout) :: pjp(mx)
403 : integer, intent(inout) :: ix(mx)
404 : real(r8), intent(inout) :: dxm(mx)
405 : real(r8), intent(inout) :: dx(mx)
406 : real(r8), intent(inout) :: dxp(mx)
407 : real(r8), intent(inout) :: dxpp(mx)
408 :
409 : integer :: jsave,j,jj,ii
410 :
411 :
412 :
413 : ! linearly interpolate in y
414 :
415 0 : IF (intpol(1) == 1) THEN
416 :
417 : ! linear in x
418 :
419 : jsave = -1
420 0 : DO jj=1,my
421 0 : j = jy(jj)
422 0 : IF (j == jsave) THEN
423 :
424 : ! j pointer has not moved since last pass (no updates or interpolation)
425 :
426 0 : ELSE IF (j == jsave+1) THEN
427 :
428 : ! update j and interpolate j+1
429 :
430 0 : DO ii=1,mx
431 0 : pj(ii) = pjp(ii)
432 : END DO
433 0 : CALL lint1(nx,p(1,j+1),mx,pjp,ix,dx)
434 : ELSE
435 :
436 : ! interpolate j,j+1in pj,pjp on xx mesh
437 :
438 0 : CALL lint1(nx,p(1,j),mx,pj,ix,dx)
439 0 : CALL lint1(nx,p(1,j+1),mx,pjp,ix,dx)
440 : END IF
441 :
442 : ! save j pointer for next pass
443 :
444 0 : jsave = j
445 :
446 : ! linearly interpolate q(ii,jj) from pjp,pj in y direction
447 :
448 0 : DO ii=1,mx
449 0 : q(ii,jj) = pj(ii)+dy(jj)*(pjp(ii)-pj(ii))
450 : END DO
451 : END DO
452 :
453 : ELSE
454 :
455 : ! cubic in x
456 :
457 : jsave = -1
458 0 : DO jj=1,my
459 0 : j = jy(jj)
460 0 : IF (j == jsave) THEN
461 :
462 : ! j pointer has not moved since last pass (no updates or interpolation)
463 :
464 0 : ELSE IF (j == jsave+1) THEN
465 :
466 : ! update j and interpolate j+1
467 :
468 0 : DO ii=1,mx
469 0 : pj(ii) = pjp(ii)
470 : END DO
471 0 : CALL cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
472 : ELSE
473 :
474 : ! interpolate j,j+1 in pj,pjp on xx mesh
475 :
476 0 : CALL cubt1(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp)
477 0 : CALL cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
478 : END IF
479 :
480 : ! save j pointer for next pass
481 :
482 0 : jsave = j
483 :
484 : ! linearly interpolate q(ii,jj) from pjp,pj in y direction
485 :
486 0 : DO ii=1,mx
487 0 : q(ii,jj) = pj(ii)+dy(jj)*(pjp(ii)-pj(ii))
488 : END DO
489 : END DO
490 : ! write(6,"('lint2: potm = ',/,(6e12.4))") p
491 : ! write(6,"('lint2: potout3 = ',/,(6e12.4))") q
492 : RETURN
493 : END IF
494 : end subroutine lint2
495 :
496 0 : subroutine cubt2(nx,ny,p,mx,my,q,intpol,jy,dym,dy,dyp, &
497 0 : dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
498 :
499 : use shr_kind_mod,only: r8 => shr_kind_r8
500 :
501 : implicit none
502 :
503 : integer, intent(in) :: nx
504 : integer, intent(in) :: ny
505 : real(r8), intent(inout) :: p(nx,ny)
506 : integer, intent(in) :: mx
507 : integer, intent(in) :: my
508 : real(r8), intent(out) :: q(mx,my)
509 : integer, intent(in) :: intpol(2)
510 : integer, intent(in) :: jy(my)
511 : real(r8), intent(in) :: dym(my)
512 : real(r8), intent(in) :: dy(my)
513 : real(r8), intent(in) :: dyp(my)
514 : real(r8), intent(inout) :: dypp(my)
515 : real(r8), intent(out) :: pjm(mx)
516 : real(r8), intent(inout) :: pj(mx)
517 : real(r8), intent(inout) :: pjp(mx)
518 : real(r8), intent(inout) :: pjpp(mx)
519 : integer, intent(inout) :: ix(mx)
520 : real(r8), intent(inout) :: dxm(mx)
521 : real(r8), intent(inout) :: dx(mx)
522 : real(r8), intent(inout) :: dxp(mx)
523 : real(r8), intent(inout) :: dxpp(mx)
524 :
525 : integer :: jsave,j,jj,ii
526 :
527 :
528 :
529 :
530 0 : IF (intpol(1) == 1) THEN
531 :
532 : ! linear in x
533 :
534 : jsave = -3
535 0 : DO jj=1,my
536 :
537 : ! load closest four j lines containing interpolate on xx mesh
538 : ! for j-1,j,j+1,j+2 in pjm,pj,pjp,pjpp
539 :
540 0 : j = jy(jj)
541 0 : IF (j == jsave) THEN
542 :
543 : ! j pointer has not moved since last pass (no updates or interpolation)
544 :
545 0 : ELSE IF (j == jsave+1) THEN
546 :
547 : ! update j-1,j,j+1 and interpolate j+2
548 :
549 0 : DO ii=1,mx
550 0 : pjm(ii) = pj(ii)
551 0 : pj(ii) = pjp(ii)
552 0 : pjp(ii) = pjpp(ii)
553 : END DO
554 0 : CALL lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
555 0 : ELSE IF (j == jsave+2) THEN
556 :
557 : ! update j-1,j and interpolate j+1,j+2
558 :
559 0 : DO ii=1,mx
560 0 : pjm(ii) = pjp(ii)
561 0 : pj(ii) = pjpp(ii)
562 : END DO
563 0 : CALL lint1(nx,p(1,j+1),mx,pjp,ix,dx)
564 0 : CALL lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
565 0 : ELSE IF (j == jsave+3) THEN
566 :
567 : ! update j-1 and interpolate j,j+1,j+2
568 :
569 0 : DO ii=1,mx
570 0 : pjm(ii) = pjpp(ii)
571 : END DO
572 0 : CALL lint1(nx,p(1,j),mx,pj,ix,dx)
573 0 : CALL lint1(nx,p(1,j+1),mx,pjp,ix,dx)
574 0 : CALL lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
575 : ELSE
576 :
577 : ! interpolate all four j-1,j,j+1,j+2
578 :
579 0 : CALL lint1(nx,p(1,j-1),mx,pjm,ix,dx)
580 0 : CALL lint1(nx,p(1,j),mx,pj,ix,dx)
581 0 : CALL lint1(nx,p(1,j+1),mx,pjp,ix,dx)
582 0 : CALL lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
583 : END IF
584 :
585 : ! save j pointer for next pass
586 :
587 0 : jsave = j
588 :
589 : ! cubically interpolate q(ii,jj) from pjm,pj,pjp,pjpp in y direction
590 :
591 0 : DO ii=1,mx
592 0 : q(ii,jj) = dym(jj)*pjm(ii)+dy(jj)*pj(ii)+dyp(jj)*pjp(ii)+ &
593 0 : dypp(jj)*pjpp(ii)
594 : END DO
595 : END DO
596 : RETURN
597 :
598 : ELSE
599 :
600 : ! cubic in x
601 :
602 : jsave = -3
603 0 : DO jj=1,my
604 :
605 : ! load closest four j lines containing interpolate on xx mesh
606 : ! for j-1,j,j+1,j+2 in pjm,pj,pjp,pjpp
607 :
608 0 : j = jy(jj)
609 0 : IF (j == jsave) THEN
610 :
611 : ! j pointer has not moved since last pass (no updates or interpolation)
612 :
613 0 : ELSE IF (j == jsave+1) THEN
614 :
615 : ! update j-1,j,j+1 and interpolate j+2
616 :
617 0 : DO ii=1,mx
618 0 : pjm(ii) = pj(ii)
619 0 : pj(ii) = pjp(ii)
620 0 : pjp(ii) = pjpp(ii)
621 : END DO
622 0 : CALL cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
623 0 : ELSE IF (j == jsave+2) THEN
624 :
625 : ! update j-1,j and interpolate j+1,j+2
626 :
627 0 : DO ii=1,mx
628 0 : pjm(ii) = pjp(ii)
629 0 : pj(ii) = pjpp(ii)
630 : END DO
631 0 : CALL cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
632 0 : CALL cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
633 0 : ELSE IF (j == jsave+3) THEN
634 :
635 : ! update j-1 and interpolate j,j+1,j+2
636 :
637 0 : DO ii=1,mx
638 0 : pjm(ii) = pjpp(ii)
639 : END DO
640 0 : CALL cubt1(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp)
641 0 : CALL cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
642 0 : CALL cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
643 : ELSE
644 :
645 : ! interpolate all four j-1,j,j+1,j+2
646 :
647 0 : CALL cubt1(nx,p(1,j-1),mx,pjm,ix,dxm,dx,dxp,dxpp)
648 0 : CALL cubt1(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp)
649 0 : CALL cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
650 0 : CALL cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
651 : END IF
652 :
653 : ! save j pointer for next pass
654 :
655 0 : jsave = j
656 :
657 : ! cubically interpolate q(ii,jj) from pjm,pj,pjp,pjpp in y direction
658 :
659 0 : DO ii=1,mx
660 0 : q(ii,jj) = dym(jj)*pjm(ii)+dy(jj)*pj(ii)+dyp(jj)*pjp(ii)+ &
661 0 : dypp(jj)*pjpp(ii)
662 : END DO
663 : END DO
664 : RETURN
665 : END IF
666 : end subroutine cubt2
667 : !------------------------------------------------------------------------------
668 :
669 0 : subroutine lint1(nx,p,mx,q,ix,dx)
670 :
671 : use shr_kind_mod,only: r8 => shr_kind_r8
672 :
673 : implicit none
674 :
675 : integer, intent(in) :: nx
676 : real(r8), intent(in) :: p(nx)
677 : integer, intent(in) :: mx
678 : real(r8), intent(out) :: q(mx)
679 : integer, intent(in) :: ix(mx)
680 : real(r8), intent(in) :: dx(mx)
681 :
682 : integer :: ii,i
683 :
684 :
685 : ! linearly interpolate p on x onto q on xx
686 :
687 0 : DO ii=1,mx
688 0 : i = ix(ii)
689 0 : q(ii) = p(i)+dx(ii)*(p(i+1)-p(i))
690 : END DO
691 0 : RETURN
692 : end subroutine lint1
693 :
694 0 : subroutine cubt1(nx,p,mx,q,ix,dxm,dx,dxp,dxpp)
695 :
696 : use shr_kind_mod,only: r8 => shr_kind_r8
697 :
698 : implicit none
699 :
700 : integer, intent(in) :: nx
701 : real(r8), intent(in) :: p(nx)
702 : integer, intent(in) :: mx
703 : real(r8), intent(out) :: q(mx)
704 : integer, intent(in) :: ix(mx)
705 : real(r8), intent(in) :: dxm(mx)
706 : real(r8), intent(in) :: dx(mx)
707 : real(r8), intent(in) :: dxp(mx)
708 : real(r8), intent(inout) :: dxpp(mx)
709 :
710 : integer :: i,ii
711 :
712 :
713 : ! cubically interpolate p on x to q on xx
714 :
715 0 : DO ii=1,mx
716 0 : i = ix(ii)
717 0 : q(ii) = dxm(ii)*p(i-1)+dx(ii)*p(i)+dxp(ii)*p(i+1) &
718 0 : +dxpp(ii)*p(i+2)
719 : END DO
720 0 : RETURN
721 : end subroutine cubt1
722 :
723 0 : subroutine cubnmx(nx,x,mx,xx,ix,dxm,dx,dxp,dxpp)
724 :
725 : use shr_kind_mod,only: r8 => shr_kind_r8
726 :
727 : implicit none
728 :
729 : integer, intent(in) :: nx
730 : real(r8), intent(in) :: x(*)
731 : integer, intent(in) :: mx
732 : real(r8), intent(in) :: xx(*)
733 : integer, intent(out) :: ix(*)
734 : real(r8), intent(out) :: dxm(*)
735 : real(r8), intent(out) :: dx(*)
736 : real(r8), intent(out) :: dxp(*)
737 : real(r8), intent(out) :: dxpp(*)
738 :
739 : integer :: i,ii,isrt
740 :
741 0 : isrt = 1
742 0 : DO ii=1,mx
743 :
744 : ! set i in [2,nx-2] closest s.t.
745 : ! x(i-1),x(i),x(i+1),x(i+2) can interpolate xx(ii)
746 :
747 0 : DO i=isrt,nx-1
748 0 : IF (x(i+1) >= xx(ii)) THEN
749 0 : ix(ii) = MIN0(nx-2,MAX0(2,i))
750 0 : isrt = ix(ii)
751 0 : GO TO 3
752 : END IF
753 : END DO
754 0 : 3 CONTINUE
755 : END DO
756 :
757 : ! set cubic scale terms
758 :
759 0 : DO ii=1,mx
760 0 : i = ix(ii)
761 0 : dxm(ii) = (xx(ii)-x(i))*(xx(ii)-x(i+1))*(xx(ii)-x(i+2))/ &
762 0 : ((x(i-1)-x(i))*(x(i-1)-x(i+1))*(x(i-1)-x(i+2)))
763 : dx(ii) = (xx(ii)-x(i-1))*(xx(ii)-x(i+1))*(xx(ii)-x(i+2))/ &
764 0 : ((x(i)-x(i-1))*(x(i)-x(i+1))*(x(i)-x(i+2)))
765 : dxp(ii) = (xx(ii)-x(i-1))*(xx(ii)-x(i))*(xx(ii)-x(i+2))/ &
766 0 : ((x(i+1)-x(i-1))*(x(i+1)-x(i))*(x(i+1)-x(i+2)))
767 : dxpp(ii) = (xx(ii)-x(i-1))*(xx(ii)-x(i))*(xx(ii)-x(i+1))/ &
768 0 : ((x(i+2)-x(i-1))*(x(i+2)-x(i))*(x(i+2)-x(i+1)))
769 : END DO
770 0 : RETURN
771 : end subroutine cubnmx
772 : !------------------------------------------------------------------------------
773 :
774 0 : subroutine linmx(nx,x,mx,xx,ix,dx)
775 :
776 : ! set x grid pointers for xx grid and interpolation scale terms
777 :
778 : use shr_kind_mod,only: r8 => shr_kind_r8
779 :
780 : implicit none
781 :
782 : integer, intent(in) :: nx
783 : real(r8), intent(in) :: x(*)
784 : integer, intent(in) :: mx
785 : real(r8), intent(in) :: xx(*)
786 : integer, intent(out) :: ix(*)
787 : real(r8), intent(out) :: dx(*)
788 :
789 : integer :: isrt,ii,i
790 :
791 0 : isrt = 1
792 0 : do ii=1,mx
793 :
794 : ! find x(i) s.t. x(i) < xx(ii) <= x(i+1)
795 :
796 0 : DO i=isrt,nx-1
797 0 : IF (x(i+1) >= xx(ii)) THEN
798 0 : isrt = i
799 0 : ix(ii) = i
800 0 : GO TO 3
801 : END IF
802 : END DO
803 0 : 3 CONTINUE
804 : END DO
805 :
806 : ! set linear scale term
807 :
808 0 : DO ii=1,mx
809 0 : i = ix(ii)
810 0 : dx(ii) = (xx(ii)-x(i))/(x(i+1)-x(i))
811 : END DO
812 0 : RETURN
813 : end subroutine linmx
814 : !------------------------------------------------------------------------------
815 : end module rgrd_mod
|
