Line data Source code
1 : !-------------------------------------------------------------------------------
2 : !$Id$
3 : !===============================================================================
4 : module interpolation
5 :
6 : use clubb_precision, only: &
7 : core_rknd ! Variable(s)
8 :
9 : implicit none
10 :
11 : private ! Default Scope
12 :
13 : public :: lin_interpolate_two_points, binary_search, zlinterp_fnc, &
14 : lin_interpolate_on_grid, linear_interp_factor, mono_cubic_interp, plinterp_fnc, &
15 : pvertinterp
16 :
17 : contains
18 :
19 : !-------------------------------------------------------------------------------
20 0 : function lin_interpolate_two_points( height_int, height_high, height_low, &
21 : var_high, var_low )
22 :
23 : ! Description:
24 : ! This function computes a linear interpolation of the value of variable.
25 : ! Given two known values of a variable at two height values, the value
26 : ! of that variable at a height between those two height levels (rather
27 : ! than a height outside of those two height levels) is computed.
28 : !
29 : ! Here is a diagram:
30 : !
31 : ! ################################ Height high, know variable value
32 : !
33 : !
34 : !
35 : ! -------------------------------- Height to be interpolated to; linear interpolation
36 : !
37 : !
38 : !
39 : !
40 : !
41 : ! ################################ Height low, know variable value
42 : !
43 : !
44 : ! FORMULA:
45 : !
46 : ! variable(@ Height interpolation) =
47 : !
48 : ! [ (variable(@ Height high) - variable(@ Height low)) / (Height high - Height low) ]
49 : ! * (Height interpolation - Height low) + variable(@ Height low)
50 :
51 : ! Comments from WRF-HOC, Brian Griffin.
52 :
53 : ! References:
54 : ! None
55 : !-------------------------------------------------------------------------------
56 :
57 : use clubb_precision, only: &
58 : core_rknd ! Variable(s)
59 :
60 : use constants_clubb, only: fstderr ! Constant
61 :
62 : implicit none
63 :
64 : ! Input Variables
65 :
66 : real( kind = core_rknd ), intent(in) :: &
67 : height_int, & ! Height to be interpolated to [m]
68 : height_high, & ! Height above the interpolation [m]
69 : height_low, & ! Height below the interpolation [m]
70 : var_high, & ! Variable above the interpolation [units vary]
71 : var_low ! Variable below the interpolation [units vary]
72 :
73 : ! Output Variables
74 : real( kind = core_rknd ) :: lin_interpolate_two_points
75 :
76 : ! Check for valid input
77 0 : if ( abs(height_low - height_high) < 1.0e-12_core_rknd ) then
78 0 : write(fstderr,*) "lin_interpolate_two_points: height_high and height_low cannot be equal."
79 0 : error stop
80 : end if
81 :
82 : ! Compute linear interpolation
83 :
84 : lin_interpolate_two_points = ( ( height_int - height_low )/( height_high - height_low ) ) &
85 0 : * ( var_high - var_low ) + var_low
86 :
87 : return
88 : end function lin_interpolate_two_points
89 :
90 : !-------------------------------------------------------------------------------------------------
91 0 : elemental real( kind = core_rknd ) function linear_interp_factor( factor, var_high, var_low )
92 : ! Description:
93 : ! Determines the coefficient for a linear interpolation
94 : !
95 : ! References:
96 : ! None
97 : !-------------------------------------------------------------------------------------------------
98 :
99 : use clubb_precision, only: &
100 : core_rknd ! Variable(s)
101 :
102 : implicit none
103 :
104 : real( kind = core_rknd ), intent(in) :: &
105 : factor, & ! Factor [units vary]
106 : var_high, & ! Variable above the interpolation [units vary]
107 : var_low ! Variable below the interpolation [units vary]
108 :
109 0 : linear_interp_factor = factor * ( var_high - var_low ) + var_low
110 :
111 : return
112 : end function linear_interp_factor
113 : !-------------------------------------------------------------------------------------------------
114 0 : function mono_cubic_interp &
115 : ( z_in, km1, k00, kp1, kp2, zm1, z00, zp1, zp2, fm1, f00, fp1, fp2 ) result ( f_out )
116 :
117 : ! Description:
118 : ! Steffen's monotone cubic interpolation method
119 : ! Returns monotone cubic interpolated value between x00 and xp1
120 :
121 : ! Original Author:
122 : ! Takanobu Yamaguchi
123 : ! tak.yamaguchi@noaa.gov
124 : !
125 : ! This version has been modified slightly for CLUBB's coding standards and
126 : ! adds the 3/2 from eqn 21. -dschanen 26 Oct 2011
127 : ! We have also added a quintic polynomial option.
128 : !
129 : ! References:
130 : ! M. Steffen, Astron. Astrophys. 239, 443-450 (1990)
131 : !-------------------------------------------------------------------------------------------------
132 :
133 : use constants_clubb, only: &
134 : eps ! Constant
135 :
136 : use clubb_precision, only: &
137 : core_rknd ! Constant
138 :
139 : use model_flags, only: &
140 : l_quintic_poly_interp ! Variable(s)
141 :
142 : implicit none
143 :
144 : ! External
145 : intrinsic :: sign, abs, min
146 :
147 : ! Input Variables
148 : real( kind = core_rknd ), intent(in) :: &
149 : z_in ! The altitude to be interpolated to [m]
150 :
151 : ! k-levels; their meaning depends on whether we're extrapolating or interpolating
152 : integer, intent(in) :: &
153 : km1, k00, kp1, kp2
154 :
155 : real( kind = core_rknd ), intent(in) :: &
156 : zm1, z00, zp1, zp2, & ! The altitudes for km1, k00, kp1, kp2 [m]
157 : fm1, f00, fp1, fp2 ! The field at km1, k00, kp1, and kp2 [units vary]
158 :
159 : ! Output Variables
160 : real( kind = core_rknd ) :: f_out ! The interpolated field
161 :
162 : ! Local Variables
163 : real( kind = core_rknd ) :: &
164 : hm1, h00, hp1, &
165 : sm1, s00, sp1, &
166 : p00, pp1, &
167 : dfdx00, dfdxp1, &
168 : c1, c2, c3, c4, &
169 : w00, wp1, &
170 : coef1, coef2, &
171 : zprime, beta, alpha, zn
172 :
173 : ! ---- Begin Code ----
174 :
175 0 : coef1 = 1.0_core_rknd
176 0 : coef2 = 1.0_core_rknd
177 :
178 0 : if ( km1 <= k00 ) then
179 0 : hm1 = z00 - zm1
180 0 : h00 = zp1 - z00
181 0 : hp1 = zp2 - zp1
182 :
183 0 : if ( km1 == k00 ) then
184 0 : s00 = ( fp1 - f00 ) / ( zp1 - z00 )
185 0 : sp1 = ( fp2 - fp1 ) / ( zp2 - zp1 )
186 0 : dfdx00 = s00
187 0 : pp1 = ( s00 * hp1 + sp1 * h00 ) / ( h00 + hp1 )
188 : dfdxp1 = coef1*( sign( 1.0_core_rknd, s00 ) + sign( 1.0_core_rknd, sp1 ) ) &
189 0 : * min( abs( s00 ), abs( sp1 ), coef2*0.5_core_rknd*abs( pp1 ) )
190 :
191 0 : else if ( kp1 == kp2 ) then
192 0 : sm1 = ( f00 - fm1 ) / ( z00 - zm1 )
193 0 : s00 = ( fp1 - f00 ) / ( zp1 - z00 )
194 0 : p00 = ( sm1 * h00 + s00 * hm1 ) / ( hm1 + h00 )
195 : dfdx00 = coef1*( sign( 1.0_core_rknd, sm1 ) + sign( 1.0_core_rknd, s00 ) ) &
196 0 : * min( abs( sm1 ), abs( s00 ), coef2*0.5_core_rknd*abs( p00 ) )
197 0 : dfdxp1 = s00
198 :
199 : else
200 0 : sm1 = ( f00 - fm1 ) / ( z00 - zm1 )
201 0 : s00 = ( fp1 - f00 ) / ( zp1 - z00 )
202 0 : sp1 = ( fp2 - fp1 ) / ( zp2 - zp1 )
203 0 : p00 = ( sm1 * h00 + s00 * hm1 ) / ( hm1 + h00 )
204 0 : pp1 = ( s00 * hp1 + sp1 * h00 ) / ( h00 + hp1 )
205 : dfdx00 = coef1*( sign( 1.0_core_rknd, sm1 ) + sign( 1.0_core_rknd, s00 ) ) &
206 0 : * min( abs( sm1 ), abs( s00 ), coef2*0.5_core_rknd*abs( p00 ) )
207 : dfdxp1 = coef1*( sign( 1.0_core_rknd, s00 ) + sign( 1.0_core_rknd, sp1 ) ) &
208 0 : * min( abs( s00 ), abs( sp1 ), coef2*0.5_core_rknd*abs( pp1 ) )
209 :
210 : end if
211 :
212 0 : c1 = ( dfdx00 + dfdxp1 - 2._core_rknd * s00 ) / ( h00 ** 2 )
213 0 : c2 = ( 3._core_rknd * s00 - 2._core_rknd * dfdx00 - dfdxp1 ) / h00
214 0 : c3 = dfdx00
215 0 : c4 = f00
216 :
217 : if ( .not. l_quintic_poly_interp ) then
218 :
219 : ! Old formula
220 : !f_out = c1 * ( (z_in - z00)**3 ) + c2 * ( (z_in - z00)**2 ) + c3 * (z_in - z00) + c4
221 :
222 : ! Faster nested multiplication
223 0 : zprime = z_in - z00
224 0 : f_out = c4 + zprime*( c3 + zprime*( c2 + ( zprime*c1 ) ) )
225 :
226 : else
227 :
228 : ! Use a quintic polynomial interpolation instead instead of the Steffen formula.
229 : ! Unlike the formula above, this formula does not guarantee monotonicity.
230 :
231 : beta = 120._core_rknd * ( (fp1-f00) - 0.5_core_rknd * h00 * (dfdx00 + dfdxp1) )
232 :
233 : ! Prevent an underflow by using a linear interpolation
234 : if ( abs( beta ) < eps ) then
235 : f_out = lin_interpolate_two_points( z00, zp1, zm1, &
236 : fp1, fm1 )
237 :
238 : else
239 : alpha = (6._core_rknd/beta) * h00 * (dfdxp1-dfdx00) + 0.5_core_rknd
240 : zn = (z_in-z00)/h00
241 :
242 : f_out = ( &
243 : (( (beta/20._core_rknd)*zn - (beta*(1._core_rknd+alpha) &
244 : / 12._core_rknd)) * zn + (beta*alpha/6._core_rknd)) &
245 : * zn**2 + dfdx00*h00 &
246 : ) * zn + f00
247 : end if ! beta < eps
248 : end if ! ~quintic_polynomial
249 :
250 : else
251 : ! Linear extrapolation
252 0 : wp1 = ( z_in - z00 ) / ( zp1 - z00 )
253 0 : w00 = 1._core_rknd - wp1
254 0 : f_out = wp1 * fp1 + w00 * f00
255 :
256 : end if
257 :
258 : return
259 : end function mono_cubic_interp
260 :
261 : !-------------------------------------------------------------------------------
262 0 : integer function binary_search( n, array, var ) &
263 : result( i )
264 :
265 : ! Description:
266 : ! This subroutine performs a binary search to find the closest value greater
267 : ! than or equal to var in the array. This function returns the index of the
268 : ! closest value of array that is greater than or equal to var. It returns a
269 : ! value of -1 if var is outside the bounds of array.
270 : !
271 : !-----------------------------------------------------------------------
272 :
273 : use clubb_precision, only: &
274 : core_rknd ! Variable(s)
275 :
276 : use constants_clubb, only: &
277 : fstderr ! Variable(s)
278 :
279 : use error_code, only: &
280 : clubb_at_least_debug_level ! Error indicator
281 :
282 : implicit none
283 :
284 : ! Input Variables
285 :
286 : ! Size of the array
287 : integer, intent(in) :: n
288 :
289 : ! The array being searched (must be sorted from least value to greatest
290 : ! value).
291 : real( kind = core_rknd ), dimension(n), intent(in) :: array
292 :
293 : ! The value being searched for
294 : real( kind = core_rknd ), intent(in) :: var
295 :
296 : ! Bounds of the search
297 : integer :: high
298 : integer :: low
299 :
300 : ! The initial value of low has been changed from 1 to 2 due to a problem
301 : ! that was occuring when var was close to the lower bound.
302 : !
303 : ! The lowest value in the array (which is sorted by increasing values) is
304 : ! found at index 1, while the highest value in the array is found at
305 : ! index n. Unless the value of var exactly corresponds with one of the
306 : ! values found in the array, or unless the value of var is found outside of
307 : ! the array, the value of var will be found between two levels of the array.
308 : ! In this scenario, the output of function binary_search is the index of the
309 : ! HIGHER level. For example, if the value of var is found between array(1)
310 : ! and array(2), the output of function binary_search will be 2.
311 : !
312 : ! Therefore, the lowest index of a HIGHER level in an interpolation is 2.
313 : ! Thus, the initial value of low has been changed to 2. This will prevent
314 : ! the value of variable "i" below from becoming 1. If the value of "i"
315 : ! becomes 1, the code below tries to access array(0) (which is array(i-1)
316 : ! when i = 1) and produces an error.
317 :
318 0 : low = 2
319 :
320 0 : high = n
321 :
322 : ! Ensure variable is within range of array and array has more than 1 element
323 0 : if ( var < array(1) .or. var > array(n) .or. n < 2 ) then
324 0 : i = -1
325 : return
326 : end if
327 :
328 : ! Special case to cover var == array(1)
329 0 : if ( var >= array(1) .and. var <= array(2) ) then
330 0 : i = 2
331 : return
332 : end if
333 :
334 :
335 0 : do while( low <= high )
336 :
337 0 : i = (low + high) / 2
338 :
339 0 : if ( var > array(i-1) .and. var <= array(i) ) then
340 :
341 : return
342 :
343 0 : elseif ( var < array(i) ) then
344 :
345 : high = i - 1
346 :
347 0 : elseif ( var > array(i) ) then
348 :
349 0 : low = i + 1
350 :
351 : endif
352 :
353 : enddo
354 :
355 : ! Code should not get to this point, but return -1 to be safe
356 0 : if ( clubb_at_least_debug_level( 0 ) ) then
357 0 : write(fstderr,*) "Logic error in function binary_search."
358 : end if
359 :
360 0 : i = -1
361 : return
362 :
363 : end function binary_search
364 :
365 : !-------------------------------------------------------------------------------
366 0 : function plinterp_fnc( dim_out, dim_src, grid_out, &
367 0 : grid_src, var_src ) &
368 0 : result( var_out )
369 : ! Description:
370 : ! Do a linear interpolation in the vertical with pressures. Assumes
371 : ! values that are less than lowest source point are zero and above the
372 : ! highest source point are zero. Also assumes altitude increases linearly.
373 : ! This function just calls zlinterp_fnc, but negates grid_out and grid_src.
374 :
375 : ! References:
376 : ! function LIN_INT from WRF-HOC
377 : !-----------------------------------------------------------------------
378 :
379 : use clubb_precision, only: &
380 : core_rknd ! Variable(s)
381 :
382 : implicit none
383 :
384 : ! Input variables
385 : integer, intent(in) :: dim_out, dim_src
386 :
387 : real( kind = core_rknd ), dimension(dim_src), intent(in) :: &
388 : grid_src, & ! [m]
389 : var_src ! [units vary]
390 :
391 : real( kind = core_rknd ), dimension(dim_out), intent(in) :: &
392 : grid_out ! [m]
393 :
394 : ! Output variable
395 : real( kind = core_rknd ), dimension(dim_out) :: &
396 : var_out ! [units vary]
397 :
398 : ! ---- Begin Code ----
399 :
400 : var_out = zlinterp_fnc( dim_out, dim_src, -grid_out, &
401 0 : -grid_src, var_src )
402 :
403 0 : return
404 0 : end function plinterp_fnc
405 : !-------------------------------------------------------------------------------
406 0 : function zlinterp_fnc( dim_out, dim_src, grid_out, &
407 0 : grid_src, var_src ) &
408 0 : result( var_out )
409 : ! Description:
410 : ! Do a linear interpolation in the vertical. Assumes values that
411 : ! are less than lowest source point are zero and above the highest
412 : ! source point are zero. Also assumes altitude increases linearly.
413 :
414 : ! References:
415 : ! function LIN_INT from WRF-HOC
416 : !-----------------------------------------------------------------------
417 :
418 : use clubb_precision, only: &
419 : core_rknd ! Variable(s)
420 :
421 : implicit none
422 :
423 : ! Input variables
424 : integer, intent(in) :: dim_out, dim_src
425 :
426 : real( kind = core_rknd ), dimension(dim_src), intent(in) :: &
427 : grid_src, & ! [m]
428 : var_src ! [units vary]
429 :
430 : real( kind = core_rknd ), dimension(dim_out), intent(in) :: &
431 : grid_out ! [m]
432 :
433 : ! Output variable
434 : real( kind = core_rknd ), dimension(dim_out) :: &
435 : var_out ! [units vary]
436 :
437 : ! Local variables
438 : integer :: k, kint, km1
439 :
440 : ! integer :: tst, kp1
441 :
442 : ! ---- Begin Code ----
443 :
444 0 : k = 1
445 :
446 0 : do kint = 1, dim_out, 1
447 :
448 : ! Set to 0 if we're below the input data's lowest point
449 0 : if ( grid_out(kint) < grid_src(1) ) then
450 0 : var_out(kint) = 0.0_core_rknd
451 0 : cycle
452 : end if
453 :
454 : ! Increment k until the level is correct
455 : ! do while ( grid_out(kint) > grid_src(k)
456 : ! . .and. k < dim_src )
457 : ! k = k + 1
458 : ! end do
459 :
460 : ! Changed so a binary search is used instead of a sequential search
461 : ! tst = binary_search(dim_src, grid_src, grid_out(kint))
462 0 : k = binary_search(dim_src, grid_src, grid_out(kint))
463 : ! Joshua Fasching April 2008
464 :
465 : ! print *, "k = ", k
466 : ! print *, "tst = ", tst
467 : ! print *, "dim_src = ", dim_src
468 : ! print *,"------------------------------"
469 :
470 : ! If the increment leads to a level above the data, set this
471 : ! point and all those above it to zero
472 : !if( k > dim_src ) then
473 0 : if ( k == -1 ) then
474 0 : var_out(kint:dim_out) = 0.0_core_rknd
475 : exit
476 : end if
477 :
478 0 : km1 = max( 1, k-1 )
479 : !kp1 = min( k+1, dim_src )
480 :
481 : ! Interpolate
482 0 : var_out(kint) = lin_interpolate_two_points( grid_out(kint), grid_src(k), &
483 0 : grid_src(km1), var_src(k), var_src(km1) )
484 :
485 : ! ( var_src(k) - var_src(km1) ) / &
486 : ! ( grid_src(k) - grid_src(km1) ) &
487 : ! * ( grid_out(kint) - grid_src(km1) ) + var_src(km1) &
488 : ! Changed to use a standard function for interpolation
489 :
490 : !! Note this ends up changing the results slightly because
491 : !the placement of variables has been changed.
492 :
493 : ! Joshua Fasching April 2008
494 :
495 : end do ! kint = 1..dim_out
496 :
497 0 : return
498 0 : end function zlinterp_fnc
499 :
500 : !-------------------------------------------------------------------------------
501 17870112 : subroutine pvertinterp( nz, ngrdcol, &
502 17870112 : p_mid, p_out, input_var, &
503 17870112 : interp_var )
504 :
505 : implicit none
506 :
507 : !------------------------ Input Variables ------------------------
508 : integer , intent(in) :: &
509 : nz, &
510 : ngrdcol
511 :
512 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
513 : p_mid, & ! input level pressure levels
514 : input_var ! input array
515 :
516 : real( kind = core_rknd ), intent(in) :: &
517 : p_out ! output pressure level
518 :
519 : !------------------------ Output Variables ------------------------
520 : real( kind = core_rknd ), dimension(ngrdcol), intent(out) :: &
521 : interp_var ! output array (interpolated)
522 :
523 : !------------------------ Local Variables ------------------------
524 : integer :: &
525 : i, k, & ! Loop indices
526 : k_upper ! Level indices for interpolation
527 :
528 : real( kind = core_rknd ) :: &
529 : dpu, & ! upper level pressure difference
530 : dpl ! lower level pressure difference
531 :
532 : logical :: &
533 : l_found, & ! true if input levels found
534 : l_error ! true if error
535 :
536 : !------------------------ Begin Code ------------------------
537 :
538 : ! Initialize index array and logical flags
539 17870112 : l_error = .false.
540 :
541 : ! If we've fallen through the k=1,nz-1 loop, we cannot interpolate and
542 : ! must extrapolate from the bottom or top data level for at least some
543 : ! of the longitude points.
544 : !$acc parallel loop gang vector default(present)
545 298389312 : do i = 1, ngrdcol
546 :
547 298389312 : if ( p_out >= p_mid(i,1) ) then
548 :
549 54547740 : interp_var(i) = input_var(i,1)
550 :
551 225971460 : elseif ( p_out <= p_mid(i,nz) ) then
552 :
553 0 : interp_var(i) = input_var(i,nz)
554 :
555 : else
556 :
557 : l_found = .false.
558 2966159100 : k_upper = 1
559 :
560 : ! Store level indices for interpolation.
561 : ! If all indices for this level have been found,
562 : ! do the interpolation
563 2966159100 : do k = 1, nz-1
564 2966159100 : if ( p_mid(i,k) > p_out .and. p_out >= p_mid(i,k+1) ) then
565 : l_found = .true.
566 : k_upper = k
567 : exit
568 : end if
569 : end do
570 :
571 : if ( .not. l_found ) then
572 225971460 : l_error = .true.
573 : end if
574 :
575 225971460 : dpu = p_mid(i,k_upper) - p_out
576 225971460 : dpl = p_out - p_mid(i,k_upper+1)
577 : interp_var(i) = ( input_var(i,k_upper)*dpl + input_var(i,k_upper+1)*dpu ) &
578 225971460 : / ( dpl + dpu )
579 : end if
580 :
581 : end do
582 : !$acc end parallel loop
583 :
584 17870112 : return
585 :
586 : end subroutine pvertinterp
587 :
588 : !-------------------------------------------------------------------------------
589 0 : subroutine lin_interpolate_on_grid &
590 0 : ( nparam, xlist, tlist, xvalue, &
591 : tvalue )
592 :
593 : ! Description:
594 : ! Linear interpolation for 25 June 1996 altocumulus case.
595 :
596 : ! For example, to interpolate between two temperatures in space, put
597 : ! your spatial coordinates in x-list and your temperature values in
598 : ! tlist. The point in question should have its spatial value stored
599 : ! in xvalue, and tvalue will be the temperature at that point.
600 :
601 : ! Author: Michael Falk for COAMPS.
602 : !-------------------------------------------------------------------------------
603 :
604 : use constants_clubb, only: fstderr ! Constant
605 :
606 : use clubb_precision, only: &
607 : core_rknd ! Variable(s)
608 :
609 : use error_code, only: &
610 : clubb_at_least_debug_level ! Error indicator
611 :
612 : implicit none
613 :
614 : ! Input Variables
615 : integer, intent(in) :: nparam ! Number of parameters in xlist and tlist
616 :
617 : ! Input/Output Variables
618 : real( kind = core_rknd ), intent(inout), dimension(nparam) :: &
619 : xlist, & ! List of x-values (independent variable)
620 : tlist ! List of t-values (dependent variable)
621 :
622 : real( kind = core_rknd ), intent(in) :: &
623 : xvalue ! x-value at which to interpolate
624 :
625 : real( kind = core_rknd ), intent(inout) :: &
626 : tvalue ! t-value solved by interpolation
627 :
628 : ! Local variables
629 : integer :: &
630 : i ! Loop control variable
631 :
632 : integer :: &
633 : bottombound, & ! Index of the smaller value in the linear interpolation
634 : topbound ! Index of the larger value in the linear interpolation
635 :
636 : !-------------------------------------------------------------------------------
637 : !
638 : ! Assure that the elements are in order so that the interpolation is between
639 : ! the two closest points to the point in question.
640 : !
641 : !-------------------------------------------------------------------------------
642 :
643 0 : if ( clubb_at_least_debug_level( 0 ) ) then
644 0 : do i=2,nparam
645 0 : if ( xlist(i) <= xlist(i-1) ) then
646 0 : write(fstderr,*) "xlist must be sorted for lin_interpolate_on_grid."
647 0 : error stop
648 : end if
649 : end do
650 : end if
651 : !-------------------------------------------------------------------------------
652 : !
653 : ! If the point in question is larger than the largest x-value or
654 : ! smaller than the smallest x-value, crash.
655 : !
656 : !-------------------------------------------------------------------------------
657 :
658 0 : if ( clubb_at_least_debug_level( 0 ) ) then
659 0 : if ( (xvalue < xlist(1)) .or. (xvalue > xlist(nparam)) ) then
660 0 : write(fstderr,*) "lin_interpolate_on_grid: Value out of range"
661 0 : error stop
662 : end if
663 : end if
664 :
665 : !-------------------------------------------------------------------------------
666 : !
667 : ! Find the correct top and bottom bounds, do the interpolation, return c
668 : ! the value.
669 : !
670 : !-------------------------------------------------------------------------------
671 :
672 : topbound = -1
673 : bottombound = -1
674 :
675 0 : do i=2,nparam
676 0 : if ( (xvalue >= xlist(i-1)) .and. (xvalue <= xlist(i)) ) then
677 0 : bottombound = i-1
678 0 : topbound = i
679 : end if
680 : end do
681 :
682 0 : if ( clubb_at_least_debug_level( 1 ) .and. (topbound == -1 .or. bottombound == -1) ) then
683 0 : write(fstderr,*) "Sanity check failed! xlist is not properly sorted"
684 0 : write(fstderr,*) "in lin_interpolate_on_grid."
685 : end if
686 :
687 : tvalue = &
688 0 : lin_interpolate_two_points( xvalue, xlist(topbound), xlist(bottombound), &
689 0 : tlist(topbound), tlist(bottombound) )
690 :
691 0 : return
692 : end subroutine lin_interpolate_on_grid
693 :
694 : end module interpolation
|
