mo_inter.F90 coverage: 0.00 %func 0.00 %block
1)
2) module mo_inter
3)
4) use shr_kind_mod, only : r8 => shr_kind_r8
5) use cam_logfile, only : iulog
6)
7) implicit none
8)
9) integer :: nintervals
10) integer, allocatable :: xi(:), xcnt(:)
11) real(r8), allocatable :: xfrac(:,:)
12)
13) contains
14)
15) subroutine inter2( ng, xg, yg, n, x, y, ierr )
16) !-----------------------------------------------------------------------------
17) ! purpose:
18) ! map input data given on single, discrete points onto a set of target
19) ! bins.
20) ! the original input data are given on single, discrete points of an
21) ! arbitrary grid and are being linearly interpolated onto a specified set
22) ! of target bins. in general, this is the case for most of the weighting
23) ! functions (action spectra, molecular cross section, and quantum yield
24) ! data), which have to be matched onto the specified wavelength intervals.
25) ! the average value in each target bin is found by averaging the trapezoi-
26) ! dal area underneath the input data curve (constructed by linearly connec-
27) ! ting the discrete input values).
28) ! some caution should be used near the endpoints of the grids. if the
29) ! input data set does not span the range of the target grid, an error
30) ! message is printed and the execution is stopped, as extrapolation of the
31) ! data is not permitted.
32) ! if the input data does not encompass the target grid, use addpnt to
33) ! expand the input array.
34) !-----------------------------------------------------------------------------
35) ! parameters:
36) ! ng - integer, number of bins + 1 in the target grid (i)
37) ! xg - real, target grid (e.g., wavelength grid); bin i is defined (i)
38) ! as [xg(i),xg(i+1)] (i = 1..ng-1)
39) ! yg - real, y-data re-gridded onto xg, yg(i) specifies the value for (o)
40) ! bin i (i = 1..ng-1)
41) ! n - integer, number of points in input grid (i)
42) ! x - real, grid on which input data are defined (i)
43) ! y - real, input y-data (i)
44) !-----------------------------------------------------------------------------
45)
46) implicit none
47)
48) !-----------------------------------------------------------------------------
49) ! ... dummy arguments
50) !-----------------------------------------------------------------------------
51) integer, intent(in) :: ng, n
52) integer, intent(out) :: ierr
53) real(r8), intent(in) :: x(n)
54) real(r8), intent(in) :: y(n)
55) real(r8), intent(in) :: xg(ng)
56) real(r8), intent(out) :: yg(ng)
57)
58) !-----------------------------------------------------------------------------
59) ! ... local variables
60) !-----------------------------------------------------------------------------
61) integer :: ngintv
62) integer :: i, k, jstart
63) real(r8) :: area, xgl, xgu
64) real(r8) :: darea, slope
65) real(r8) :: a1, a2, b1, b2
66)
67) ierr = 0
68) !-----------------------------------------------------------------------------
69) ! ... test for correct ordering of data, by increasing value of x
70) !-----------------------------------------------------------------------------
71) do i = 2,n
72) if( x(i) <= x(i-1) ) then
73) ierr = 1
74) write(iulog,*) 'inter2: x coord not monotonically increasing'
75) return
76) end if
77) end do
78)
79) do i = 2,ng
80) if( xg(i) <= xg(i-1) ) then
81) ierr = 2
82) write(iulog,*) 'inter2: xg coord not monotonically increasing'
83) return
84) end if
85) end do
86)
87) !-----------------------------------------------------------------------------
88) ! ... check for xg-values outside the x-range
89) !-----------------------------------------------------------------------------
90) if( x(1) > xg(1) .or. x(n) < xg(ng) ) then
91) write(iulog,*) 'inter2: data does not span grid'
92) write(iulog,*) ' use addpnt to expand data and re-run'
93) ! call endrun
94) end if
95)
96) !-----------------------------------------------------------------------------
97) ! ... find the integral of each grid interval and use this to
98) ! calculate the average y value for the interval
99) ! xgl and xgu are the lower and upper limits of the grid interval
100) !-----------------------------------------------------------------------------
101) jstart = 1
102) ngintv = ng - 1
103) do i = 1,ngintv
104) !-----------------------------------------------------------------------------
105) ! ... initalize:
106) !-----------------------------------------------------------------------------
107) area = 0._r8
108) xgl = xg(i)
109) xgu = xg(i+1)
110) !-----------------------------------------------------------------------------
111) ! ... discard data before the first grid interval and after the
112) ! last grid interval
113) ! for internal grid intervals, start calculating area by interpolating
114) ! between the last point which lies in the previous interval and the
115) ! first point inside the current interval
116) !-----------------------------------------------------------------------------
117) k = jstart
118) if( k <= n-1 ) then
119) !-----------------------------------------------------------------------------
120) ! ... if both points are before the first grid, go to the next point
121) !-----------------------------------------------------------------------------
122) do
123) if( x(k+1) <= xgl ) then
124) jstart = k - 1
125) k = k+1
126) if( k <= n-1 ) then
127) cycle
128) else
129) exit
130) end if
131) else
132) exit
133) end if
134) end do
135) !-----------------------------------------------------------------------------
136) ! ... if the last point is beyond the end of the grid,
137) ! complete and go to the next grid
138) !-----------------------------------------------------------------------------
139) do
140) if( k <= n-1 .and. x(k) < xgu ) then
141) jstart = k-1
142) !-----------------------------------------------------------------------------
143) ! ... compute x-coordinates of increment
144) !-----------------------------------------------------------------------------
145) a1 = max( x(k),xgl )
146) a2 = min( x(k+1),xgu )
147) !-----------------------------------------------------------------------------
148) ! ... if points coincide, contribution is zero
149) !-----------------------------------------------------------------------------
150) if( x(k+1) == x(k) ) then
151) darea = 0._r8
152) else
153) slope = (y(k+1) - y(k))/(x(k+1) - x(k))
154) b1 = y(k) + slope*(a1 - x(k))
155) b2 = y(k) + slope*(a2 - x(k))
156) darea = .5_r8*(a2 - a1)*(b2 + b1)
157) end if
158) !-----------------------------------------------------------------------------
159) ! ... find the area under the trapezoid from a1 to a2
160) !-----------------------------------------------------------------------------
161) area = area + darea
162) k = k+1
163) cycle
164) else
165) exit
166) end if
167) end do
168) end if
169) !-----------------------------------------------------------------------------
170) ! ... calculate the average y after summing the areas in the interval
171) !-----------------------------------------------------------------------------
172) yg(i) = area/(xgu - xgl)
173) end do
174)
175) end subroutine inter2
176)
177) subroutine inter_inti( ng, xg, n, x )
178) !-----------------------------------------------------------------------------
179) ! ... initialization
180) !-----------------------------------------------------------------------------
181)
182) use cam_abortutils, only : endrun
183)
184) implicit none
185)
186) !-----------------------------------------------------------------------------
187) ! ... dummy arguments
188) !-----------------------------------------------------------------------------
189) integer, intent(in) :: ng, n
190) real(r8), intent(in) :: xg(ng)
191) real(r8), intent(in) :: x(n)
192)
193) !-----------------------------------------------------------------------------
194) ! ... local variables
195) !-----------------------------------------------------------------------------
196) integer :: i, ii, iil, astat
197) integer :: ndim(1)
198) #ifdef DEBUG
199) write(iulog,*) 'inter3: diagnostics; ng,n = ',ng,n
200) write(iulog,'('' xg '' )')
201) write(iulog,'(1p,5e21.13)') xg(:ng)
202) write(iulog,'('' x '' )')
203) write(iulog,'(1p,5e21.13)') x(:n)
204) #endif
205) allocate( xi(ng), xcnt(ng-1), stat=astat )
206) if( astat /= 0 ) then
207) write(iulog,*) 'inter_inti: failed to allocate wrk arrays; error = ',astat
208) call endrun
209) else
210) xi(:) = 0
211) xcnt(:) = 0
212) end if
213) iil = 1
214) do i = 1,ng
215) do ii = iil,n-1
216) if( xg(i) < x(ii) ) then
217) xi(i) = ii - 1
218) iil = ii
219) exit
220) end if
221) end do
222) end do
223) nintervals = count( xi(:) /= 0 )
224) if( nintervals == 0 ) then
225) write(iulog,*) 'inter_inti: wavelength grids do not overlap'
226) call endrun
227) else
228) nintervals = nintervals - 1
229) end if
230) xcnt(1:nintervals) = xi(2:nintervals+1) - xi(1:nintervals) + 1
231) ndim(:) = maxval( xcnt(1:nintervals) )
232) allocate( xfrac(ndim(1),nintervals),stat=astat )
233) if( astat /= 0 ) then
234) write(iulog,*) 'inter_inti: failed to allocate wrk array; error = ',astat
235) call endrun
236) else
237) xfrac(:,:) = 1._r8
238) end if
239) do i = 1,nintervals
240) iil = xi(i)
241) xfrac(1,i) = (min( x(iil+1),xg(i+1) ) - xg(i))/(x(iil+1) - x(iil))
242) if( xcnt(i) > 1 ) then
243) iil = xi(i) + xcnt(i) - 1
244) xfrac(xcnt(i),i) = (xg(i+1) - x(iil))/(x(iil+1) - x(iil))
245) end if
246) end do
247) write(iulog,*) 'inter_inti: diagnostics; ng,n,nintervals = ',ng,n,nintervals
248) write(iulog,'('' xi'')')
249) write(iulog,'(10i4)') xi(1:nintervals+1)
250) write(iulog,'('' xcnt'')')
251) write(iulog,'(10i4)') xcnt(1:nintervals)
252) write(iulog,'('' xfrac'')')
253) do i = 1,nintervals
254) write(iulog,'(1p,5e21.13)') xfrac(1:xcnt(i),i)
255) end do
256)
257) end subroutine inter_inti
258)
259) subroutine inter3( ng, xg, yg, n, x, y )
260) !-----------------------------------------------------------------------------
261) ! purpose:
262) ! map input data given on a set of bins onto a different set of target
263) ! bins.
264) ! the input data are given on a set of bins (representing the integral
265) ! of the input quantity over the range of each bin) and are being matched
266) ! onto another set of bins (target grid). a typical example would be an
267) ! input data set spcifying the extra-terrestrial flux on wavelength inter-
268) ! vals, that has to be matched onto the working wavelength grid.
269) ! the resulting area in a given bin of the target grid is calculated by
270) ! simply adding all fractional areas of the input data that cover that
271) ! particular target bin.
272) ! some caution should be used near the endpoints of the grids. if the
273) ! input data do not span the full range of the target grid, the area in
274) ! the "missing" bins will be assumed to be zero. if the input data extend
275) ! beyond the upper limit of the target grid, the user has the option to
276) ! integrate the "overhang" data and fold the remaining area back into the
277) ! last target bin. using this option is recommended when re-gridding
278) ! vertical profiles that directly affect the total optical depth of the
279) ! model atmosphere.
280) !-----------------------------------------------------------------------------
281) ! parameters:
282) ! ng - integer, number of bins + 1 in the target grid (i)
283) ! xg - real, target grid (e.g. working wavelength grid); bin i (i)
284) ! is defined as [xg(i),xg(i+1)] (i = 1..ng-1)
285) ! yg - real, y-data re-gridded onto xg; yg(i) specifies the (o)
286) ! y-value for bin i (i = 1..ng-1)
287) ! n - integer, number of bins + 1 in the input grid (i)
288) ! x - real, input grid (e.g. data wavelength grid); bin i is (i)
289) ! defined as [x(i),x(i+1)] (i = 1..n-1)
290) ! y - real, input y-data on grid x; y(i) specifies the (i)
291) ! y-value for bin i (i = 1..n-1)
292) !-----------------------------------------------------------------------------
293)
294) implicit none
295)
296) !-----------------------------------------------------------------------------
297) ! ... dummy arguments
298) !-----------------------------------------------------------------------------
299) integer, intent(in) :: n, ng
300) real(r8), intent(in) :: xg(ng)
301) real(r8), intent(in) :: x(n)
302) real(r8), intent(in) :: y(n)
303) real(r8), intent(out) :: yg(ng)
304)
305) !-----------------------------------------------------------------------------
306) ! ... local variables
307) !-----------------------------------------------------------------------------
308) integer :: i, ii, iil
309)
310) !-----------------------------------------------------------------------------
311) ! ... do interpolation
312) !-----------------------------------------------------------------------------
313) yg(:) = 0._r8
314) do i = 1,nintervals
315) iil = xi(i)
316) ii = xcnt(i)
317) if( ii == 1 ) then
318) yg(i) = xfrac(1,i)*y(iil)
319) else
320) yg(i) = dot_product( xfrac(1:ii,i),y(iil:iil+ii-1) )
321) end if
322) end do
323) #ifdef DEBUG
324) write(iulog,'('' y '')')
325) write(iulog,'(1p,5e21.13)') y
326) write(iulog,'('' yg '')')
327) write(iulog,'(1p,5e21.13)') yg(1:nintervals)
328) #endif
329)
330) end subroutine inter3
331)
332) end module mo_inter
