Line data Source code
1 : module interpolate_data
2 : ! Description:
3 : ! Routines for interpolation of data in latitude, longitude and time.
4 : !
5 : ! Author: Gathered from a number of places and put into the current format by Jim Edwards
6 : !
7 : ! Modules Used:
8 : !
9 : use shr_kind_mod, only: r8 => shr_kind_r8
10 : use cam_abortutils, only: endrun
11 : use cam_logfile, only: iulog
12 : implicit none
13 : private
14 : !
15 : ! Public Methods:
16 : !
17 :
18 : public :: interp_type, lininterp, vertinterp, bilin, get_timeinterp_factors
19 : public :: lininterp_init, lininterp_finish
20 : type interp_type
21 : real(r8), pointer :: wgts(:)
22 : real(r8), pointer :: wgtn(:)
23 : integer, pointer :: jjm(:)
24 : integer, pointer :: jjp(:)
25 : end type interp_type
26 : interface lininterp
27 : module procedure lininterp_original
28 : module procedure lininterp_full1d
29 : module procedure lininterp1d
30 : module procedure lininterp2d2d
31 : module procedure lininterp2d1d
32 : module procedure lininterp3d2d
33 : end interface
34 :
35 : integer, parameter, public :: extrap_method_zero = 0
36 : integer, parameter, public :: extrap_method_bndry = 1
37 : integer, parameter, public :: extrap_method_cycle = 2
38 :
39 : contains
40 1536 : subroutine lininterp_full1d (arrin, yin, nin, arrout, yout, nout)
41 : integer, intent(in) :: nin, nout
42 : real(r8), intent(in) :: arrin(nin), yin(nin), yout(nout)
43 : real(r8), intent(out) :: arrout(nout)
44 : type (interp_type) :: interp_wgts
45 :
46 1536 : call lininterp_init(yin, nin, yout, nout, extrap_method_bndry, interp_wgts)
47 1536 : call lininterp1d(arrin, nin, arrout, nout, interp_wgts)
48 1536 : call lininterp_finish(interp_wgts)
49 :
50 1536 : end subroutine lininterp_full1d
51 :
52 645180856 : subroutine lininterp_init(yin, nin, yout, nout, extrap_method, interp_wgts, &
53 : cyclicmin, cyclicmax)
54 : !
55 : ! Description:
56 : ! Initialize a variable of type(interp_type) with weights for linear interpolation.
57 : ! this variable can then be used in calls to lininterp1d and lininterp2d.
58 : ! yin is a 1d array of length nin of locations to interpolate from - this array must
59 : ! be monotonic but can be increasing or decreasing
60 : ! yout is a 1d array of length nout of locations to interpolate to, this array need
61 : ! not be ordered
62 : ! extrap_method determines how to handle yout points beyond the bounds of yin
63 : ! if 0 set values outside output grid to 0
64 : ! if 1 set to boundary value
65 : ! if 2 set to cyclic boundaries
66 : ! optional values cyclicmin and cyclicmax can be used to set the bounds of the
67 : ! cyclic mapping - these default to 0 and 360.
68 : !
69 :
70 : integer, intent(in) :: nin
71 : integer, intent(in) :: nout
72 : real(r8), intent(in) :: yin(:) ! input mesh
73 : real(r8), intent(in) :: yout(:) ! output mesh
74 : integer, intent(in) :: extrap_method ! if 0 set values outside output grid to 0
75 : ! if 1 set to boundary value
76 : ! if 2 set to cyclic boundaries
77 : real(r8), intent(in), optional :: cyclicmin, cyclicmax
78 :
79 : type (interp_type), intent(out) :: interp_wgts
80 : real(r8) :: cmin, cmax
81 : real(r8) :: extrap
82 : real(r8) :: dyinwrap
83 : real(r8) :: ratio
84 : real(r8) :: avgdyin
85 : integer :: i, j, icount
86 : integer :: jj
87 645180856 : real(r8), pointer :: wgts(:)
88 645180856 : real(r8), pointer :: wgtn(:)
89 645180856 : integer, pointer :: jjm(:)
90 645180856 : integer, pointer :: jjp(:)
91 : logical :: increasing
92 : !
93 : ! Check validity of input coordinate arrays: must be monotonically increasing,
94 : ! and have a total of at least 2 elements
95 : !
96 645180856 : if (nin.lt.2) then
97 0 : call endrun('LININTERP: Must have at least 2 input points for interpolation')
98 : end if
99 645180856 : if(present(cyclicmin)) then
100 352944 : cmin=cyclicmin
101 : else
102 : cmin=0_r8
103 : end if
104 645180856 : if(present(cyclicmax)) then
105 352944 : cmax=cyclicmax
106 : else
107 : cmax=360_r8
108 : end if
109 645180856 : if(cmax<=cmin) then
110 0 : call endrun('LININTERP: cyclic min value must be < max value')
111 : end if
112 645180856 : increasing=.true.
113 645180856 : icount = 0
114 >12269*10^7 : do j=1,nin-1
115 >12269*10^7 : if (yin(j).gt.yin(j+1)) icount = icount + 1
116 : end do
117 645180856 : if(icount.eq.nin-1) then
118 133410616 : increasing = .false.
119 : icount=0
120 : endif
121 511770240 : if (icount.gt.0) then
122 0 : call endrun('LININTERP: Non-monotonic input coordinate array found')
123 : end if
124 : allocate(interp_wgts%jjm(nout), &
125 : interp_wgts%jjp(nout), &
126 : interp_wgts%wgts(nout), &
127 4516265992 : interp_wgts%wgtn(nout))
128 :
129 645180856 : jjm => interp_wgts%jjm
130 645180856 : jjp => interp_wgts%jjp
131 645180856 : wgts => interp_wgts%wgts
132 645180856 : wgtn => interp_wgts%wgtn
133 :
134 : !
135 : ! Initialize index arrays for later checking
136 : !
137 1305065840 : jjm = 0
138 1305065840 : jjp = 0
139 :
140 645180856 : extrap = 0._r8
141 645180856 : select case (extrap_method)
142 : case (extrap_method_zero)
143 : !
144 : ! For values which extend beyond boundaries, set weights
145 : ! such that values will be 0.
146 : !
147 0 : do j=1,nout
148 0 : if(increasing) then
149 0 : if (yout(j).lt.yin(1)) then
150 0 : jjm(j) = 1
151 0 : jjp(j) = 1
152 0 : wgts(j) = 0._r8
153 0 : wgtn(j) = 0._r8
154 0 : extrap = extrap + 1._r8
155 0 : else if (yout(j).gt.yin(nin)) then
156 0 : jjm(j) = nin
157 0 : jjp(j) = nin
158 0 : wgts(j) = 0._r8
159 0 : wgtn(j) = 0._r8
160 0 : extrap = extrap + 1._r8
161 : end if
162 : else
163 0 : if (yout(j).gt.yin(1)) then
164 0 : jjm(j) = 1
165 0 : jjp(j) = 1
166 0 : wgts(j) = 0._r8
167 0 : wgtn(j) = 0._r8
168 0 : extrap = extrap + 1._r8
169 0 : else if (yout(j).lt.yin(nin)) then
170 0 : jjm(j) = nin
171 0 : jjp(j) = nin
172 0 : wgts(j) = 0._r8
173 0 : wgtn(j) = 0._r8
174 0 : extrap = extrap + 1._r8
175 : end if
176 : end if
177 : end do
178 : case (extrap_method_bndry)
179 : !
180 : ! For values which extend beyond boundaries, set weights
181 : ! such that values will just be copied.
182 : !
183 1299172496 : do j=1,nout
184 1299172496 : if(increasing) then
185 520933968 : if (yout(j).le.yin(1)) then
186 8704902 : jjm(j) = 1
187 8704902 : jjp(j) = 1
188 8704902 : wgts(j) = 1._r8
189 8704902 : wgtn(j) = 0._r8
190 8704902 : extrap = extrap + 1._r8
191 512229066 : else if (yout(j).gt.yin(nin)) then
192 48578098 : jjm(j) = nin
193 48578098 : jjp(j) = nin
194 48578098 : wgts(j) = 1._r8
195 48578098 : wgtn(j) = 0._r8
196 48578098 : extrap = extrap + 1._r8
197 : end if
198 : else
199 133410616 : if (yout(j).gt.yin(1)) then
200 0 : jjm(j) = 1
201 0 : jjp(j) = 1
202 0 : wgts(j) = 1._r8
203 0 : wgtn(j) = 0._r8
204 0 : extrap = extrap + 1._r8
205 133410616 : else if (yout(j).le.yin(nin)) then
206 0 : jjm(j) = nin
207 0 : jjp(j) = nin
208 0 : wgts(j) = 1._r8
209 0 : wgtn(j) = 0._r8
210 0 : extrap = extrap + 1._r8
211 : end if
212 : end if
213 : end do
214 : case (extrap_method_cycle)
215 : !
216 : ! For values which extend beyond boundaries, set weights
217 : ! for circular boundaries
218 : !
219 352944 : dyinwrap = yin(1) + (cmax-cmin) - yin(nin)
220 352944 : avgdyin = abs(yin(nin)-yin(1))/(nin-1._r8)
221 352944 : ratio = dyinwrap/avgdyin
222 352944 : if (ratio < 0.9_r8 .or. ratio > 1.1_r8) then
223 0 : write(iulog,*) 'Lininterp: Bad dyinwrap value =',dyinwrap,&
224 0 : ' avg=', avgdyin, yin(1),yin(nin)
225 0 : call endrun('interpolate_data')
226 : end if
227 :
228 651074200 : do j=1,nout
229 5893344 : if(increasing) then
230 5540400 : if (yout(j) <= yin(1)) then
231 0 : jjm(j) = nin
232 0 : jjp(j) = 1
233 0 : wgts(j) = (yin(1)-yout(j))/dyinwrap
234 0 : wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap
235 5540400 : else if (yout(j) > yin(nin)) then
236 31464 : jjm(j) = nin
237 31464 : jjp(j) = 1
238 31464 : wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap
239 31464 : wgtn(j) = (yout(j)-yin(nin))/dyinwrap
240 : end if
241 : else
242 0 : if (yout(j) > yin(1)) then
243 0 : jjm(j) = nin
244 0 : jjp(j) = 1
245 0 : wgts(j) = (yin(1)-yout(j))/dyinwrap
246 0 : wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap
247 0 : else if (yout(j) <= yin(nin)) then
248 0 : jjm(j) = nin
249 0 : jjp(j) = 1
250 0 : wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap
251 0 : wgtn(j) = (yout(j)+(cmax-cmin)-yin(nin))/dyinwrap
252 : end if
253 :
254 : endif
255 : end do
256 : end select
257 :
258 : !
259 : ! Loop though output indices finding input indices and weights
260 : !
261 645180856 : if(increasing) then
262 1038244608 : do j=1,nout
263 44740200290 : do jj=1,nin-1
264 44228430050 : if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then
265 469159904 : jjm(j) = jj
266 469159904 : jjp(j) = jj + 1
267 469159904 : wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj))
268 469159904 : wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj))
269 469159904 : exit
270 : end if
271 : end do
272 : end do
273 : else
274 266821232 : do j=1,nout
275 3815618918 : do jj=1,nin-1
276 3682208302 : if (yout(j).le.yin(jj) .and. yout(j).gt.yin(jj+1)) then
277 133410616 : jjm(j) = jj
278 133410616 : jjp(j) = jj + 1
279 133410616 : wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj))
280 133410616 : wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj))
281 133410616 : exit
282 : end if
283 : end do
284 : end do
285 : end if
286 :
287 : !
288 : ! Check that interp/extrap points have been found for all outputs
289 : !
290 645180856 : icount = 0
291 1305065840 : do j=1,nout
292 659884984 : if (jjm(j).eq.0 .or. jjp(j).eq.0) icount = icount + 1
293 659884984 : ratio=wgts(j)+wgtn(j)
294 1305065840 : if((ratio<0.9_r8.or.ratio>1.1_r8).and.extrap_method.ne.0) then
295 0 : write(iulog,*) j, wgts(j),wgtn(j),jjm(j),jjp(j), increasing,extrap_method
296 0 : call endrun('Bad weight computed in LININTERP_init')
297 : end if
298 : end do
299 645180856 : if (icount.gt.0) then
300 0 : call endrun('LININTERP: Point found without interp indices')
301 : end if
302 :
303 645180856 : end subroutine lininterp_init
304 :
305 17628179928 : subroutine lininterp1d (arrin, n1, arrout, m1, interp_wgts)
306 : !-----------------------------------------------------------------------
307 : !
308 : ! Purpose: Do a linear interpolation from input mesh to output
309 : ! mesh with weights as set in lininterp_init.
310 : !
311 : !
312 : ! Author: Jim Edwards
313 : !
314 : !-----------------------------------------------------------------------
315 : !-----------------------------------------------------------------------
316 : implicit none
317 : !-----------------------------------------------------------------------
318 : !
319 : ! Arguments
320 : !
321 : integer, intent(in) :: n1 ! number of input latitudes
322 : integer, intent(in) :: m1 ! number of output latitudes
323 :
324 : real(r8), intent(in) :: arrin(n1) ! input array of values to interpolate
325 : type(interp_type), intent(in) :: interp_wgts
326 : real(r8), intent(out) :: arrout(m1) ! interpolated array
327 :
328 : !
329 : ! Local workspace
330 : !
331 : integer j ! latitude indices
332 17628179928 : integer, pointer :: jjm(:)
333 17628179928 : integer, pointer :: jjp(:)
334 :
335 17628179928 : real(r8), pointer :: wgts(:)
336 17628179928 : real(r8), pointer :: wgtn(:)
337 :
338 :
339 17628179928 : jjm => interp_wgts%jjm
340 17628179928 : jjp => interp_wgts%jjp
341 17628179928 : wgts => interp_wgts%wgts
342 17628179928 : wgtn => interp_wgts%wgtn
343 :
344 : !
345 : ! Do the interpolation
346 : !
347 35462180784 : do j=1,m1
348 35462180784 : arrout(j) = arrin(jjm(j))*wgts(j) + arrin(jjp(j))*wgtn(j)
349 : end do
350 :
351 17628179928 : return
352 17628179928 : end subroutine lininterp1d
353 :
354 0 : subroutine lininterp2d2d(arrin, n1, n2, arrout, m1, m2, wgt1, wgt2)
355 : implicit none
356 : !-----------------------------------------------------------------------
357 : !
358 : ! Arguments
359 : !
360 : integer, intent(in) :: n1, n2, m1, m2
361 : real(r8), intent(in) :: arrin(n1,n2) ! input array of values to interpolate
362 : type(interp_type), intent(in) :: wgt1, wgt2
363 : real(r8), intent(out) :: arrout(m1,m2) ! interpolated array
364 : !
365 : ! locals
366 : !
367 : integer i,j ! indices
368 0 : integer, pointer :: iim(:), jjm(:)
369 0 : integer, pointer :: iip(:), jjp(:)
370 :
371 0 : real(r8), pointer :: wgts1(:), wgts2(:)
372 0 : real(r8), pointer :: wgtn1(:), wgtn2(:)
373 :
374 0 : real(r8) :: arrtmp(n1,m2)
375 :
376 :
377 0 : jjm => wgt2%jjm
378 0 : jjp => wgt2%jjp
379 0 : wgts2 => wgt2%wgts
380 0 : wgtn2 => wgt2%wgtn
381 :
382 0 : iim => wgt1%jjm
383 0 : iip => wgt1%jjp
384 0 : wgts1 => wgt1%wgts
385 0 : wgtn1 => wgt1%wgtn
386 :
387 0 : do j=1,m2
388 0 : do i=1,n1
389 0 : arrtmp(i,j) = arrin(i,jjm(j))*wgts2(j) + arrin(i,jjp(j))*wgtn2(j)
390 : end do
391 : end do
392 :
393 0 : do j=1,m2
394 0 : do i=1,m1
395 0 : arrout(i,j) = arrtmp(iim(i),j)*wgts1(i) + arrtmp(iip(i),j)*wgtn1(i)
396 : end do
397 : end do
398 :
399 0 : end subroutine lininterp2d2d
400 3869167928 : subroutine lininterp2d1d(arrin, n1, n2, arrout, m1, wgt1, wgt2, fldname)
401 : implicit none
402 : !-----------------------------------------------------------------------
403 : !
404 : ! Arguments
405 : !
406 : integer, intent(in) :: n1, n2, m1
407 : real(r8), intent(in) :: arrin(n1,n2) ! input array of values to interpolate
408 : type(interp_type), intent(in) :: wgt1, wgt2
409 : real(r8), intent(out) :: arrout(m1) ! interpolated array
410 : character(len=*), intent(in), optional :: fldname(:)
411 : !
412 : ! locals
413 : !
414 : integer i ! indices
415 3869167928 : integer, pointer :: iim(:), jjm(:)
416 3869167928 : integer, pointer :: iip(:), jjp(:)
417 :
418 3869167928 : real(r8), pointer :: wgts(:), wgte(:)
419 3869167928 : real(r8), pointer :: wgtn(:), wgtw(:)
420 :
421 3869167928 : jjm => wgt2%jjm
422 3869167928 : jjp => wgt2%jjp
423 3869167928 : wgts => wgt2%wgts
424 3869167928 : wgtn => wgt2%wgtn
425 :
426 3869167928 : iim => wgt1%jjm
427 3869167928 : iip => wgt1%jjp
428 3869167928 : wgtw => wgt1%wgts
429 3869167928 : wgte => wgt1%wgtn
430 :
431 7742158192 : do i=1,m1
432 15491961056 : arrout(i) = arrin(iim(i),jjm(i))*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i))*wgte(i)*wgts(i) + &
433 23234119248 : arrin(iim(i),jjp(i))*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i))*wgte(i)*wgtn(i)
434 : end do
435 :
436 :
437 3869167928 : end subroutine lininterp2d1d
438 160992 : subroutine lininterp3d2d(arrin, n1, n2, n3, arrout, m1, len1, wgt1, wgt2)
439 : implicit none
440 : !-----------------------------------------------------------------------
441 : !
442 : ! Arguments
443 : !
444 : integer, intent(in) :: n1, n2, n3, m1, len1 ! m1 is to len1 as ncols is to pcols
445 : real(r8), intent(in) :: arrin(n1,n2,n3) ! input array of values to interpolate
446 : type(interp_type), intent(in) :: wgt1, wgt2
447 : real(r8), intent(out) :: arrout(len1, n3) ! interpolated array
448 :
449 : !
450 : ! locals
451 : !
452 : integer i, k ! indices
453 160992 : integer, pointer :: iim(:), jjm(:)
454 160992 : integer, pointer :: iip(:), jjp(:)
455 :
456 160992 : real(r8), pointer :: wgts(:), wgte(:)
457 160992 : real(r8), pointer :: wgtn(:), wgtw(:)
458 :
459 160992 : jjm => wgt2%jjm
460 160992 : jjp => wgt2%jjp
461 160992 : wgts => wgt2%wgts
462 160992 : wgtn => wgt2%wgtn
463 :
464 160992 : iim => wgt1%jjm
465 160992 : iip => wgt1%jjp
466 160992 : wgtw => wgt1%wgts
467 160992 : wgte => wgt1%wgtn
468 :
469 4346784 : do k=1,n3
470 70053984 : do i=1,m1
471 328536000 : arrout(i,k) = arrin(iim(i),jjm(i),k)*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i),k)*wgte(i)*wgts(i) + &
472 398428992 : arrin(iim(i),jjp(i),k)*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i),k)*wgte(i)*wgtn(i)
473 : end do
474 : end do
475 :
476 160992 : end subroutine lininterp3d2d
477 :
478 :
479 :
480 :
481 645180856 : subroutine lininterp_finish(interp_wgts)
482 : type(interp_type) :: interp_wgts
483 :
484 0 : deallocate(interp_wgts%jjm, &
485 0 : interp_wgts%jjp, &
486 0 : interp_wgts%wgts, &
487 645180856 : interp_wgts%wgtn)
488 :
489 : nullify(interp_wgts%jjm, &
490 645180856 : interp_wgts%jjp, &
491 645180856 : interp_wgts%wgts, &
492 645180856 : interp_wgts%wgtn)
493 645180856 : end subroutine lininterp_finish
494 :
495 0 : subroutine lininterp_original (arrin, yin, nlev, nlatin, arrout, &
496 0 : yout, nlatout)
497 : !-----------------------------------------------------------------------
498 : !
499 : ! Purpose: Do a linear interpolation from input mesh defined by yin to output
500 : ! mesh defined by yout. Where extrapolation is necessary, values will
501 : ! be copied from the extreme edge of the input grid. Vectorization is over
502 : ! the vertical (nlev) dimension.
503 : !
504 : ! Method: Check validity of input, then determine weights, then do the N-S interpolation.
505 : !
506 : ! Author: Jim Rosinski
507 : ! Modified: Jim Edwards so that there is no requirement of monotonacity on the yout array
508 : !
509 : !-----------------------------------------------------------------------
510 : implicit none
511 : !-----------------------------------------------------------------------
512 : !
513 : ! Arguments
514 : !
515 : integer, intent(in) :: nlev ! number of vertical levels
516 : integer, intent(in) :: nlatin ! number of input latitudes
517 : integer, intent(in) :: nlatout ! number of output latitudes
518 :
519 : real(r8), intent(in) :: arrin(nlev,nlatin) ! input array of values to interpolate
520 : real(r8), intent(in) :: yin(nlatin) ! input mesh
521 : real(r8), intent(in) :: yout(nlatout) ! output mesh
522 :
523 : real(r8), intent(out) :: arrout(nlev,nlatout) ! interpolated array
524 : !
525 : ! Local workspace
526 : !
527 : integer j, jj ! latitude indices
528 : integer jjprev ! latitude indices
529 : integer k ! level index
530 : integer icount ! number of values
531 :
532 : real(r8) extrap ! percent grid non-overlap
533 : !
534 : ! Dynamic
535 : !
536 0 : integer :: jjm(nlatout)
537 0 : integer :: jjp(nlatout)
538 :
539 0 : real(r8) :: wgts(nlatout)
540 0 : real(r8) :: wgtn(nlatout)
541 : !
542 : ! Check validity of input coordinate arrays: must be monotonically increasing,
543 : ! and have a total of at least 2 elements
544 : !
545 0 : if (nlatin.lt.2) then
546 0 : call endrun('LININTERP: Must have at least 2 input points for interpolation')
547 : end if
548 :
549 0 : icount = 0
550 0 : do j=1,nlatin-1
551 0 : if (yin(j).gt.yin(j+1)) icount = icount + 1
552 : end do
553 :
554 :
555 0 : if (icount.gt.0) then
556 0 : call endrun('LININTERP: Non-monotonic coordinate array(s) found')
557 : end if
558 : !
559 : ! Initialize index arrays for later checking
560 : !
561 0 : do j=1,nlatout
562 0 : jjm(j) = 0
563 0 : jjp(j) = 0
564 : end do
565 : !
566 : ! For values which extend beyond N and S boundaries, set weights
567 : ! such that values will just be copied.
568 : !
569 : extrap = 0._r8
570 :
571 0 : do j=1,nlatout
572 0 : if (yout(j).le.yin(1)) then
573 0 : jjm(j) = 1
574 0 : jjp(j) = 1
575 0 : wgts(j) = 1._r8
576 0 : wgtn(j) = 0._r8
577 : extrap=extrap+1._r8
578 0 : else if (yout(j).gt.yin(nlatin)) then
579 0 : jjm(j) = nlatin
580 0 : jjp(j) = nlatin
581 0 : wgts(j) = 1._r8
582 0 : wgtn(j) = 0._r8
583 : extrap=extrap+1._r8
584 : endif
585 : end do
586 :
587 : !
588 : ! Loop though output indices finding input indices and weights
589 : !
590 0 : do j=1,nlatout
591 0 : do jj=1,nlatin-1
592 0 : if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then
593 0 : jjm(j) = jj
594 0 : jjp(j) = jj + 1
595 0 : wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj))
596 0 : wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj))
597 0 : exit
598 : end if
599 : end do
600 : end do
601 : !
602 : ! Check that interp/extrap points have been found for all outputs
603 : !
604 : icount = 0
605 0 : do j=1,nlatout
606 0 : if (jjm(j).eq.0 .or. jjp(j).eq.0) then
607 0 : icount = icount + 1
608 : end if
609 : end do
610 0 : if (icount.gt.0) then
611 0 : call endrun('LININTERP: Point found without interp indices')
612 : end if
613 : !
614 : ! Do the interpolation
615 : !
616 0 : do j=1,nlatout
617 0 : do k=1,nlev
618 0 : arrout(k,j) = arrin(k,jjm(j))*wgts(j) + arrin(k,jjp(j))*wgtn(j)
619 : end do
620 : end do
621 :
622 0 : return
623 : end subroutine lininterp_original
624 :
625 :
626 0 : subroutine bilin (arrin, xin, yin, nlondin, nlonin, &
627 0 : nlevdin, nlev, nlatin, arrout, xout, &
628 0 : yout, nlondout, nlonout, nlevdout, nlatout)
629 : !-----------------------------------------------------------------------
630 : !
631 : ! Purpose:
632 : !
633 : ! Do a bilinear interpolation from input mesh defined by xin, yin to output
634 : ! mesh defined by xout, yout. Circularity is assumed in the x-direction so
635 : ! input x-grid must be in degrees east and must start from Greenwich. When
636 : ! extrapolation is necessary in the N-S direction, values will be copied
637 : ! from the extreme edge of the input grid. Vectorization is over the
638 : ! longitude dimension. Input grid is assumed rectangular. Output grid
639 : ! is assumed ragged, i.e. xout is a 2-d mesh.
640 : !
641 : ! Method: Interpolate first in longitude, then in latitude.
642 : !
643 : ! Author: Jim Rosinski
644 : !
645 : !-----------------------------------------------------------------------
646 : use shr_kind_mod, only: r8 => shr_kind_r8
647 : use cam_abortutils, only: endrun
648 : !-----------------------------------------------------------------------
649 : implicit none
650 : !-----------------------------------------------------------------------
651 : !
652 : ! Input arguments
653 : !
654 : integer, intent(in) :: nlondin ! longitude dimension of input grid
655 : integer, intent(in) :: nlonin ! number of real longitudes (input)
656 : integer, intent(in) :: nlevdin ! vertical dimension of input grid
657 : integer, intent(in) :: nlev ! number of vertical levels
658 : integer, intent(in) :: nlatin ! number of input latitudes
659 : integer, intent(in) :: nlatout ! number of output latitudes
660 : integer, intent(in) :: nlondout ! longitude dimension of output grid
661 : integer, intent(in) :: nlonout(nlatout) ! number of output longitudes per lat
662 : integer, intent(in) :: nlevdout ! vertical dimension of output grid
663 :
664 : real(r8), intent(in) :: arrin(nlondin,nlevdin,nlatin) ! input array of values to interpolate
665 : real(r8), intent(in) :: xin(nlondin) ! input x mesh
666 : real(r8), intent(in) :: yin(nlatin) ! input y mesh
667 : real(r8), intent(in) :: xout(nlondout,nlatout) ! output x mesh
668 : real(r8), intent(in) :: yout(nlatout) ! output y mesh
669 : !
670 : ! Output arguments
671 : !
672 : real(r8), intent(out) :: arrout(nlondout,nlevdout,nlatout) ! interpolated array
673 : !
674 : ! Local workspace
675 : !
676 : integer :: i, ii, iw, ie, iiprev ! longitude indices
677 : integer :: j, jj, js, jn, jjprev ! latitude indices
678 : integer :: k ! level index
679 : integer :: icount ! number of bad values
680 :
681 : real(r8) :: extrap ! percent grid non-overlap
682 : real(r8) :: dxinwrap ! delta-x on input grid for 2-pi
683 : real(r8) :: avgdxin ! avg input delta-x
684 : real(r8) :: ratio ! compare dxinwrap to avgdxin
685 : real(r8) :: sum ! sum of weights (used for testing)
686 : !
687 : ! Dynamic
688 : !
689 0 : integer :: iim(nlondout) ! interpolation index to the left
690 0 : integer :: iip(nlondout) ! interpolation index to the right
691 0 : integer :: jjm(nlatout) ! interpolation index to the south
692 0 : integer :: jjp(nlatout) ! interpolation index to the north
693 :
694 0 : real(r8) :: wgts(nlatout) ! interpolation weight to the north
695 0 : real(r8) :: wgtn(nlatout) ! interpolation weight to the north
696 0 : real(r8) :: wgte(nlondout) ! interpolation weight to the north
697 0 : real(r8) :: wgtw(nlondout) ! interpolation weight to the north
698 0 : real(r8) :: igrid(nlonin) ! interpolation weight to the north
699 : !
700 : ! Check validity of input coordinate arrays: must be monotonically increasing,
701 : ! and have a total of at least 2 elements
702 : !
703 0 : if (nlonin < 2 .or. nlatin < 2) then
704 0 : call endrun ('BILIN: Must have at least 2 input points for interpolation')
705 : end if
706 :
707 0 : if (xin(1) < 0._r8 .or. xin(nlonin) > 360._r8) then
708 0 : call endrun ('BILIN: Input x-grid must be between 0 and 360')
709 : end if
710 :
711 0 : icount = 0
712 0 : do i=1,nlonin-1
713 0 : if (xin(i) >= xin(i+1)) icount = icount + 1
714 : end do
715 :
716 0 : do j=1,nlatin-1
717 0 : if (yin(j) >= yin(j+1)) icount = icount + 1
718 : end do
719 :
720 0 : do j=1,nlatout-1
721 0 : if (yout(j) >= yout(j+1)) icount = icount + 1
722 : end do
723 :
724 0 : do j=1,nlatout
725 0 : do i=1,nlonout(j)-1
726 0 : if (xout(i,j) >= xout(i+1,j)) icount = icount + 1
727 : end do
728 : end do
729 :
730 0 : if (icount > 0) then
731 0 : call endrun ('BILIN: Non-monotonic coordinate array(s) found')
732 : end if
733 :
734 0 : if (yout(nlatout) <= yin(1) .or. yout(1) >= yin(nlatin)) then
735 0 : call endrun ('BILIN: No overlap in y-coordinates')
736 : end if
737 :
738 0 : do j=1,nlatout
739 0 : if (xout(1,j) < 0._r8 .or. xout(nlonout(j),j) > 360._r8) then
740 0 : call endrun ('BILIN: Output x-grid must be between 0 and 360')
741 : end if
742 :
743 0 : if (xout(nlonout(j),j) <= xin(1) .or. &
744 0 : xout(1,j) >= xin(nlonin)) then
745 0 : call endrun ('BILIN: No overlap in x-coordinates')
746 : end if
747 : end do
748 : !
749 : ! Initialize index arrays for later checking
750 : !
751 0 : do j=1,nlatout
752 0 : jjm(j) = 0
753 0 : jjp(j) = 0
754 : end do
755 : !
756 : ! For values which extend beyond N and S boundaries, set weights
757 : ! such that values will just be copied.
758 : !
759 0 : do js=1,nlatout
760 0 : if (yout(js) > yin(1)) exit
761 0 : jjm(js) = 1
762 0 : jjp(js) = 1
763 0 : wgts(js) = 1._r8
764 0 : wgtn(js) = 0._r8
765 : end do
766 :
767 0 : do jn=nlatout,1,-1
768 0 : if (yout(jn) <= yin(nlatin)) exit
769 0 : jjm(jn) = nlatin
770 0 : jjp(jn) = nlatin
771 0 : wgts(jn) = 1._r8
772 0 : wgtn(jn) = 0._r8
773 : end do
774 : !
775 : ! Loop though output indices finding input indices and weights
776 : !
777 0 : jjprev = 1
778 0 : do j=js,jn
779 0 : do jj=jjprev,nlatin-1
780 0 : if (yout(j) > yin(jj) .and. yout(j) <= yin(jj+1)) then
781 0 : jjm(j) = jj
782 0 : jjp(j) = jj + 1
783 0 : wgts(j) = (yin(jj+1) - yout(j)) / (yin(jj+1) - yin(jj))
784 0 : wgtn(j) = (yout(j) - yin(jj)) / (yin(jj+1) - yin(jj))
785 0 : goto 30
786 : end if
787 : end do
788 0 : call endrun ('BILIN: Failed to find interp values')
789 0 : 30 jjprev = jj
790 : end do
791 :
792 0 : dxinwrap = xin(1) + 360._r8 - xin(nlonin)
793 : !
794 : ! Check for sane dxinwrap values. Allow to differ no more than 10% from avg
795 : !
796 0 : avgdxin = (xin(nlonin)-xin(1))/(nlonin-1._r8)
797 0 : ratio = dxinwrap/avgdxin
798 0 : if (ratio < 0.9_r8 .or. ratio > 1.1_r8) then
799 0 : write(iulog,*)'BILIN: Insane dxinwrap value =',dxinwrap,' avg=', avgdxin
800 0 : call endrun
801 : end if
802 : !
803 : ! Check that interp/extrap points have been found for all outputs, and that
804 : ! they are reasonable (i.e. within 32-bit roundoff)
805 : !
806 0 : icount = 0
807 0 : do j=1,nlatout
808 0 : if (jjm(j) == 0 .or. jjp(j) == 0) icount = icount + 1
809 0 : sum = wgts(j) + wgtn(j)
810 0 : if (sum < 0.99999_r8 .or. sum > 1.00001_r8) icount = icount + 1
811 0 : if (wgts(j) < 0._r8 .or. wgts(j) > 1._r8) icount = icount + 1
812 0 : if (wgtn(j) < 0._r8 .or. wgtn(j) > 1._r8) icount = icount + 1
813 : end do
814 :
815 0 : if (icount > 0) then
816 0 : call endrun ('BILIN: something bad in latitude indices or weights')
817 : end if
818 : !
819 : ! Do the bilinear interpolation
820 : !
821 0 : do j=1,nlatout
822 : !
823 : ! Initialize index arrays for later checking
824 : !
825 0 : do i=1,nlondout
826 0 : iim(i) = 0
827 0 : iip(i) = 0
828 : end do
829 : !
830 : ! For values which extend beyond E and W boundaries, set weights such that
831 : ! values will be interpolated between E and W edges of input grid.
832 : !
833 0 : do iw=1,nlonout(j)
834 0 : if (xout(iw,j) > xin(1)) exit
835 0 : iim(iw) = nlonin
836 0 : iip(iw) = 1
837 0 : wgtw(iw) = (xin(1) - xout(iw,j)) /dxinwrap
838 0 : wgte(iw) = (xout(iw,j)+360._r8 - xin(nlonin))/dxinwrap
839 : end do
840 :
841 0 : do ie=nlonout(j),1,-1
842 0 : if (xout(ie,j) <= xin(nlonin)) exit
843 0 : iim(ie) = nlonin
844 0 : iip(ie) = 1
845 0 : wgtw(ie) = (xin(1)+360._r8 - xout(ie,j)) /dxinwrap
846 0 : wgte(ie) = (xout(ie,j) - xin(nlonin))/dxinwrap
847 : end do
848 : !
849 : ! Loop though output indices finding input indices and weights
850 : !
851 : iiprev = 1
852 0 : do i=iw,ie
853 0 : do ii=iiprev,nlonin-1
854 0 : if (xout(i,j) > xin(ii) .and. xout(i,j) <= xin(ii+1)) then
855 0 : iim(i) = ii
856 0 : iip(i) = ii + 1
857 0 : wgtw(i) = (xin(ii+1) - xout(i,j)) / (xin(ii+1) - xin(ii))
858 0 : wgte(i) = (xout(i,j) - xin(ii)) / (xin(ii+1) - xin(ii))
859 0 : goto 60
860 : end if
861 : end do
862 0 : call endrun ('BILIN: Failed to find interp values')
863 0 : 60 iiprev = ii
864 : end do
865 :
866 0 : icount = 0
867 0 : do i=1,nlonout(j)
868 0 : if (iim(i) == 0 .or. iip(i) == 0) icount = icount + 1
869 0 : sum = wgtw(i) + wgte(i)
870 0 : if (sum < 0.99999_r8 .or. sum > 1.00001_r8) icount = icount + 1
871 0 : if (wgtw(i) < 0._r8 .or. wgtw(i) > 1._r8) icount = icount + 1
872 0 : if (wgte(i) < 0._r8 .or. wgte(i) > 1._r8) icount = icount + 1
873 : end do
874 :
875 0 : if (icount > 0) then
876 0 : write(iulog,*)'BILIN: j=',j,' Something bad in longitude indices or weights'
877 0 : call endrun
878 : end if
879 : !
880 : ! Do the interpolation, 1st in longitude then latitude
881 : !
882 0 : do k=1,nlev
883 0 : do i=1,nlonin
884 0 : igrid(i) = arrin(i,k,jjm(j))*wgts(j) + arrin(i,k,jjp(j))*wgtn(j)
885 : end do
886 :
887 0 : do i=1,nlonout(j)
888 0 : arrout(i,k,j) = igrid(iim(i))*wgtw(i) + igrid(iip(i))*wgte(i)
889 : end do
890 : end do
891 : end do
892 :
893 :
894 0 : return
895 : end subroutine bilin
896 :
897 : !=========================================================================================
898 :
899 5969088 : subroutine vertinterp(ncol, ncold, nlev, pin, pout, arrin, arrout, &
900 0 : extrapolate, ln_interp, ps, phis, tbot)
901 :
902 : ! Vertically interpolate input array to output pressure level. The
903 : ! interpolation is linear in either pressure (the default) or ln pressure.
904 : !
905 : ! If above the top boundary then return top boundary value.
906 : !
907 : ! If below bottom boundary then use bottom boundary value, or optionally
908 : ! for T or Z use the extrapolation algorithm from ECMWF (which was taken
909 : ! from the NCL code base).
910 :
911 : use physconst, only: rair, rga
912 :
913 : !------------------------------Arguments--------------------------------
914 : integer, intent(in) :: ncol ! number active columns
915 : integer, intent(in) :: ncold ! column dimension
916 : integer, intent(in) :: nlev ! vertical dimension
917 : real(r8), intent(in) :: pin(ncold,nlev) ! input pressure levels
918 : real(r8), intent(in) :: pout ! output pressure level
919 : real(r8), intent(in) :: arrin(ncold,nlev) ! input array
920 : real(r8), intent(out) :: arrout(ncold) ! output array (interpolated)
921 :
922 : character(len=1), optional, intent(in) :: extrapolate ! either 'T' or 'Z'
923 : logical, optional, intent(in) :: ln_interp ! set true to interolate
924 : ! in ln(pressure)
925 : real(r8), optional, intent(in) :: ps(ncold) ! surface pressure
926 : real(r8), optional, intent(in) :: phis(ncold) ! surface geopotential
927 : real(r8), optional, intent(in) :: tbot(ncold) ! temperature at bottom level
928 :
929 : !---------------------------Local variables-----------------------------
930 : real(r8) :: alpha
931 : logical :: linear_interp
932 : logical :: do_extrapolate_T, do_extrapolate_Z
933 :
934 : integer :: i,k ! indices
935 5969088 : integer :: kupper(ncold) ! Level indices for interpolation
936 : real(r8) :: dpu ! upper level pressure difference
937 : real(r8) :: dpl ! lower level pressure difference
938 5969088 : logical :: found(ncold) ! true if input levels found
939 : logical :: error ! error flag
940 : !----------------------------------------------------------------------------
941 :
942 2984544 : alpha = 0.0065_r8*rair*rga
943 :
944 2984544 : do_extrapolate_T = .false.
945 2984544 : do_extrapolate_Z = .false.
946 2984544 : if (present(extrapolate)) then
947 :
948 0 : if (extrapolate == 'T') do_extrapolate_T = .true.
949 0 : if (extrapolate == 'Z') do_extrapolate_Z = .true.
950 :
951 0 : if (.not. do_extrapolate_T .and. .not. do_extrapolate_Z) then
952 0 : call endrun ('VERTINTERP: extrapolate must be T or Z')
953 : end if
954 0 : if (.not. present(ps)) then
955 0 : call endrun ('VERTINTERP: ps required for extrapolation')
956 : end if
957 0 : if (.not. present(phis)) then
958 0 : call endrun ('VERTINTERP: phis required for extrapolation')
959 : end if
960 0 : if (do_extrapolate_Z) then
961 0 : if (.not. present(tbot)) then
962 0 : call endrun ('VERTINTERP: extrapolate must be T or Z')
963 : end if
964 : end if
965 : end if
966 :
967 2984544 : linear_interp = .true.
968 2984544 : if (present(ln_interp)) then
969 0 : if (ln_interp) linear_interp = .false.
970 : end if
971 :
972 49834944 : do i = 1, ncol
973 46850400 : found(i) = .false.
974 49834944 : kupper(i) = 1
975 : end do
976 2984544 : error = .false.
977 :
978 : ! Find indices of upper pressure bound
979 280547136 : do k = 1, nlev - 1
980 4637634336 : do i = 1, ncol
981 4634649792 : if ((.not. found(i)) .and. pin(i,k)<pout .and. pout<=pin(i,k+1)) then
982 46850400 : found(i) = .true.
983 46850400 : kupper(i) = k
984 : end if
985 : end do
986 : end do
987 :
988 49834944 : do i = 1, ncol
989 :
990 49834944 : if (pout <= pin(i,1)) then
991 :
992 : ! if above top pressure level then use value at top (no interpolation)
993 0 : arrout(i) = arrin(i,1)
994 :
995 46850400 : else if (pout >= pin(i,nlev)) then
996 :
997 0 : if (do_extrapolate_T) then
998 : ! use ECMWF algorithm to extrapolate T
999 0 : arrout(i) = extrapolate_T()
1000 :
1001 0 : else if (do_extrapolate_Z) then
1002 : ! use ECMWF algorithm to extrapolate Z
1003 0 : arrout(i) = extrapolate_Z()
1004 :
1005 : else
1006 : ! use bottom value
1007 0 : arrout(i) = arrin(i,nlev)
1008 : end if
1009 :
1010 46850400 : else if (found(i)) then
1011 46850400 : if (linear_interp) then
1012 : ! linear interpolation in p
1013 46850400 : dpu = pout - pin(i,kupper(i))
1014 46850400 : dpl = pin(i,kupper(i)+1) - pout
1015 46850400 : arrout(i) = (arrin(i,kupper(i))*dpl + arrin(i,kupper(i)+1)*dpu)/(dpl + dpu)
1016 : else
1017 : ! linear interpolation in ln(p)
1018 0 : arrout(i) = arrin(i,kupper(i)) + (arrin(i,kupper(i)+1) - arrin(i,kupper(i)))* &
1019 0 : log(pout/pin(i,kupper(i))) / &
1020 0 : log(pin(i,kupper(i)+1)/pin(i,kupper(i)))
1021 : end if
1022 : else
1023 : error = .true.
1024 : end if
1025 : end do
1026 :
1027 : ! Error check
1028 5969088 : if (error) then
1029 0 : call endrun ('VERTINTERP: ERROR FLAG')
1030 : end if
1031 :
1032 : contains
1033 :
1034 : !======================================================================================
1035 :
1036 0 : real(r8) function extrapolate_T()
1037 :
1038 : ! local variables
1039 : real(r8) :: sixth
1040 : real(r8) :: tstar
1041 : real(r8) :: hgt
1042 : real(r8) :: alnp
1043 : real(r8) :: t0
1044 : real(r8) :: tplat
1045 : real(r8) :: tprime0
1046 : !-------------------------------------------------------------------------
1047 :
1048 0 : sixth = 1._r8 / 6._r8
1049 0 : tstar = arrin(i,nlev)*(1._r8 + alpha*(ps(i)/pin(i,nlev) - 1._r8))
1050 0 : hgt = phis(i)*rga
1051 :
1052 0 : if (hgt < 2000._r8) then
1053 0 : alnp = alpha*log(pout/ps(i))
1054 : else
1055 0 : t0 = tstar + 0.0065_r8*hgt
1056 0 : tplat = min(t0, 298._r8)
1057 :
1058 0 : if (hgt <= 2500._r8) then
1059 : tprime0 = 0.002_r8*((2500._r8 - hgt)*t0 + &
1060 0 : (hgt - 2000._r8)*tplat)
1061 : else
1062 : tprime0 = tplat
1063 : end if
1064 :
1065 0 : if (tprime0 < tstar) then
1066 : alnp = 0._r8
1067 : else
1068 0 : alnp = rair*(tprime0 - tstar)/phis(i) * log(pout/ps(i))
1069 : end if
1070 :
1071 : end if
1072 :
1073 0 : extrapolate_T = tstar*(1._r8 + alnp + 0.5_r8*alnp**2 + sixth*alnp**3)
1074 :
1075 2984544 : end function extrapolate_T
1076 :
1077 : !======================================================================================
1078 :
1079 0 : real(r8) function extrapolate_Z()
1080 :
1081 : ! local variables
1082 : real(r8) :: sixth
1083 : real(r8) :: hgt
1084 : real(r8) :: tstar
1085 : real(r8) :: t0
1086 : real(r8) :: alph
1087 : real(r8) :: alnp
1088 : !-------------------------------------------------------------------------
1089 :
1090 0 : sixth = 1._r8 / 6._r8
1091 0 : hgt = phis(i)*rga
1092 0 : tstar = tbot(i)*(1._r8 + alpha*(ps(i)/pin(i,nlev) - 1._r8))
1093 0 : t0 = tstar + 0.0065_r8*hgt
1094 :
1095 0 : if (tstar <= 290.5_r8 .and. t0 > 290.5_r8) then
1096 0 : alph = rair/phis(i) * (290.5_r8 - tstar)
1097 :
1098 0 : else if (tstar > 290.5_r8 .and. t0 > 290.5_r8) then
1099 0 : alph = 0._r8
1100 0 : tstar = 0.5_r8*(290.5_r8 + tstar)
1101 :
1102 : else
1103 : alph = alpha
1104 :
1105 : end if
1106 :
1107 0 : if (tstar < 255._r8) then
1108 0 : tstar = 0.5_r8*(tstar + 255._r8)
1109 : end if
1110 0 : alnp = alph*log(pout/ps(i))
1111 :
1112 : extrapolate_Z = hgt - rair*tstar*rga*log(pout/ps(i))* &
1113 0 : (1._r8 + 0.5_r8*alnp + sixth*alnp**2)
1114 :
1115 0 : end function extrapolate_Z
1116 :
1117 : end subroutine vertinterp
1118 :
1119 : !=========================================================================================
1120 :
1121 0 : subroutine get_timeinterp_factors (cycflag, np1, cdayminus, cdayplus, cday, &
1122 : fact1, fact2, str)
1123 : !---------------------------------------------------------------------------
1124 : !
1125 : ! Purpose: Determine time interpolation factors (normally for a boundary dataset)
1126 : ! for linear interpolation.
1127 : !
1128 : ! Method: Assume 365 days per year. Output variable fact1 will be the weight to
1129 : ! apply to data at calendar time "cdayminus", and fact2 the weight to apply
1130 : ! to data at time "cdayplus". Combining these values will produce a result
1131 : ! valid at time "cday". Output arguments fact1 and fact2 will be between
1132 : ! 0 and 1, and fact1 + fact2 = 1 to roundoff.
1133 : !
1134 : ! Author: Jim Rosinski
1135 : !
1136 : !---------------------------------------------------------------------------
1137 : implicit none
1138 : !
1139 : ! Arguments
1140 : !
1141 : logical, intent(in) :: cycflag ! flag indicates whether dataset is being cycled yearly
1142 :
1143 : integer, intent(in) :: np1 ! index points to forward time slice matching cdayplus
1144 :
1145 : real(r8), intent(in) :: cdayminus ! calendar day of rearward time slice
1146 : real(r8), intent(in) :: cdayplus ! calendar day of forward time slice
1147 : real(r8), intent(in) :: cday ! calenar day to be interpolated to
1148 : real(r8), intent(out) :: fact1 ! time interpolation factor to apply to rearward time slice
1149 : real(r8), intent(out) :: fact2 ! time interpolation factor to apply to forward time slice
1150 :
1151 : character(len=*), intent(in) :: str ! string to be added to print in case of error (normally the callers name)
1152 : !
1153 : ! Local workspace
1154 : !
1155 : real(r8) :: deltat ! time difference (days) between cdayminus and cdayplus
1156 : real(r8), parameter :: daysperyear = 365._r8 ! number of days in a year
1157 : !
1158 : ! Initial sanity checks
1159 : !
1160 0 : if (np1 == 1 .and. .not. cycflag) then
1161 0 : call endrun ('GETFACTORS:'//str//' cycflag false and forward month index = Jan. not allowed')
1162 : end if
1163 :
1164 0 : if (np1 < 1) then
1165 0 : call endrun ('GETFACTORS:'//str//' input arg np1 must be > 0')
1166 : end if
1167 :
1168 0 : if (cycflag) then
1169 0 : if ((cday < 1._r8) .or. (cday > (daysperyear+1._r8))) then
1170 0 : write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday
1171 0 : call endrun ()
1172 : end if
1173 : else
1174 0 : if (cday < 1._r8) then
1175 0 : write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday
1176 0 : call endrun ()
1177 : end if
1178 : end if
1179 : !
1180 : ! Determine time interpolation factors. Account for December-January
1181 : ! interpolation if dataset is being cycled yearly.
1182 : !
1183 0 : if (cycflag .and. np1 == 1) then ! Dec-Jan interpolation
1184 0 : deltat = cdayplus + daysperyear - cdayminus
1185 0 : if (cday > cdayplus) then ! We are in December
1186 0 : fact1 = (cdayplus + daysperyear - cday)/deltat
1187 0 : fact2 = (cday - cdayminus)/deltat
1188 : else ! We are in January
1189 0 : fact1 = (cdayplus - cday)/deltat
1190 0 : fact2 = (cday + daysperyear - cdayminus)/deltat
1191 : end if
1192 : else
1193 0 : deltat = cdayplus - cdayminus
1194 0 : fact1 = (cdayplus - cday)/deltat
1195 0 : fact2 = (cday - cdayminus)/deltat
1196 : end if
1197 :
1198 0 : if (.not. valid_timeinterp_factors (fact1, fact2)) then
1199 0 : write(iulog,*) 'GETFACTORS: ', str, ' bad fact1 and/or fact2=', fact1, fact2
1200 0 : call endrun ()
1201 : end if
1202 :
1203 0 : return
1204 : end subroutine get_timeinterp_factors
1205 :
1206 0 : logical function valid_timeinterp_factors (fact1, fact2)
1207 : !---------------------------------------------------------------------------
1208 : !
1209 : ! Purpose: check sanity of time interpolation factors to within 32-bit roundoff
1210 : !
1211 : !---------------------------------------------------------------------------
1212 : implicit none
1213 :
1214 : real(r8), intent(in) :: fact1, fact2 ! time interpolation factors
1215 :
1216 0 : valid_timeinterp_factors = .true.
1217 :
1218 : ! The fact1 .ne. fact1 and fact2 .ne. fact2 comparisons are to detect NaNs.
1219 : if (abs(fact1+fact2-1._r8) > 1.e-6_r8 .or. &
1220 : fact1 > 1.000001_r8 .or. fact1 < -1.e-6_r8 .or. &
1221 : fact2 > 1.000001_r8 .or. fact2 < -1.e-6_r8 .or. &
1222 0 : fact1 .ne. fact1 .or. fact2 .ne. fact2) then
1223 :
1224 0 : valid_timeinterp_factors = .false.
1225 : end if
1226 :
1227 : return
1228 : end function valid_timeinterp_factors
1229 :
1230 0 : end module interpolate_data
|
