* 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 20, 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 20 CONTINUE *_______________________________________________________________________ RETURN END *=============================================================================* 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 10, i = 2, n IF (x(i) .LE. x(i-1)) THEN ierr = 1 WRITE(*,*)'data not sorted' RETURN ENDIF 10 CONTINUE DO i = 2, ng IF (xg(i) .LE. xg(i-1)) THEN ierr = 2 WRITE(0,*) '>>> ERROR (inter2) <<< 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 WRITE(0,*) '>>> ERROR (inter2) <<< Data do not span '// > 'grid. ' WRITE(0,*) ' Use ADDPNT to '// > 'expand data and re-run.' STOP 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 50, 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) 50 CONTINUE *_______________________________________________________________________ RETURN END *=============================================================================* 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 WRITE(0,*) '>>> ERROR (inter3) <<< Value for FOLDIN invalid. ' WRITE(0,*) ' Must be 0 or 1' STOP ENDIF * do interpolation jstart = 1 DO 30, 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 = AMAX1(x(j),xg(i)) a2 = AMIN1(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 30 CONTINUE * 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 *_______________________________________________________________________ RETURN END *=============================================================================* 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 WRITE(0,*) '>>> ERROR (inter3) <<< Value for FOLDIN invalid. ' WRITE(0,*) ' Must be 0 or 1' STOP ENDIF * do interpolation jstart = 1 DO 30, 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 = AMAX1(x(j),xg(i)) a2 = AMIN1(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 30 CONTINUE * 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 *_______________________________________________________________________ RETURN END *=============================================================================* SUBROUTINE addpnt ( x, y, ld, n, xnew, ynew ) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Add a point to a set of data pairs . 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 ld, n REAL x(ld), y(ld) REAL xnew, ynew * local variables INTEGER insert INTEGER i *----------------------------------------------------------------------- * check n>> ERROR (ADDPNT) <<< Cannot expand array ' WRITE(0,*) ' All elements used.' STOP 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 WRITE(0,*) '>>> ERROR (ADDPNT) <<< x-data must be '// > 'in ascending order!' STOP ELSE IF (xnew .GT. x(i)) insert = i + 1 ENDIF i = i+1 GOTO 10 ENDIF * if needs to be appended at the end, just do so, * otherwise, insert 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 zero1(x,m) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Initialize all elements of a floating point vector with zero. =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= X - REAL, vector to be initialized (O)=* *= M - INTEGER, number of elements in X (I)=* *-----------------------------------------------------------------------------* IMPLICIT NONE INTEGER i, m REAL x(m) DO 1 i = 1, m x(i) = 0. 1 CONTINUE RETURN END *=============================================================================* SUBROUTINE zero2(x,m,n) *-----------------------------------------------------------------------------* *= PURPOSE: =* *= Initialize all elements of a 2D floating point array with zero. =* *-----------------------------------------------------------------------------* *= PARAMETERS: =* *= X - REAL, array to be initialized (O)=* *= M - INTEGER, number of elements along the first dimension of X, (I)=* *= exactly as specified in the calling program =* *= N - INTEGER, number of elements along the second dimension of X, (I)=* *= exactly as specified in the calling program =* *-----------------------------------------------------------------------------* IMPLICIT none * m,n : dimensions of x, exactly as specified in the calling program INTEGER i, j, m, n REAL x(m,n) DO 1 j = 1, n DO 2 i = 1, m x(i,j) = 0. 2 CONTINUE 1 CONTINUE RETURN END