! This file contains the following subroutines, related to interpolations
! of input data, addition of points to arrays, and zeroing of arrays:
! inter1
! inter2
! inter3
! inter4
! addpnt
! zero1
! zero2
!=============================================================================*
SUBROUTINE inter1(ng,xg,yg, n,x,y)
!-----------------------------------------------------------------------------*
!= PURPOSE: =*
!= Map input data given on single, discrete points, onto a discrete target =*
!= grid. =*
!= The original input data are given on single, discrete points of an =*
!= arbitrary grid and are being linearly interpolated onto a specified =*
!= discrete target grid. A typical example would be the re-gridding of a =*
!= given data set for the vertical temperature profile to match the speci- =*
!= fied altitude grid. =*
!= Some caution should be used near the end points of the grids. If the =*
!= input data set does not span the range of the target grid, the remaining =*
!= points will be set to zero, as extrapolation is not permitted. =*
!= If the input data does not encompass the target grid, use ADDPNT to =*
!= expand the input array. =*
!-----------------------------------------------------------------------------*
!= PARAMETERS: =*
!= NG - INTEGER, number of points in the target grid (I)=*
!= XG - REAL, target grid (e.g. altitude grid) (I)=*
!= YG - REAL, y-data re-gridded onto XG (O)=*
!= N - INTEGER, number of points in the input data set (I)=*
!= X - REAL, grid on which input data are defined (I)=*
!= Y - REAL, input y-data (I)=*
!-----------------------------------------------------------------------------*
IMPLICIT NONE
! input:
INTEGER n, ng
REAL xg(ng)
REAL x(n), y(n)
! output:
REAL yg(ng)
! local:
REAL slope
INTEGER jsave, i, j
jsave = 1
DO i = 1, ng
yg(i) = 0.
j = jsave
10 CONTINUE
IF ((x(j) .GT. xg(i)) .OR. (xg(i) .GE. x(j+1))) THEN
j = j+1
IF (j .LE. n-1) GOTO 10
! ---- end of loop 10 ----
ELSE
slope = (y(j+1)-y(j)) / (x(j+1)-x(j))
yg(i) = y(j) + slope * (xg(i) - x(j))
jsave = j
ENDIF
ENDDO
END SUBROUTINE inter1
!=============================================================================*
SUBROUTINE inter2(ng,xg,yg,n,x,y,ierr)
!-----------------------------------------------------------------------------*
!= PURPOSE: =*
!= Map input data given on single, discrete points onto a set of target =*
!= bins. =*
!= The original input data are given on single, discrete points of an =*
!= arbitrary grid and are being linearly interpolated onto a specified set =*
!= of target bins. In general, this is the case for most of the weighting =*
!= functions (action spectra, molecular cross section, and quantum yield =*
!= data), which have to be matched onto the specified wavelength intervals. =*
!= The average value in each target bin is found by averaging the trapezoi- =*
!= dal area underneath the input data curve (constructed by linearly connec-=*
!= ting the discrete input values). =*
!= Some caution should be used near the endpoints of the grids. If the =*
!= input data set does not span the range of the target grid, an error =*
!= message is printed and the execution is stopped, as extrapolation of the =*
!= data is not permitted. =*
!= If the input data does not encompass the target grid, use ADDPNT to =*
!= expand the input array. =*
!-----------------------------------------------------------------------------*
!= PARAMETERS: =*
!= NG - INTEGER, number of bins + 1 in the target grid (I)=*
!= XG - REAL, target grid (e.g., wavelength grid); bin i is defined (I)=*
!= as [XG(i),XG(i+1)] (i = 1..NG-1) =*
!= YG - REAL, y-data re-gridded onto XG, YG(i) specifies the value for (O)=*
!= bin i (i = 1..NG-1) =*
!= N - INTEGER, number of points in input grid (I)=*
!= X - REAL, grid on which input data are defined (I)=*
!= Y - REAL, input y-data (I)=*
!-----------------------------------------------------------------------------*
IMPLICIT NONE
! input:
INTEGER ng, n
REAL x(n), y(n), xg(ng)
! output:
REAL yg(ng)
! local:
REAL area, xgl, xgu
REAL darea, slope
REAL a1, a2, b1, b2
INTEGER ngintv
INTEGER i, k, jstart
INTEGER ierr
ierr = 0
! test for correct ordering of data, by increasing value of x
DO i = 2, n
IF (x(i) .LE. x(i-1)) THEN
ierr = 1
call wrf_debug( 0,'inter2: ERROR <<< x-grid not sorted' )
RETURN
ENDIF
ENDDO
DO i = 2, ng
IF (xg(i) .LE. xg(i-1)) THEN
ierr = 2
call wrf_debug( 0,'inter2: ERROR <<< xg-grid not sorted!' )
RETURN
ENDIF
ENDDO
! check for xg-values outside the x-range
IF ( (x(1) .GT. xg(1)) .OR. (x(n) .LT. xg(ng)) ) THEN
call wrf_error_fatal( 'inter2: <<< Data do not span grid; Use ADDPNT to expand data and re-run.' )
ENDIF
! find the integral of each grid interval and use this to
! calculate the average y value for the interval
! xgl and xgu are the lower and upper limits of the grid interval
jstart = 1
ngintv = ng - 1
DO i = 1,ngintv
! initalize:
area = 0.0
xgl = xg(i)
xgu = xg(i+1)
! discard data before the first grid interval and after the
! last grid interval
! for internal grid intervals, start calculating area by interpolating
! between the last point which lies in the previous interval and the
! first point inside the current interval
k = jstart
IF (k .LE. n-1) THEN
! if both points are before the first grid, go to the next point
30 CONTINUE
IF (x(k+1) .LE. xgl) THEN
jstart = k - 1
k = k+1
IF (k .LE. n-1) GO TO 30
ENDIF
! if the last point is beyond the end of the grid, complete and go to the next
! grid
40 CONTINUE
IF ((k .LE. n-1) .AND. (x(k) .LT. xgu)) THEN
jstart = k-1
! compute x-coordinates of increment
a1 = MAX(x(k),xgl)
a2 = MIN(x(k+1),xgu)
! if points coincide, contribution is zero
IF (x(k+1).EQ.x(k)) THEN
darea = 0.e0
ELSE
slope = (y(k+1) - y(k))/(x(k+1) - x(k))
b1 = y(k) + slope*(a1 - x(k))
b2 = y(k) + slope*(a2 - x(k))
darea = (a2 - a1)*(b2 + b1)/2.
ENDIF
! find the area under the trapezoid from a1 to a2
area = area + darea
! go to next point
k = k+1
GO TO 40
ENDIF
ENDIF
! calculate the average y after summing the areas in the interval
yg(i) = area/(xgu - xgl)
ENDDO
END SUBROUTINE inter2
!=============================================================================*
SUBROUTINE inter3(ng,xg,yg, n,x,y, FoldIn)
!-----------------------------------------------------------------------------*
!= PURPOSE: =*
!= Map input data given on a set of bins onto a different set of target =*
!= bins. =*
!= The input data are given on a set of bins (representing the integral =*
!= of the input quantity over the range of each bin) and are being matched =*
!= onto another set of bins (target grid). A typical example would be an =*
!= input data set spcifying the extra-terrestrial flux on wavelength inter- =*
!= vals, that has to be matched onto the working wavelength grid. =*
!= The resulting area in a given bin of the target grid is calculated by =*
!= simply adding all fractional areas of the input data that cover that =*
!= particular target bin. =*
!= Some caution should be used near the endpoints of the grids. If the =*
!= input data do not span the full range of the target grid, the area in =*
!= the "missing" bins will be assumed to be zero. If the input data extend =*
!= beyond the upper limit of the target grid, the user has the option to =*
!= integrate the "overhang" data and fold the remaining area back into the =*
!= last target bin. Using this option is recommended when re-gridding =*
!= vertical profiles that directly affect the total optical depth of the =*
!= model atmosphere. =*
!-----------------------------------------------------------------------------*
!= PARAMETERS: =*
!= NG - INTEGER, number of bins + 1 in the target grid (I)=*
!= XG - REAL, target grid (e.g. working wavelength grid); bin i (I)=*
!= is defined as [XG(i),XG(i+1)] (i = 1..NG-1) =*
!= YG - REAL, y-data re-gridded onto XG; YG(i) specifies the (O)=*
!= y-value for bin i (i = 1..NG-1) =*
!= N - INTEGER, number of bins + 1 in the input grid (I)=*
!= X - REAL, input grid (e.g. data wavelength grid); bin i is (I)=*
!= defined as [X(i),X(i+1)] (i = 1..N-1) =*
!= Y - REAL, input y-data on grid X; Y(i) specifies the (I)=*
!= y-value for bin i (i = 1..N-1) =*
!= FoldIn - Switch for folding option of "overhang" data (I)=*
!= FoldIn = 0 -> No folding of "overhang" data =*
!= FoldIn = 1 -> Integerate "overhang" data and fold back into =*
!= last target bin =*
!-----------------------------------------------------------------------------*
IMPLICIT NONE
! input:
INTEGER n, ng
REAL xg(ng)
REAL x(n), y(n)
INTEGER FoldIn
! output:
REAL yg(ng)
! local:
REAL a1, a2, sum
REAL tail
INTEGER jstart, i, j, k
! check whether flag given is legal
IF ((FoldIn .NE. 0) .AND. (FoldIn .NE. 1)) THEN
call wrf_error_fatal( 'inter3: ERROR <<< Value for FOLDIN invalid. Must be 0 or 1' )
ENDIF
! do interpolation
jstart = 1
DO i = 1, ng - 1
yg(i) = 0.
sum = 0.
j = jstart
IF (j .LE. n-1) THEN
20 CONTINUE
IF (x(j+1) .LT. xg(i)) THEN
jstart = j
j = j+1
IF (j .LE. n-1) GO TO 20
ENDIF
25 CONTINUE
IF ((x(j) .LE. xg(i+1)) .AND. (j .LE. n-1)) THEN
a1 = MAX(x(j),xg(i))
a2 = MIN(x(j+1),xg(i+1))
sum = sum + y(j) * (a2-a1)/(x(j+1)-x(j))
j = j+1
GO TO 25
ENDIF
yg(i) = sum
ENDIF
ENDDO
! if wanted, integrate data "overhang" and fold back into last bin
IF (FoldIn .EQ. 1) THEN
j = j-1
a1 = xg(ng) ! upper limit of last interpolated bin
a2 = x(j+1) ! upper limit of last input bin considered
! do folding only if grids don't match up and there is more input
IF ((a2 .GT. a1) .OR. (j+1 .LT. n)) THEN
tail = y(j) * (a2-a1)/(x(j+1)-x(j))
DO k = j+1, n-1
tail = tail + y(k) * (x(k+1)-x(k))
ENDDO
yg(ng-1) = yg(ng-1) + tail
ENDIF
ENDIF
END SUBROUTINE inter3
!=============================================================================*
SUBROUTINE inter4(ng,xg,yg, n,x,y, FoldIn)
!-----------------------------------------------------------------------------*
!= PURPOSE: =*
!= Map input data given on a set of bins onto a different set of target =*
!= bins. =*
!= The input data are given on a set of bins (representing the integral =*
!= of the input quantity over the range of each bin) and are being matched =*
!= onto another set of bins (target grid). A typical example would be an =*
!= input data set spcifying the extra-terrestrial flux on wavelength inter- =*
!= vals, that has to be matched onto the working wavelength grid. =*
!= The resulting area in a given bin of the target grid is calculated by =*
!= simply adding all fractional areas of the input data that cover that =*
!= particular target bin. =*
!= Some caution should be used near the endpoints of the grids. If the =*
!= input data do not span the full range of the target grid, the area in =*
!= the "missing" bins will be assumed to be zero. If the input data extend =*
!= beyond the upper limit of the target grid, the user has the option to =*
!= integrate the "overhang" data and fold the remaining area back into the =*
!= last target bin. Using this option is recommended when re-gridding =*
!= vertical profiles that directly affect the total optical depth of the =*
!= model atmosphere. =*
!-----------------------------------------------------------------------------*
!= PARAMETERS: =*
!= NG - INTEGER, number of bins + 1 in the target grid (I)=*
!= XG - REAL, target grid (e.g. working wavelength grid); bin i (I)=*
!= is defined as [XG(i),XG(i+1)] (i = 1..NG-1) =*
!= YG - REAL, y-data re-gridded onto XG; YG(i) specifies the (O)=*
!= y-value for bin i (i = 1..NG-1) =*
!= N - INTEGER, number of bins + 1 in the input grid (I)=*
!= X - REAL, input grid (e.g. data wavelength grid); bin i is (I)=*
!= defined as [X(i),X(i+1)] (i = 1..N-1) =*
!= Y - REAL, input y-data on grid X; Y(i) specifies the (I)=*
!= y-value for bin i (i = 1..N-1) =*
!= FoldIn - Switch for folding option of "overhang" data (I)=*
!= FoldIn = 0 -> No folding of "overhang" data =*
!= FoldIn = 1 -> Integerate "overhang" data and fold back into =*
!= last target bin =*
!-----------------------------------------------------------------------------*
IMPLICIT NONE
! input:
INTEGER n, ng
REAL xg(ng)
REAL x(n), y(n)
INTEGER FoldIn
! output:
REAL yg(ng)
! local:
REAL a1, a2, sum
REAL tail
INTEGER jstart, i, j, k
! check whether flag given is legal
IF ((FoldIn .NE. 0) .AND. (FoldIn .NE. 1)) THEN
call wrf_error_fatal( 'inter3: ERROR <<< Value for FOLDIN invalid. Must be 0 or 1' )
ENDIF
! do interpolation
jstart = 1
DO i = 1, ng - 1
yg(i) = 0.
sum = 0.
j = jstart
IF (j .LE. n-1) THEN
20 CONTINUE
IF (x(j+1) .LT. xg(i)) THEN
jstart = j
j = j+1
IF (j .LE. n-1) GO TO 20
ENDIF
25 CONTINUE
IF ((x(j) .LE. xg(i+1)) .AND. (j .LE. n-1)) THEN
a1 = MAX(x(j),xg(i))
a2 = MIN(x(j+1),xg(i+1))
sum = sum + y(j) * (a2-a1)
j = j+1
GO TO 25
ENDIF
yg(i) = sum /(xg(i+1)-xg(i))
ENDIF
ENDDO
! if wanted, integrate data "overhang" and fold back into last bin
IF (FoldIn .EQ. 1) THEN
j = j-1
a1 = xg(ng) ! upper limit of last interpolated bin
a2 = x(j+1) ! upper limit of last input bin considered
! do folding only if grids don't match up and there is more input
IF ((a2 .GT. a1) .OR. (j+1 .LT. n)) THEN
tail = y(j) * (a2-a1)/(x(j+1)-x(j))
DO k = j+1, n-1
tail = tail + y(k) * (x(k+1)-x(k))
ENDDO
yg(ng-1) = yg(ng-1) + tail
ENDIF
ENDIF
END SUBROUTINE inter4
!=============================================================================*
SUBROUTINE addpnt ( x, y, ld, n, xnew, ynew )
!-----------------------------------------------------------------------------*
!= PURPOSE: =*
!= Add a point <xnew,ynew> to a set of data pairs <x,y>. x must be in =*
!= ascending order =*
!-----------------------------------------------------------------------------*
!= PARAMETERS: =*
!= X - REAL vector of length LD, x-coordinates (IO)=*
!= Y - REAL vector of length LD, y-values (IO)=*
!= LD - INTEGER, dimension of X, Y exactly as declared in the calling (I)=*
!= program =*
!= N - INTEGER, number of elements in X, Y. On entry, it must be: (IO)=*
!= N < LD. On exit, N is incremented by 1. =*
!= XNEW - REAL, x-coordinate at which point is to be added (I)=*
!= YNEW - REAL, y-value of point to be added (I)=*
!-----------------------------------------------------------------------------*
IMPLICIT NONE
! calling parameters
INTEGER, intent(in) :: ld
INTEGER, intent(inout) :: n
REAL, intent(inout) :: x(ld), y(ld)
REAL, intent(in) :: xnew, ynew
! local variables
INTEGER insert
INTEGER i
CHARACTER(len=256) :: emsg
! check n<ld to make sure x will hold another point
IF (n .GE. ld) THEN
call wrf_error_fatal( 'addpnt: ERROR <<< Cannot expand array All elements used.' )
ENDIF
insert = 1
i = 2
! check whether x is already sorted.
! also, use this loop to find the point at which xnew needs to be inserted
! into vector x, if x is sorted.
10 CONTINUE
IF (i .LT. n) THEN
IF (x(i) .LT. x(i-1)) THEN
call wrf_error_fatal( 'addpnt: ERROR <<< x-data must be in ascending order!' )
ELSE
IF (xnew .GT. x(i-1)) insert = i
ENDIF
i = i+1
GOTO 10
ENDIF
! if <xnew,ynew> needs to be appended at the end, just do so,
! otherwise, insert <xnew,ynew> at position INSERT
IF ( xnew .GT. x(n) ) THEN
x(n+1) = xnew
y(n+1) = ynew
ELSE
! shift all existing points one index up
DO i = n, insert, -1
x(i+1) = x(i)
y(i+1) = y(i)
ENDDO
! insert new point
x(insert) = xnew
y(insert) = ynew
ENDIF
! increase total number of elements in x, y
n = n+1
END SUBROUTINE addpnt