chem/numer.F

   1 ! This file contains the following subroutines, related to interpolations
   2 ! of input data, addition of points to arrays, and zeroing of arrays:
   3 !     inter1
   4 !     inter2
   5 !     inter3
   6 !     inter4
   7 !     addpnt
   8 !     zero1
   9 !     zero2
  10 !=============================================================================*
  11
  12       SUBROUTINE inter1(ng,xg,yg, n,x,y)
  13 !-----------------------------------------------------------------------------*
  14 !=  PURPOSE:                                                                 =*
  15 !=  Map input data given on single, discrete points, onto a discrete target  =*
  16 !=  grid.                                                                    =*
  17 !=  The original input data are given on single, discrete points of an       =*
  18 !=  arbitrary grid and are being linearly interpolated onto a specified      =*
  19 !=  discrete target grid.  A typical example would be the re-gridding of a   =*
  20 !=  given data set for the vertical temperature profile to match the speci-  =*
  21 !=  fied altitude grid.                                                      =*
  22 !=  Some caution should be used near the end points of the grids.  If the    =*
  23 !=  input data set does not span the range of the target grid, the remaining =*
  24 !=  points will be set to zero, as extrapolation is not permitted.           =*
  25 !=  If the input data does not encompass the target grid, use ADDPNT to      =*
  26 !=  expand the input array.                                                  =*
  27 !-----------------------------------------------------------------------------*
  28 !=  PARAMETERS:                                                              =*
  29 !=  NG  - INTEGER, number of points in the target grid                    (I)=*
  30 !=  XG  - REAL, target grid (e.g. altitude grid)                          (I)=*
  31 !=  YG  - REAL, y-data re-gridded onto XG                                 (O)=*
  32 !=  N   - INTEGER, number of points in the input data set                 (I)=*
  33 !=  X   - REAL, grid on which input data are defined                      (I)=*
  34 !=  Y   - REAL, input y-data                                              (I)=*
  35 !-----------------------------------------------------------------------------*
  36
  37       IMPLICIT NONE
  38
  39 ! input:
  40       INTEGER n, ng
  41       REAL xg(ng)
  42       REAL x(n), y(n)
  43
  44 ! output:
  45       REAL yg(ng)
  46
  47 ! local:
  48       REAL slope
  49       INTEGER jsave, i, j
  50
  51       jsave = 1
  52       DO i = 1, ng
  53          yg(i) = 0.
  54          j = jsave
  55    10    CONTINUE
  56          IF ((x(j) .GT. xg(i)) .OR. (xg(i) .GE. x(j+1))) THEN
  57             j = j+1
  58             IF (j .LE. n-1) GOTO 10
  59 !        ---- end of loop 10 ----
  60          ELSE
  61             slope = (y(j+1)-y(j)) / (x(j+1)-x(j))
  62             yg(i) = y(j) + slope * (xg(i) - x(j))
  63             jsave = j
  64          ENDIF
  65       ENDDO
  66
  67       END SUBROUTINE inter1
  68
  69 !=============================================================================*
  70
  71       SUBROUTINE inter2(ng,xg,yg,n,x,y,ierr)
  72
  73 !-----------------------------------------------------------------------------*
  74 !=  PURPOSE:                                                                 =*
  75 !=  Map input data given on single, discrete points onto a set of target     =*
  76 !=  bins.                                                                    =*
  77 !=  The original input data are given on single, discrete points of an       =*
  78 !=  arbitrary grid and are being linearly interpolated onto a specified set  =*
  79 !=  of target bins.  In general, this is the case for most of the weighting  =*
  80 !=  functions (action spectra, molecular cross section, and quantum yield    =*
  81 !=  data), which have to be matched onto the specified wavelength intervals. =*
  82 !=  The average value in each target bin is found by averaging the trapezoi- =*
  83 !=  dal area underneath the input data curve (constructed by linearly connec-=*
  84 !=  ting the discrete input values).                                         =*
  85 !=  Some caution should be used near the endpoints of the grids.  If the     =*
  86 !=  input data set does not span the range of the target grid, an error      =*
  87 !=  message is printed and the execution is stopped, as extrapolation of the =*
  88 !=  data is not permitted.                                                   =*
  89 !=  If the input data does not encompass the target grid, use ADDPNT to      =*
  90 !=  expand the input array.                                                  =*
  91 !-----------------------------------------------------------------------------*
  92 !=  PARAMETERS:                                                              =*
  93 !=  NG  - INTEGER, number of bins + 1 in the target grid                  (I)=*
  94 !=  XG  - REAL, target grid (e.g., wavelength grid);  bin i is defined    (I)=*
  95 !=        as [XG(i),XG(i+1)] (i = 1..NG-1)                                   =*
  96 !=  YG  - REAL, y-data re-gridded onto XG, YG(i) specifies the value for  (O)=*
  97 !=        bin i (i = 1..NG-1)                                                =*
  98 !=  N   - INTEGER, number of points in input grid                         (I)=*
  99 !=  X   - REAL, grid on which input data are defined                      (I)=*
 100 !=  Y   - REAL, input y-data                                              (I)=*
 101 !-----------------------------------------------------------------------------*
 102
 103       IMPLICIT NONE
 104
 105 ! input:
 106       INTEGER ng, n
 107       REAL x(n), y(n), xg(ng)
 108
 109 ! output:
 110       REAL yg(ng)
 111
 112 ! local:
 113       REAL area, xgl, xgu
 114       REAL darea, slope
 115       REAL a1, a2, b1, b2
 116       INTEGER ngintv
 117       INTEGER i, k, jstart
 118       INTEGER ierr
 119
 120       ierr = 0
 121
 122 !  test for correct ordering of data, by increasing value of x
 123
 124       DO i = 2, n
 125         IF (x(i) .LE. x(i-1)) THEN
 126           ierr = 1
 127           call wrf_debug( 0,'inter2: ERROR <<< x-grid not sorted' )
 128           RETURN
 129         ENDIF
 130       ENDDO
 131
 132       DO i = 2, ng
 133         IF (xg(i) .LE. xg(i-1)) THEN
 134           ierr = 2
 135           call wrf_debug( 0,'inter2: ERROR <<< xg-grid not sorted!' )
 136           RETURN
 137         ENDIF
 138       ENDDO
 139
 140 ! check for xg-values outside the x-range
 141
 142       IF ( (x(1) .GT. xg(1)) .OR. (x(n) .LT. xg(ng)) ) THEN
 143         call wrf_error_fatal( 'inter2: <<<  Data do not span grid; Use ADDPNT to expand data and re-run.' )
 144       ENDIF
 145
 146 !  find the integral of each grid interval and use this to
 147 !  calculate the average y value for the interval
 148 !  xgl and xgu are the lower and upper limits of the grid interval
 149
 150       jstart = 1
 151       ngintv = ng - 1
 152       DO i = 1,ngintv
 153
 154 ! initalize:
 155             area = 0.0
 156             xgl = xg(i)
 157             xgu = xg(i+1)
 158
 159 !  discard data before the first grid interval and after the
 160 !  last grid interval
 161 !  for internal grid intervals, start calculating area by interpolating
 162 !  between the last point which lies in the previous interval and the
 163 !  first point inside the current interval
 164
 165             k = jstart
 166             IF (k .LE. n-1) THEN
 167
 168 !  if both points are before the first grid, go to the next point
 169    30         CONTINUE
 170               IF (x(k+1) .LE. xgl) THEN
 171                 jstart = k - 1
 172                 k = k+1
 173                 IF (k .LE. n-1) GO TO 30
 174               ENDIF
 175
 176 !  if the last point is beyond the end of the grid, complete and go to the next
 177 !  grid
 178    40         CONTINUE
 179               IF ((k .LE. n-1) .AND. (x(k) .LT. xgu)) THEN
 180                 jstart = k-1
 181 ! compute x-coordinates of increment
 182                 a1 = MAX(x(k),xgl)
 183                 a2 = MIN(x(k+1),xgu)
 184 !  if points coincide, contribution is zero
 185                 IF (x(k+1).EQ.x(k)) THEN
 186                   darea = 0.e0
 187                 ELSE
 188                   slope = (y(k+1) - y(k))/(x(k+1) - x(k))
 189                   b1 = y(k) + slope*(a1 - x(k))
 190                   b2 = y(k) + slope*(a2 - x(k))
 191                   darea = (a2 - a1)*(b2 + b1)/2.
 192                 ENDIF
 193
 194 !  find the area under the trapezoid from a1 to a2
 195                 area = area + darea
 196 ! go to next point
 197                 k = k+1
 198                 GO TO 40
 199               ENDIF
 200             ENDIF
 201
 202 !  calculate the average y after summing the areas in the interval
 203             yg(i) = area/(xgu - xgl)
 204
 205       ENDDO
 206
 207       END SUBROUTINE inter2
 208
 209 !=============================================================================*
 210
 211       SUBROUTINE inter3(ng,xg,yg, n,x,y, FoldIn)
 212 !-----------------------------------------------------------------------------*
 213 !=  PURPOSE:                                                                 =*
 214 !=  Map input data given on a set of bins onto a different set of target     =*
 215 !=  bins.                                                                    =*
 216 !=  The input data are given on a set of bins (representing the integral     =*
 217 !=  of the input quantity over the range of each bin) and are being matched  =*
 218 !=  onto another set of bins (target grid).  A typical example would be an   =*
 219 !=  input data set spcifying the extra-terrestrial flux on wavelength inter- =*
 220 !=  vals, that has to be matched onto the working wavelength grid.           =*
 221 !=  The resulting area in a given bin of the target grid is calculated by    =*
 222 !=  simply adding all fractional areas of the input data that cover that     =*
 223 !=  particular target bin.                                                   =*
 224 !=  Some caution should be used near the endpoints of the grids.  If the     =*
 225 !=  input data do not span the full range of the target grid, the area in    =*
 226 !=  the "missing" bins will be assumed to be zero.  If the input data extend =*
 227 !=  beyond the upper limit of the target grid, the user has the option to    =*
 228 !=  integrate the "overhang" data and fold the remaining area back into the  =*
 229 !=  last target bin.  Using this option is recommended when re-gridding      =*
 230 !=  vertical profiles that directly affect the total optical depth of the    =*
 231 !=  model atmosphere.                                                        =*
 232 !-----------------------------------------------------------------------------*
 233 !=  PARAMETERS:                                                              =*
 234 !=  NG     - INTEGER, number of bins + 1 in the target grid               (I)=*
 235 !=  XG     - REAL, target grid (e.g. working wavelength grid);  bin i     (I)=*
 236 !=           is defined as [XG(i),XG(i+1)] (i = 1..NG-1)                     =*
 237 !=  YG     - REAL, y-data re-gridded onto XG;  YG(i) specifies the        (O)=*
 238 !=           y-value for bin i (i = 1..NG-1)                                 =*
 239 !=  N      - INTEGER, number of bins + 1 in the input grid                (I)=*
 240 !=  X      - REAL, input grid (e.g. data wavelength grid);  bin i is      (I)=*
 241 !=           defined as [X(i),X(i+1)] (i = 1..N-1)                           =*
 242 !=  Y      - REAL, input y-data on grid X;  Y(i) specifies the            (I)=*
 243 !=           y-value for bin i (i = 1..N-1)                                  =*
 244 !=  FoldIn - Switch for folding option of "overhang" data                 (I)=*
 245 !=           FoldIn = 0 -> No folding of "overhang" data                     =*
 246 !=           FoldIn = 1 -> Integerate "overhang" data and fold back into     =*
 247 !=                         last target bin                                   =*
 248 !-----------------------------------------------------------------------------*
 249
 250       IMPLICIT NONE
 251
 252 ! input:
 253       INTEGER n, ng
 254       REAL xg(ng)
 255       REAL x(n), y(n)
 256
 257       INTEGER FoldIn
 258
 259 ! output:
 260       REAL yg(ng)
 261
 262 ! local:
 263       REAL a1, a2, sum
 264       REAL tail
 265       INTEGER jstart, i, j, k
 266
 267 ! check whether flag given is legal
 268       IF ((FoldIn .NE. 0) .AND. (FoldIn .NE. 1)) THEN
 269          call wrf_error_fatal( 'inter3: ERROR <<<  Value for FOLDIN invalid. Must be 0 or 1' )
 270       ENDIF
 271
 272 ! do interpolation
 273
 274       jstart = 1
 275       DO i = 1, ng - 1
 276          yg(i) = 0.
 277          sum = 0.
 278          j = jstart
 279          IF (j .LE. n-1) THEN
 280    20      CONTINUE
 281
 282            IF (x(j+1) .LT. xg(i)) THEN
 283               jstart = j
 284               j = j+1
 285               IF (j .LE. n-1) GO TO 20
 286            ENDIF
 287
 288    25      CONTINUE
 289
 290            IF ((x(j) .LE. xg(i+1)) .AND. (j .LE. n-1)) THEN
 291               a1 = MAX(x(j),xg(i))
 292               a2 = MIN(x(j+1),xg(i+1))
 293               sum = sum + y(j) * (a2-a1)/(x(j+1)-x(j))
 294               j = j+1
 295               GO TO 25
 296            ENDIF
 297            yg(i) = sum
 298          ENDIF
 299       ENDDO
 300
 301
 302 ! if wanted, integrate data "overhang" and fold back into last bin
 303
 304       IF (FoldIn .EQ. 1) THEN
 305          j = j-1
 306          a1 = xg(ng)            ! upper limit of last interpolated bin
 307          a2 = x(j+1)            ! upper limit of last input bin considered
 308 !        do folding only if grids don't match up and there is more input
 309          IF ((a2 .GT. a1) .OR. (j+1 .LT. n)) THEN
 310             tail = y(j) * (a2-a1)/(x(j+1)-x(j))
 311             DO k = j+1, n-1
 312                tail = tail + y(k) * (x(k+1)-x(k))
 313             ENDDO
 314             yg(ng-1) = yg(ng-1) + tail
 315          ENDIF
 316       ENDIF
 317
 318       END SUBROUTINE inter3
 319
 320 !=============================================================================*
 321
 322       SUBROUTINE inter4(ng,xg,yg, n,x,y, FoldIn)
 323 !-----------------------------------------------------------------------------*
 324 !=  PURPOSE:                                                                 =*
 325 !=  Map input data given on a set of bins onto a different set of target     =*
 326 !=  bins.                                                                    =*
 327 !=  The input data are given on a set of bins (representing the integral     =*
 328 !=  of the input quantity over the range of each bin) and are being matched  =*
 329 !=  onto another set of bins (target grid).  A typical example would be an   =*
 330 !=  input data set spcifying the extra-terrestrial flux on wavelength inter- =*
 331 !=  vals, that has to be matched onto the working wavelength grid.           =*
 332 !=  The resulting area in a given bin of the target grid is calculated by    =*
 333 !=  simply adding all fractional areas of the input data that cover that     =*
 334 !=  particular target bin.                                                   =*
 335 !=  Some caution should be used near the endpoints of the grids.  If the     =*
 336 !=  input data do not span the full range of the target grid, the area in    =*
 337 !=  the "missing" bins will be assumed to be zero.  If the input data extend =*
 338 !=  beyond the upper limit of the target grid, the user has the option to    =*
 339 !=  integrate the "overhang" data and fold the remaining area back into the  =*
 340 !=  last target bin.  Using this option is recommended when re-gridding      =*
 341 !=  vertical profiles that directly affect the total optical depth of the    =*
 342 !=  model atmosphere.                                                        =*
 343 !-----------------------------------------------------------------------------*
 344 !=  PARAMETERS:                                                              =*
 345 !=  NG     - INTEGER, number of bins + 1 in the target grid               (I)=*
 346 !=  XG     - REAL, target grid (e.g. working wavelength grid);  bin i     (I)=*
 347 !=           is defined as [XG(i),XG(i+1)] (i = 1..NG-1)                     =*
 348 !=  YG     - REAL, y-data re-gridded onto XG;  YG(i) specifies the        (O)=*
 349 !=           y-value for bin i (i = 1..NG-1)                                 =*
 350 !=  N      - INTEGER, number of bins + 1 in the input grid                (I)=*
 351 !=  X      - REAL, input grid (e.g. data wavelength grid);  bin i is      (I)=*
 352 !=           defined as [X(i),X(i+1)] (i = 1..N-1)                           =*
 353 !=  Y      - REAL, input y-data on grid X;  Y(i) specifies the            (I)=*
 354 !=           y-value for bin i (i = 1..N-1)                                  =*
 355 !=  FoldIn - Switch for folding option of "overhang" data                 (I)=*
 356 !=           FoldIn = 0 -> No folding of "overhang" data                     =*
 357 !=           FoldIn = 1 -> Integerate "overhang" data and fold back into     =*
 358 !=                         last target bin                                   =*
 359 !-----------------------------------------------------------------------------*
 360
 361       IMPLICIT NONE
 362
 363 ! input:
 364       INTEGER n, ng
 365       REAL xg(ng)
 366       REAL x(n), y(n)
 367
 368       INTEGER FoldIn
 369
 370 ! output:
 371       REAL yg(ng)
 372
 373 ! local:
 374       REAL a1, a2, sum
 375       REAL tail
 376       INTEGER jstart, i, j, k
 377
 378 ! check whether flag given is legal
 379       IF ((FoldIn .NE. 0) .AND. (FoldIn .NE. 1)) THEN
 380          call wrf_error_fatal( 'inter3: ERROR <<<  Value for FOLDIN invalid. Must be 0 or 1' )
 381       ENDIF
 382
 383 ! do interpolation
 384
 385       jstart = 1
 386       DO i = 1, ng - 1
 387          yg(i) = 0.
 388          sum = 0.
 389          j = jstart
 390          IF (j .LE. n-1) THEN
 391    20      CONTINUE
 392              IF (x(j+1) .LT. xg(i)) THEN
 393                 jstart = j
 394                 j = j+1
 395                 IF (j .LE. n-1) GO TO 20
 396              ENDIF
 397    25      CONTINUE
 398            IF ((x(j) .LE. xg(i+1)) .AND. (j .LE. n-1)) THEN
 399               a1 = MAX(x(j),xg(i))
 400               a2 = MIN(x(j+1),xg(i+1))
 401               sum = sum + y(j) * (a2-a1)
 402               j = j+1
 403               GO TO 25
 404            ENDIF
 405            yg(i) = sum /(xg(i+1)-xg(i))
 406         ENDIF
 407       ENDDO
 408
 409
 410 ! if wanted, integrate data "overhang" and fold back into last bin
 411
 412       IF (FoldIn .EQ. 1) THEN
 413          j = j-1
 414          a1 = xg(ng)     ! upper limit of last interpolated bin
 415          a2 = x(j+1)     ! upper limit of last input bin considered
 416 !        do folding only if grids don't match up and there is more input
 417          IF ((a2 .GT. a1) .OR. (j+1 .LT. n)) THEN
 418            tail = y(j) * (a2-a1)/(x(j+1)-x(j))
 419            DO k = j+1, n-1
 420               tail = tail + y(k) * (x(k+1)-x(k))
 421            ENDDO
 422            yg(ng-1) = yg(ng-1) + tail
 423          ENDIF
 424       ENDIF
 425
 426       END SUBROUTINE inter4
 427
 428 !=============================================================================*
 429
 430       SUBROUTINE addpnt ( x, y, ld, n, xnew, ynew )
 431 !-----------------------------------------------------------------------------*
 432 !=  PURPOSE:                                                                 =*
 433 !=  Add a point <xnew,ynew> to a set of data pairs <x,y>.  x must be in      =*
 434 !=  ascending order                                                          =*
 435 !-----------------------------------------------------------------------------*
 436 !=  PARAMETERS:                                                              =*
 437 !=  X    - REAL vector of length LD, x-coordinates                       (IO)=*
 438 !=  Y    - REAL vector of length LD, y-values                            (IO)=*
 439 !=  LD   - INTEGER, dimension of X, Y exactly as declared in the calling  (I)=*
 440 !=         program                                                           =*
 441 !=  N    - INTEGER, number of elements in X, Y.  On entry, it must be:   (IO)=*
 442 !=         N < LD.  On exit, N is incremented by 1.                          =*
 443 !=  XNEW - REAL, x-coordinate at which point is to be added               (I)=*
 444 !=  YNEW - REAL, y-value of point to be added                             (I)=*
 445 !-----------------------------------------------------------------------------*
 446
 447       IMPLICIT NONE
 448
 449 ! calling parameters
 450
 451       INTEGER, intent(in) :: ld
 452       INTEGER, intent(inout) :: n
 453       REAL, intent(inout)    :: x(ld), y(ld)
 454       REAL, intent(in) :: xnew, ynew
 455
 456 ! local variables
 457
 458       INTEGER insert
 459       INTEGER i
 460       CHARACTER(len=256) :: emsg
 461
 462 ! check n<ld to make sure x will hold another point
 463       IF (n .GE. ld) THEN
 464          call wrf_error_fatal( 'addpnt: ERROR <<<  Cannot expand array All elements used.' )
 465       ENDIF
 466
 467       insert = 1
 468       i = 2
 469
 470 ! check whether x is already sorted.
 471 ! also, use this loop to find the point at which xnew needs to be inserted
 472 ! into vector x, if x is sorted.
 473
 474  10   CONTINUE
 475       IF (i .LT. n) THEN
 476         IF (x(i) .LT. x(i-1)) THEN
 477            call wrf_error_fatal( 'addpnt: ERROR <<<  x-data must be in ascending order!' )
 478         ELSE
 479            IF (xnew .GT. x(i-1)) insert = i
 480         ENDIF
 481         i = i+1
 482         GOTO 10
 483       ENDIF
 484
 485 ! if <xnew,ynew> needs to be appended at the end, just do so,
 486 ! otherwise, insert <xnew,ynew> at position INSERT
 487
 488       IF ( xnew .GT. x(n) ) THEN
 489          x(n+1) = xnew
 490          y(n+1) = ynew
 491       ELSE
 492 ! shift all existing points one index up
 493          DO i = n, insert, -1
 494            x(i+1) = x(i)
 495            y(i+1) = y(i)
 496          ENDDO
 497 ! insert new point
 498          x(insert) = xnew
 499          y(insert) = ynew
 500       ENDIF
 501
 502 ! increase total number of elements in x, y
 503
 504       n = n+1
 505
 506       END SUBROUTINE addpnt