text:
!dis
!dis    Open Source License/Disclaimer, Forecast Systems Laboratory
!dis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
!dis
!dis    This software is distributed under the Open Source Definition,
!dis    which may be found at http://www.opensource.org/osd.html.
!dis
!dis    In particular, redistribution and use in source and binary forms,
!dis    with or without modification, are permitted provided that the
!dis    following conditions are met:
!dis
!dis    - Redistributions of source code must retain this notice, this
!dis    list of conditions and the following disclaimer.
!dis
!dis    - Redistributions in binary form must provide access to this
!dis    notice, this list of conditions and the following disclaimer, and
!dis    the underlying source code.
!dis
!dis    - All modifications to this software must be clearly documented,
!dis    and are solely the responsibility of the agent making the
!dis    modifications.
!dis
!dis    - If significant modifications or enhancements are made to this
!dis    software, the FSL Software Policy Manager
!dis    (softwaremgr@fsl.noaa.gov) should be notified.
!dis
!dis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
!dis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
!dis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
!dis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
!dis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
!dis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
!dis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
!dis
!dis

MODULE horiz_interp

  ! Contains routines for doing horizontal interpolation
  IMPLICIT NONE
  PRIVATE

  ! Interp method codes for standard interpolation
  INTEGER, PUBLIC, PARAMETER          :: METHOD_NEAREST = 0
  INTEGER, PUBLIC, PARAMETER          :: METHOD_LINEAR  = 1
  INTEGER, PUBLIC, PARAMETER          :: METHOD_HIGHER = 2

  ! Data type codes for masked interpolation
  INTEGER, PUBLIC, PARAMETER          :: FLAG_FIELD = 0
  INTEGER, PUBLIC, PARAMETER          :: CATEGORY_FIELD = 1
  INTEGER, PUBLIC, PARAMETER          :: VALUE_FIELD = 2
  PUBLIC interpolate_standard
  PUBLIC interpolate_masked_val

CONTAINS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE interpolate_standard(nx_in, ny_in, data_in, &
                                  nx_out, ny_out, xloc, yloc, &
                                  method, data_out)
    !
    ! Performs standard interpolation (i.e., non-masked fields)
    IMPLICIT NONE

    ! Dimensions and grid we are interpolating from
    INTEGER, INTENT(IN)        :: nx_in
    INTEGER, INTENT(IN)        :: ny_in
    REAL, INTENT(IN)           :: data_in(nx_in,ny_in)

    ! Dimensions of output grid
    INTEGER, INTENT(IN)        :: nx_out
    INTEGER, INTENT(IN)        :: ny_out

    ! Pre-computed real i/j indices of input grid
    ! for each point in output grid
    REAL, INTENT(IN)           :: xloc(nx_out,ny_out)
    REAL, INTENT(IN)           :: yloc(nx_out,ny_out)

    ! Code indicating interp method (see module top)
    INTEGER, INTENT(IN)        :: method

    ! Output interpolated data
    REAL, INTENT(OUT)          :: data_out(nx_out,ny_out)

    ! Locals
    INTEGER                    :: i
    INTEGER                    :: j
    INTEGER                    :: ilo, jlo
    REAL                       :: dlon, dlat
    REAL                       :: deltalat, deltalon

    SELECT CASE(method)
      CASE(METHOD_NEAREST)
        DO j = 1, ny_out
          DO i = 1, nx_out
            ilo = NINT(xloc(i,j))
            jlo = NINT(yloc(i,j))
            data_out(i,j) = data_in( ilo, jlo )
          ENDDO
        ENDDO

      CASE(METHOD_LINEAR)
        DO j = 1, ny_out
          DO i = 1, nx_out
            ilo = MIN(FLOOR(xloc(i,j)), nx_in-1)
            jlo = MIN(FLOOR(yloc(i,j)), ny_in-1)
            dlon = xloc(i,j) - ilo
            dlat = yloc(i,j) - jlo
            deltalat = 1.
            deltalon = 1.
            data_out(i,j) = ((deltalon-dlon)*((deltalat-dlat)*data_in(ilo,jlo) &
                          + dlat*data_in(ilo,jlo+1))   &
                          + dlon*((deltalat-dlat)*data_in(ilo+1,jlo)     &
                          + dlat*data_in(ilo+1,jlo+1))) &
                          / ( deltalat * deltalon )
          ENDDO
        ENDDO
      CASE(METHOD_HIGHER)
        DO j = 1, ny_out
          DO i = 1, nx_out
            data_out(i,j) = BINT( xloc(i,j), yloc(i,j), data_in, nx_in, ny_in, 0)
          ENDDO
        ENDDO
      CASE DEFAULT
        PRINT '(A,I2)', 'Uknown interpolation method code: ', method
        STOP 'INTERPOLATE_STANDARD'
    END SELECT
    RETURN
  END SUBROUTINE interpolate_standard
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    SUBROUTINE interpolate_masked_val(nx_src, ny_src, lmask_src, data_src, &
                                 nx_out, ny_out, lmask_out, data_out, &
                                 isrcr, jsrcr, make_srcmask, &
                                 min_val, max_val, def_val, val_mask, &
                                 method)

    IMPLICIT none

    INTEGER, INTENT(IN)           :: nx_src
    INTEGER, INTENT(IN)           :: ny_src
    REAL, INTENT(INOUT)           :: lmask_src(nx_src,ny_src)
    REAL, INTENT(IN)              :: data_src(nx_src,ny_src)
    INTEGER, INTENT(IN)           :: nx_out, ny_out
    REAL, INTENT(IN)              :: lmask_out(nx_out,ny_out)
    REAL, INTENT(INOUT)             :: data_out(nx_out,ny_out)
    REAL, INTENT(IN)              :: isrcr(nx_out,ny_out)
    REAL, INTENT(IN)              :: jsrcr(nx_out,ny_out)
    LOGICAL, INTENT(IN)           :: make_srcmask
    REAL, INTENT(IN)              :: min_val, max_val, def_val, val_mask
    INTEGER, INTENT(IN)           :: method

    INTEGER                       :: i,j,k,l,ilo,jlo, ii, jj, kk, ll
    INTEGER                       :: search_rad
    INTEGER                       :: search_rad_max
    REAL, PARAMETER               :: bad_point = -99999.
    LOGICAL                       :: bad_points_flag
    LOGICAL                       :: fixed_bad
    INTEGER                       :: total_points
    INTEGER                       :: points_val_mask
    INTEGER                       :: points_skipped
    INTEGER                       :: num_to_fix
    INTEGER                       :: fixed_with_out
    INTEGER                       :: fixed_with_src
    INTEGER                       :: fixed_with_def
    INTEGER                       :: close_pts
    INTEGER                       :: ic,jc,iic,jjc
    REAL                          :: deltagx, deltagy
    REAL                          :: dix, djy
    REAL                          :: data_close(4)
    REAL                          :: distance(4)
    REAL                          :: distx,disty
    REAL                          :: sum_dist
    REAL                          :: a,b,c,d,e,f,g,h
    REAL                          :: stl(4,4)
    REAL                          :: valb
    REAL                          :: srcmask(nx_src,ny_src)

    ! Can we use the source data land mask that was input, or do we need to make it
    ! from the min_val/max_val arguments?

    IF (make_srcmask) THEN
      PRINT '(A)', 'INTERPOLATE_MASKED_VAL: Making source data land mask...'

      ! This code has the flexibility to handled either water data points (lmask =0)
      ! or land points (lmask = 1). So if we are making the landmask using valid
      ! range, but we do not know anything about the variable other than the
      ! valid mask value, we need to figure out which points are 0 and which should be
      ! 1.  To do this, initialize the array to the invalid value, which we determine
      ! from the valid value.
      IF (val_mask .EQ. 1) THEN
        lmask_src(:,:) = 0
      ELSE
        lmask_src(:,:) = 1
      ENDIF

      ! Now figure out where to set the mask using a WHERE statement
      WHERE((data_src .GE. min_val).AND.(data_src .LE. max_val)) lmask_src = val_mask
    ELSE
      PRINT '(A)', 'INTERPOLATE_MASKED: Using source landmask field.'
    ENDIF

    bad_points_flag = .false.

    ! Initialize counters
    total_points = 0
    points_val_mask = 0
    points_skipped = 0
    num_to_fix = 0
    fixed_with_out = 0
    fixed_with_src = 0
    fixed_with_def = 0

    ! Select interpolation method.  Putting the case statement out here
    ! increases the amount of replicated code but should be more efficient
    ! than checking this condition at every grid point.

    SELECT CASE(method)

      CASE(METHOD_NEAREST)
        ! Use nearest neigbor
        PRINT '(A)', 'Masked interpolation using nearest neighbor value...'
        out_j_loop_1: DO j = 1, ny_out
          out_i_loop_1: DO i = 1, nx_out

            total_points = total_points + 1
            ! We only need to process this point if the lmask_out is equal
            ! to val_mask.  For example, if one is processing soil parameters,
            ! which are only valid at land points (lmask = 1), then the user
            ! passes in 1. for the val_mask.  Any output point that is water
            ! (lmask = 0) will then be skipped.  During this loop, if we have
            ! points that cannot be properly assigned, we will mark them as bad
            ! and set the bad_points_lag.

            IF (lmask_out(i,j) .EQ. val_mask) THEN
              ! Process this point
              points_val_mask = points_val_mask + 1
              ilo = NINT(isrcr(i,j))
              jlo = NINT(jsrcr(i,j))

              ! See if this point can be used
              IF (lmask_src(ilo,jlo).EQ. lmask_out(i,j)) THEN
                data_out(i,j) = data_src(ilo,jlo)
              ELSE
                data_out(i,j) = bad_point
                num_to_fix = num_to_fix + 1
                bad_points_flag = .true.
              ENDIF
            ELSE
              ! The output grid does not require a value for this point
              ! But do not zero out in case this is a field begin
              ! done twice (once for water and once for land, e.g.
              ! SKINTEMP
              !data_out(i,j) = 0.
              points_skipped = points_skipped + 1
            ENDIF
          ENDDO out_i_loop_1
        ENDDO out_j_loop_1

      CASE (METHOD_LINEAR)
        ! Use a 4-point interpolation
        PRINT '(A)', 'Masked interpolation using 4-pt linear interpolation...'
        deltagx = 1.
        deltagy = 1.
        out_j_loop_2: DO j = 1, ny_out
          out_i_loop_2: DO i = 1, nx_out

            total_points = total_points + 1
            ! We only need to process this point if the lmask_out is equal
            ! to val_mask.  For example, if one is processing soil parameters,
            ! which are only valid at land points (lmask = 1), then the user
            ! passes in 1. for the val_mask.  Any output point that is water
            ! (lmask = 0) will then be skipped.  During this loop, if we have
            ! points that cannot be properly assigned, we will mark them as bad
            ! and set the bad_points_lag.

            IF (lmask_out(i,j) .EQ. val_mask) THEN
              ! Process this point
              points_val_mask = points_val_mask + 1
              ilo = MIN(FLOOR(isrcr(i,j)),nx_src-1)
              jlo = MIN(FLOOR(jsrcr(i,j)),ny_src-1)
              dix = isrcr(i,j) - FLOAT(ilo)
              djy = jsrcr(i,j) - FLOAT(jlo)

              ! Loop around the four surrounding points
              ! and count up the number of points we can use based
              ! on common mask value
              close_pts = 0
              sum_dist = 0.
              outer_four_j: DO jc = 0,1
                outer_four_i: DO ic = 0,1
                  iic = ilo + ic
                  jjc = jlo + jc
                  IF (lmask_src(iic,jjc).EQ. lmask_out(i,j)) THEN
                    close_pts = close_pts + 1
                    data_close(close_pts) = data_src(iic,jjc)

                    ! Compute distance to this valid point
                    ! in grid units and add to sum of distances (actually,
                    ! we are doing a weight, which is inversely proportional
                    ! to distance)
                    IF (ic .EQ. 0) THEN
                      distx = deltagx - dix
                    ELSE
                      distx =  dix
                    ENDIF
                    IF (jc .EQ. 0) THEN
                      disty = deltagy - djy
                    ELSE
                      disty = djy
                    ENDIF
                    distance(close_pts) = SQRT(distx**2+disty**2)
                    sum_dist = sum_dist + distance(close_pts)
                  ENDIF
                ENDDO outer_four_i
              ENDDO outer_four_j

              ! Did we find at least one point in the surrounding four
              ! that was usable?

              IF (close_pts .GT. 0) THEN

                ! If we have all four points, then do bilinear interpolation
                IF (close_pts .EQ. 4) THEN
                   data_out(i,j) = ((deltagx - dix)*((deltagy-djy)*data_close(1) &
                                 + djy*data_close(3)) &
                                 + dix*((deltagy-djy)*data_close(2) &
                                 + djy*data_close(4))) &
                                 / (deltagx * deltagy)
                ELSE IF ((close_pts .GT. 1).AND.(close_pts .LT. 4)) THEN

                  ! Simple distance-weighted average by computing
                  ! the sum of all distances to each point and using
                  ! each individual distance divided by the total
                  ! distance as the weighting

                  data_out(i,j) = 0.
                  DO k = 1, close_pts
                    data_out(i,j) = data_out(i,j) + &
                                    (distance(k)/sum_dist) * data_close(k)
                  ENDDO
                ELSE
                  ! Set output value = to one point we found
                  data_out(i,j) = data_close(1)
                ENDIF
              ELSE
                bad_points_flag = .true.
                data_out(i,j) = bad_point
                num_to_fix = num_to_fix + 1
              ENDIF
            ELSE
              ! The output grid does not require a value for this point
              ! But do not zero out in case this is a field begin                             ! done twice (once for water and once for land, e.g.                            ! SKINTEMP
              !data_out(i,j) = 0.
              points_skipped = points_skipped + 1
            ENDIF
          ENDDO out_i_loop_2
        ENDDO out_j_loop_2

      CASE (METHOD_HIGHER)
        ! 16-point interpolation
        out_j_loop_3: DO j = 1, ny_out
          out_i_loop_3: DO i = 1, nx_out

            total_points = total_points + 1
            ! We only need to process this point if the lmask_out is equal
            ! to val_mask.  For example, if one is processing soil parameters,
            ! which are only valid at land points (lmask = 1), then the user
            ! passes in 1. for the val_mask.  Any output point that is water
            ! (lmask = 0) will then be skipped.  During this loop, if we have
            ! points that cannot be properly assigned, we will mark them as bad
            ! and set the bad_points_lag.

            IF (lmask_out(i,j) .EQ. val_mask) THEN
              ! Process this point
              points_val_mask = points_val_mask + 1

              ! Do a 4x4 loop around in the input data around the output
              ! point to get our 16 points of influence.  Only use those
              ! that are appropriately masked
              valb = 0.
              close_pts = 0
              ilo = INT(isrcr(i,j)+0.00001)
              jlo = INT(jsrcr(i,j)+0.00001)
              dix = isrcr(i,j) - ilo
              djy = jsrcr(i,j) - jlo
              IF ( (ABS(dix).GT.0.0001).OR.(ABS(djy).GT.0.0001) ) THEN
                ! Do the interpolation loop
                stl(:,:) = 0.
                loop_16_1: DO k = 1,4
                  kk = ilo + k - 2
                  IF ((kk .LT. 1).OR.(kk .GT. nx_src)) CYCLE loop_16_1
                  loop_16_2: DO l = 1, 4
                    ll = jlo + l - 2
                    IF ((ll .LT. 1).OR.(ll .GT. ny_src)) CYCLE loop_16_2

                    ! Check land mask at this source point
                    IF (lmask_src(kk,ll).NE. val_mask) CYCLE loop_16_2
                    ! If we are here, then mask tests passed
                    stl(k,l) = data_src(kk,ll)
                    IF ( (stl(k,l) .EQ. 0.).AND.(min_val.LE.0.).AND. &
                                                (max_val.GE.0.) ) THEN
                      stl = 1.E-5
                    ENDIF
                    close_pts = close_pts + 1
                  ENDDO loop_16_2
                ENDDO loop_16_1

                ! Did we find any valid points?

                IF ( (close_pts .GT. 0).AND. ( &
                  (stl(2,2).GT.0.).AND.(stl(2,3).GT.0.).AND. &
                  (stl(3,2).GT.0.).AND.(stl(3,3).GT.0.)  ) ) THEN
                  a = oned(dix,stl(1,1),stl(2,1),stl(3,1),stl(4,1))
                  b = oned(dix,stl(1,2),stl(2,2),stl(3,2),stl(4,2))
                  c = oned(dix,stl(1,3),stl(2,3),stl(3,3),stl(4,3))
                  d = oned(dix,stl(1,4),stl(2,4),stl(3,4),stl(4,4))
                  valb = oned(djy,a,b,c,d)
                  IF (close_pts .NE. 16) THEN
                    e = oned(djy,stl(1,1),stl(1,2),stl(1,3),stl(1,4))
                    f = oned(djy,stl(2,1),stl(2,2),stl(2,3),stl(2,4))
                    g = oned(djy,stl(3,1),stl(3,2),stl(3,3),stl(3,4))
                    h = oned(djy,stl(4,1),stl(4,2),stl(4,3),stl(4,4))
                    valb = (valb+oned(dix,e,f,g,h)) * 0.5
                  ENDIF
                  data_out(i,j) = valb

                ELSE
                  bad_points_flag = .true.
                  data_out(i,j) = bad_point
                  num_to_fix = num_to_fix + 1
                ENDIF
              ELSE
                ! We are right on a source point, so try to use it
                IF (lmask_src(ilo,jlo).EQ.val_mask) THEN
                  data_out(i,j) = data_src(ilo,jlo)
                ELSE
                  bad_points_flag = .true.
                  data_out(i,j) = bad_point
                  num_to_fix = num_to_fix + 1
                ENDIF
              ENDIF
            ELSE
              ! The output grid does not require a value for this point
              !data_out(i,j) = 0.
              points_skipped = points_skipped + 1
            ENDIF
          ENDDO out_i_loop_3
        ENDDO out_j_loop_3

    END SELECT

     ! Do we need to correct bad points?

    IF (bad_points_flag) THEN

      search_rad_max = 10
      fix_bad_j: DO j = 1, ny_out
        fix_bad_i: DO i = 1, nx_out

          IF (data_out(i,j).NE. bad_point) CYCLE fix_bad_i

          ! First, search for nearest non-bad point in the output domain
          ! which is usually higher resolution.
          fixed_bad = .false.
          search_out_loop: DO search_rad = 1, search_rad_max
            search_out_j: DO ll = -(search_rad-1), (search_rad-1),1
              jj = j + ll
              IF ((jj .LT. 1).OR.(jj .GT. ny_out)) CYCLE search_out_j
              search_out_i: DO kk = -(search_rad), search_rad, 1
                 ii = i + kk
                 IF ((ii .LT. 1).OR.(ii .GT. nx_out)) CYCLE search_out_j
                 IF ((data_out(ii,jj).NE.bad_point).AND. &
                    (lmask_out(ii,jj) .EQ. val_mask) ) THEN
                  data_out(i,j) = data_out(ii,jj)
                  fixed_bad = .true.
                  fixed_with_out = fixed_with_out + 1
                  EXIT search_out_loop
                ENDIF
              ENDDO search_out_i
            ENDDO search_out_j
          ENDDO search_out_loop

          ! Did we fix the point?  If not, then do same search on src data.
          IF (.NOT. fixed_bad) THEN
            search_rad_max = 10
            search_src_loop: DO search_rad = 1, search_rad_max
              search_src_j: DO ll = -(search_rad-1), (search_rad-1),1
                jj = NINT(jsrcr(i,j)) + ll
                IF ((jj .LT. 1).OR.(jj .GT. ny_src)) CYCLE search_src_j
                search_src_i: DO kk = -(search_rad), search_rad, 1
                   ii = NINT(isrcr(i,j)) + kk
                   IF ((ii .LT. 1).OR.(ii .GT. nx_src)) CYCLE search_src_j
                   IF (lmask_src(ii,jj).EQ.val_mask) THEN
                     data_out(i,j) = data_src(ii,jj)
                     fixed_bad = .true.
                     fixed_with_src = fixed_with_src + 1
                     EXIT search_src_loop
                   ENDIF
                ENDDO search_src_i
              ENDDO search_src_j
            ENDDO search_src_loop
          ENDIF
          ! Now is the point fixed?  If not, we have to use a default value.
          IF (.NOT.fixed_bad) THEN
            fixed_with_def = fixed_with_def + 1
            data_out(i,j) = def_val
            PRINT '(A,F10.3,A,2I5)', 'INTERPOLATE_MASKED: Bogus value of ', def_val, &
                ' used at point ', i, j
          ENDIF

        ENDDO fix_bad_i
      ENDDO fix_bad_j
    ENDIF
    PRINT '(A)',     '----------------------------------------'
    PRINT '(A)',     'MASKED INTERPOLATION SUMMARY: '
    PRINT '(A,I10)', '  TOTAL POINTS IN GRID:       ', total_points
    PRINT '(A,I10)', '  POINTS NEEDING VALUES:      ', points_val_mask
    PRINT '(A,I10)', '  POINTS NOT REQUIRED:        ', points_skipped
    PRINT '(A,I10)', '  POINTS NEEDING FIX:         ', num_to_fix
    PRINT '(A,I10)', '  POINTS FIXED WITH OUT GRID: ', fixed_with_out
    PRINT '(A,I10)', '  POINTS FIXED WITH SRC GRID: ', fixed_with_src
    PRINT '(A,I10)', '  POINTS FIXED WITH DEF VAL:  ', fixed_with_def
    PRINT '(A)',     '----------------------------------------'
    RETURN
  END SUBROUTINE interpolate_masked_val

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   FUNCTION bint(xx,yy,list,iii,jjj,ibint)

      !  16 Point interpolation

      IMPLICIT NONE

      REAL :: xx , yy
      INTEGER :: ibint , iii, jjj
      REAL list(iii,jjj),stl(4,4)

      INTEGER :: ib , jb, n , i , j , k , kk , l , ll
      REAL :: bint , x , y , a , b , c , d , e , f , g , h
      ib=iii-ibint
      jb=jjj-ibint
      bint = 0.0
      n = 0
      i = INT(xx+0.00001)
      j = INT(yy+0.00001)
      x = xx - i
      y = yy-j

      IF ( ( ABS(x).GT.0.0001 ) .OR. ( abs(y).gt.0.0001 ) ) THEN
         stl(:,:)=1.E-10
         loop_1 : DO k = 1,4
            kk = i + k - 2
            IF ( ( kk .LT. 1) .OR. ( kk .GT. ib ) ) THEN
               CYCLE loop_1
            END IF
            loop_2 : DO l = 1,4
               stl(k,l) = 0.
               ll = j + l - 2
               IF ( ( ll .GT. jb ) .OR. ( ll .LT. 1 ) ) THEN
                  CYCLE loop_2
               END IF
               stl(k,l) = list(kk,ll)
               n = n + 1
               IF ( stl(k,l) .EQ. 0. ) THEN
                  stl(k,l) = 1.E-20
               END IF
            END DO loop_2
         END DO loop_1
         a = oned(x,stl(1,1),stl(2,1),stl(3,1),stl(4,1))
         b = oned(x,stl(1,2),stl(2,2),stl(3,2),stl(4,2))
         c = oned(x,stl(1,3),stl(2,3),stl(3,3),stl(4,3))
         d = oned(x,stl(1,4),stl(2,4),stl(3,4),stl(4,4))
         bint = oned(y,a,b,c,d)

         IF(n.NE.16) THEN
            e = oned(y,stl(1,1),stl(1,2),stl(1,3),stl(1,4))
            f = oned(y,stl(2,1),stl(2,2),stl(2,3),stl(2,4))
            g = oned(y,stl(3,1),stl(3,2),stl(3,3),stl(3,4))
            h = oned(y,stl(4,1),stl(4,2),stl(4,3),stl(4,4))
            bint = (bint+oned(x,e,f,g,h)) * 0.5
         END IF

      ELSE
         bint = list(i,j)
      END IF
   END FUNCTION bint

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   FUNCTION oned(x,a,b,c,d)

      IMPLICIT NONE

      REAL :: x,a,b,c,d,oned

      oned = 0.

      IF      ( x .EQ. 0. ) THEN
         oned = b
      ELSE IF ( x .EQ. 1. ) THEN
         oned = c
      END IF

      IF(b*c.NE.0.) THEN
         IF ( a*d .EQ. 0. ) THEN
            IF      ( ( a .EQ. 0 ) .AND. ( d .EQ. 0 ) ) THEN
               oned = b*(1.0-x)+c*x
            ELSE IF ( a .NE. 0. ) THEN
               oned = b+x*(0.5*(c-a)+x*(0.5*(c+a)-b))
            ELSE IF ( d .NE. 0. ) THEN
               oned = c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c))
            END IF
         ELSE
            oned = (1.0-x)*(b+x*(0.5*(c-a)+x*(0.5*(c+a)-b)))+ &
                   x*(c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c)))
         END IF
      END IF

   END FUNCTION oned
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE horiz_interp

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
      subroutine gauher2(x,w,n)
      implicit none
      integer n,maxit
      real*8 w(n),x(n)
      real*8 eps,pim4
      parameter(eps=3.d-14,pim4=.7511255444649425d0,maxit=10)
      integer i,its,j,m
      real*8 p1,p2,p3,pp,z,z1

      m=(n+1)/2
      z=0.d0
      do i=1,n
         w(i)=0.0
         x(i)=0.0
      end do

      do 13 i=1,m
        if(i.eq.1)then
          z=sqrt(float(2*n+1))-1.85575*(2*n+1)**(-.16667)
        else if(i.eq.2)then
          z=z-1.14*n**.426/z
        else if (i.eq.3)then
          z=1.86*z-.86*x(1)
        else if (i.eq.4)then
          z=1.91*z-.91*x(2)
        else
          z=2.*z-x(i-2)
        endif
        do 12 its=1,MAXIT
          p1=PIM4
          p2=0.d0
          do 11 j=1,n
            p3=p2
            p2=p1
            p1=z*sqrt(2.d0/j)*p2-sqrt(dble(j-1)/dble(j))*p3
11        continue
          pp=sqrt(2.d0*n)*p2
          z1=z
          z=z1-p1/pp
          if(abs(z-z1).le.EPS)goto 1
12      continue
        call rexit("too many iterations in gauher")
1       x(i)=z
        x(n+1-i)=-z
        w(i)=2.d0/(pp*pp)
        w(n+1-i)=w(i)
13    continue
      return
      END

ext: f

max_stars_repo_name: kdevick/densitymediation

max_issues_repo_name: kdevick/densitymediation

max_forks_repo_name: kdevick/densitymediation


text:
*deck argos
      program argos
      implicit logical(a-z)
c
c  allocate space and call the driver routine.
c
c  ifirst = first usable space in the the real*8 work array a(*).
c  lcore  = length of workspace array in working precision.
c  mem1   = additional machine-dependent memory allocation variable.
c  a(*)   = workspace array. a(ifirst : ifirst+lcore-1) is useable.
c
      integer ifirst, lcore
c
c     # mem1 and a(*) should be declared below within mdc blocks.
c
c     # local...
c     # lcored = default value for lcore.
c     # ierr   = error return code.
c
      integer lcored, ierr
c
c     # bummer error types.
      integer   wrnerr,  nfterr,  faterr
      parameter(wrnerr=0,nfterr=1,faterr=2)
c
c     # "call pbginf" should be the first executable statement.  this
c     # sets up any communications necessary for parallel execution,
c     # and does nothing in the sequential case.
c
*mdc*if unicos
*c
*c     # nonstandard unicos memory allocation.
*c     # mem1 is a pointer to the beginning of a(*).
*c
*      parameter (lcored = 1 800 000)
*      pointer (mem1, a)
*      real*8 a(2)
*c
*      call pbginf
*c
*      call getlcl( lcored, lcore, ierr )
*      if ( ierr .ne. 0 ) then
*         call bummer('argos: from getlcl: ierr=',ierr,faterr)
*      endif
*c
*      call hpalloc( mem1, lcore, ierr, 1 )
*      ifirst = 1
*mdc*elseif alliant
*c
*c     # fortran 90 style core allocation.
*c
*      integer mem1
*      parameter (lcored = 1 800 000)
*      real*8 a(:)
*c
*      call pbginf
*c
*      call getlcl( lcored, lcore, ierr )
*      if ( ierr .ne. 0 ) then
*         call bummer('argos: from getlcl: ierr=',ierr,faterr)
*      endif
*c
*      allocate( a(lcore) )
*c
*      mem1   = 0
*      ifirst = 1
*mdc*elseif sun iris convex
c
c     # nonstandard sun bsd unix memory allocation.
c     # mem1 is a pointer to the beginning of a(*).
c
      integer mem1
      parameter (lcored = 1 800 000)
      pointer (mem1, a)
      real*8 a(2)
c
c     integer  malloc
c     external malloc
c
      call pbginf
c
      call getlcl( lcored, lcore, ierr )
      if ( ierr .ne. 0 ) then
         call bummer('argos: from getlcl: ierr=',ierr,faterr)
      endif
c
c     #  malloc() argument is assumed to be in bytes.
      mem1 = malloc( 8*lcore )
      if ( mem1 .eq. 0 ) then
         call bummer('argos: from malloc: mem1=',mem1,faterr)
      endif
      ifirst = 1
*mdc*elseif (fps .and. obsolete)
*c
*c     # this code is for machines that support malloc(),
*c     # but do not support pointer syntax.
*c
*      integer mem1
*      parameter (lcored = 1 800 000)
*c
*      integer  malloc
*      external malloc
*c
*      call pbginf
*c
*      call getlcl( lcored, lcore, ierr )
*      if ( ierr .ne. 0 ) then
*         call bummer('argos: from getlcl: ierr=',ierr,faterr)
*      endif
*c
*c     #  malloc() argument is assumed to be in bytes.
*      mem1 = malloc( 8*lcore )
*      if ( mem1 .eq. 0 ) then
*         call bummer('argos: from malloc: mem1=',mem1,faterr)
*      endif
*      ifirst = 1
*mdc*elseif unix titan
*c
*c     # generic unix memory allocation using falloc() interface.
*c     # see also $COLUMBUS/special/unix/falloc.c
*c
*      integer mem1, offset
*      parameter (lcored =  1 800 000)
*      real*8 a(1)
*c
*      call pbginf
*c
*      call getlcl( lcored, lcore, ierr )
*      if ( ierr .ne. 0 ) then
*         call bummer('argos: from getlcl: ierr=',ierr,faterr)
*      endif
*c
*c     # on return from falloc(), a(offset+1:offset+lcore) is useable.
*      call falloc( lcore, 8, 0, a, mem1, offset )
*      if ( mem1 .eq. 0 ) then
*         call bummer('argos: from falloc: mem1=',mem1,faterr)
*      endif
*      ifirst = offset + 1
*mdc*else
*c
*c     # standard fortran workspace allocation.
*c     # mem1 is not used.  use blank common to reduce executable size.
*c
*      integer mem1
*      parameter ( lcored = 4 000 000 )
*      real*8    a
*      common // a( lcored )
*c
*      call pbginf
*c
*      call getlcl( lcored, lcore, ierr )
*      if ( ierr .ne. 0 ) then
*         call bummer('argos: from getlcl: ierr=',ierr,faterr)
*      elseif ( lcore .gt. lcored ) then
*c        # make sure lcore is consistent.  this check should be
*c        # redundant with logic within getlcl(), but is included
*c        # here just to be safe.
*         call bummer('argos: (lcore-lcored)=',
*     &    (lcore-lcored), faterr )
*      endif
*c
*      mem1   = 0
*      ifirst = 1
*mdc*endif
c
c     # lcore, mem1, and ifirst should now be defined.
c     # the values of lcore, mem1, and ifirst are passed to
c     # driver() to be printed.
c     # since these arguments are passed by expression, they
c     # must not be modified within driver().
c
*mdc*if (fps .and. obsolete)
*c
*c      # for machines that support malloc()-type calls, but do not
*c      # support pointer syntax,  the workaround is to pass the
*c      # address of a(*) by value.
*c
*      call driver( %val(mem1), (lcore), (mem1), (ifirst) )
*c
*mdc*else
c
      call driver( a(ifirst), (lcore), (mem1), (ifirst) )
c
*mdc*endif
c
c     # return from driver() implies successful execution.
c
c     # message-passing cleanup: stub if not in parallel
      call pend
c
      stop 'end of argos'
      end

max_stars_repo_name: MOLFDIR/MOLFDIR

max_issues_repo_name: yingxingcheng/MOLFDIR

max_forks_repo_name: yingxingcheng/MOLFDIR


text:
#ifdef TIMERS
#define NBTIMER(a) a = dnekclock()
#define STIMER(a) a = dnekclock_sync()
#define ACCUMTIMER(b,a) b = b + (dnekclock()- a )
#else
#define NBTIMER(a)
#define STIMER(a)
#define ACCUMTIMER(a,b)
#endif

c-----------------------------------------------------------------------
      subroutine cg(x,f,g,c,r,w,p,z,n,niter,flop_cg)
      include 'SIZE'
      include 'TIMER'

c     Solve Ax=f where A is SPD and is invoked by ax()
c
c     Output:  x - vector of length n
c
c     Input:   f - vector of length n
c     Input:   g - geometric factors for SEM operator
c     Input:   c - inverse of the counting matrix
c
c     Work arrays:   r,w,p,z  - vectors of length n
c
c     User-provided ax(w,z,n) returns  w := Az,
c
c     User-provided solveM(z,r,n) ) returns  z := M^-1 r,
c
      parameter (lt=lx1*ly1*lz1*lelt)
c     real ur(lt),us(lt),ut(lt)

c     parameter (lxyz=lx1*ly1*lz1)
c     real ur(lxyz),us(lxyz),ut(lxyz),wk(lxyz)

      real x(n),f(n),r(n),w(n),p(n),z(n),g(1),c(n)
      real rnorminit, fbeta, fpap, falpha, frnorm

      integer thread, numth, find, lind, fel, lel
      integer omp_get_thread_num, omp_get_num_threads
      integer fiter, tmt

      pap = 0.0

c     set machine tolerances
      one = 1.
      eps = 1.e-20
      if (one+eps .eq. one) eps = 1.e-14
      if (one+eps .eq. one) eps = 1.e-7

      rtz1=1.0
      miter = niter

c$OMP PARALLEL DEFAULT(shared) PRIVATE(thread,numth,find,lind,iter,
c$OMP&  fel,lel,rtz2,beta,alpha,alphm,rlim2,rtr0)

      thread = 0
      numth = 1
#ifdef _OPENMP
      thread = omp_get_thread_num()
      numth = omp_get_num_threads()
#endif
      tmt = thread + 1

      if (numth < nelt) then
        fel = (thread*nelt)/numth + 1
        lel = ((thread+1)*nelt)/numth
      else
        if (thread < nelt) then
          fel = thread + 1
          lel = fel
        else
          fel = nelt+1
          lel = nelt
        end if
      end if

      find = (fel-1) *(nx1*ny1*nz1)+1
      lind = lel * (nx1*ny1*nz1)

      NBTIMER(ttemp1)
      call rzeroi(x,n,find,lind)
      call copyi(r,f,n,find,lind)
      ACCUMTIMER(tzc(tmt), ttemp1)

      if (thread == 0) call mask (r)   ! Zero out Dirichlet conditions

      NBTIMER(ttemp1)
      call glsc3i(rnorminit,r,c,r,n,find,lind)
      ACCUMTIMER(tglsc3(tmt), ttemp1)

      do iter=1,miter
         NBTIMER(ttemp1)
         call solveMi(z,r,n,find,lind)    ! preconditioner here
         ACCUMTIMER(tzc(tmt), ttemp1)

         rtz2=rtz1                                                       ! OPS
         NBTIMER(ttemp1)
         call glsc3i(rtz1,r,c,z,n,find,lind)
         ACCUMTIMER(tglsc3(tmt), ttemp1)

         beta = rtz1/rtz2
         if (iter.eq.1) beta=0.0

         NBTIMER(ttemp1)
         call add2s1i(p,z,beta,n,find,lind)                              ! 2n
         ACCUMTIMER(tadd2sx(tmt), ttemp1)

         NBTIMER(ttemp1)
         call axi(w,p,g,n,fel,lel,find,lind)                             ! flopa
         ACCUMTIMER(tax(tmt), ttemp1)

         NBTIMER(ttemp1)
         call glsc3i(pap, w,c,p,n,find,lind)                             ! 3n
         ACCUMTIMER(tglsc3(tmt), ttemp1)

         alpha=rtz1/pap
         alphm=-alpha

         NBTIMER(ttemp1)
         call add2s2i(x,p,alpha,n,find,lind)                             ! 2n
         call add2s2i(r,w,alphm,n,find,lind)                             ! 2n
         ACCUMTIMER(tadd2sx(tmt), ttemp1)

         NBTIMER(ttemp1)
         call  glsc3i(rtr, r,c,r,n,find,lind)                            ! 3n
         ACCUMTIMER(tglsc3(tmt), ttemp1)

         if (iter.eq.1) rlim2 = rtr*eps**2
         if (iter.eq.1) rtr0  = rtr
         rnorm = sqrt(rtr)

      enddo

      if (thread == 0) then
        fiter = iter
        fbeta = beta
        falpha= alpha
        fpap  = pap
        frnorm = rnorm
      end if

c$OMP END PARALLEL

    6    format('cg:',i4,1p4e12.4)

      if (nid.eq.0) then
        write(6,6) 0,sqrt(rnorminit)
        write(6,6) fiter,frnorm,falpha,fbeta,fpap
      end if

      flop_cg = flop_cg + (fiter-1)*15.0*n + 3.0*n

      return
      end
c-----------------------------------------------------------------------
      subroutine solveM(z,r,n)
      real z(n),r(n)

      call copy(z,r,n)

      return
      end
c-----------------------------------------------------------------------
      subroutine axi(w,u,gxyz,n,fel,lel,find,lind) ! Matrix-vector product: w=A*u

      include 'SIZE'
      include 'TOTAL'
      include 'TIMER'

      parameter (lxyz=lx1*ly1*lz1)
      real w(nx1*ny1*nz1,nelt),u(nx1*ny1*nz1,nelt)
      real gxyz(2*ldim,nx1*ny1*nz1,nelt)
      parameter (lt=lx1*ly1*lz1*lelt)

      integer fel, lel, find, lind
      integer e,thread, tmt, omp_get_thread_num

      thread = 0
#ifdef _OPENMP
      thread = omp_get_thread_num()
#endif
      tmt = thread + 1

      NBTIMER(ttemp2)
      do e= fel, lel
         call ax_e( w(1,e),u(1,e),gxyz(1,1,e))
      enddo
      ACCUMTIMER(taxe(tmt),ttemp2)

      NBTIMER(ttemp2)
      call gs_op(gsh,w,1,1,0)  ! Gather-scatter operation  ! w   = QQ  w
      ACCUMTIMER(tgsop(tmt),ttemp2)
                                                           !            L
      NBTIMER(ttemp2)
      call add2s2i(w,u,.1,n,find,lind)
      ACCUMTIMER(tadd2s2(tmt),ttemp2)

      if (find == 1) then
        call mask(w)             ! Zero out Dirichlet conditions
        nxyz=nx1*ny1*nz1
        flop_a = flop_a + (19*nxyz+12*nx1*nxyz)*nelt
      end if

      return
      end
c-------------------------------------------------------------------------
      subroutine ax1(w,u,n)
      include 'SIZE'
      real w(n),u(n)
      real h2i

      h2i = (n+1)*(n+1)
      do i = 2,n-1
         w(i)=h2i*(2*u(i)-u(i-1)-u(i+1))
      enddo
      w(1)  = h2i*(2*u(1)-u(2  ))
      w(n)  = h2i*(2*u(n)-u(n-1))

      return
      end
c-------------------------------------------------------------------------
      subroutine ax_e(w,u,g) ! Local matrix-vector product

      include 'SIZE'
      include 'TOTAL'

      parameter (lxyz=lx1*ly1*lz1)
      real w(lxyz),u(lxyz),g(2*ldim,lxyz)

      real ur(nx1*ny1*nz1),us(nx1*ny1*nz1),ut(nx1*ny1*nz1)

      nxyz = nx1*ny1*nz1
      n    = nx1-1

      call local_grad3(ur,us,ut,u,n,dxm1,dxtm1)

      do i=1,nxyz
         wr = g(1,i)*ur(i) + g(2,i)*us(i) + g(3,i)*ut(i)
         ws = g(2,i)*ur(i) + g(4,i)*us(i) + g(5,i)*ut(i)
         wt = g(3,i)*ur(i) + g(5,i)*us(i) + g(6,i)*ut(i)
         ur(i) = wr
         us(i) = ws
         ut(i) = wt
      enddo

      call local_grad3_t(w,ur,us,ut,n,dxm1,dxtm1)

      return
      end
c-------------------------------------------------------------------------
      subroutine local_grad3(ur,us,ut,u,n,D,Dt)
c     Output: ur,us,ut         Input:u,n,D,Dt
      real ur(0:n,0:n,0:n),us(0:n,0:n,0:n),ut(0:n,0:n,0:n)
      real u (0:n,0:n,0:n)
      real D (0:n,0:n),Dt(0:n,0:n)
      integer e

      m1 = n+1
      m2 = m1*m1

      call mxm(D ,m1,u,m1,ur,m2)
      do k=0,n
         call mxm(u(0,0,k),m1,Dt,m1,us(0,0,k),m1)
      enddo
      call mxm(u,m2,Dt,m1,ut,m1)

      return
      end
c-----------------------------------------------------------------------
      subroutine local_grad3_t(u,ur,us,ut,N,D,Dt)
c     Output: ur,us,ut         Input:u,N,D,Dt
      real u (0:N,0:N,0:N)
      real ur(0:N,0:N,0:N),us(0:N,0:N,0:N),ut(0:N,0:N,0:N)
      real D (0:N,0:N),Dt(0:N,0:N)
      real w (0:N,0:N,0:N)
      integer e

      m1 = N+1
      m2 = m1*m1
      m3 = m1*m1*m1

      call mxm(Dt,m1,ur,m1,u,m2)

      do k=0,N
         call mxm(us(0,0,k),m1,D ,m1,w(0,0,k),m1)
      enddo
      call add2(u,w,m3)

      call mxm(ut,m2,D ,m1,w,m1)
      call add2(u,w,m3)

      return
      end
c-----------------------------------------------------------------------
      subroutine mask(w)   ! Zero out Dirichlet conditions
      include 'SIZE'
      real w(1)

      if (nid.eq.0) w(1) = 0.  ! suitable for solvability

      return
      end
c-----------------------------------------------------------------------

max_stars_repo_name: HPCToolkit/hpctest

max_issues_repo_name: HPCToolkit/hpctest

max_forks_repo_name: HPCToolkit/hpctest


text:
      logical function chkpra(inproc)
c-------------------------------------------------------------------------------
c chkproc checks to see if there if the user has specified an FK process in the
c procs list.
c   Outputs:
c      .TRUE.  if an FK process is specified.
c      .FALSE. if no FK process specified
c
c   Last Modified:
c      10/20/88
c-------------------------------------------------------------------------------
c
      character*(*) inproc
c
c.. The following are needed to check the user processes (set in GETPRO)
      parameter ( NAMESS  = 47)
      integer pnum                                                      ! The number of proceeses called
      character*7 pnames(NAMESS)
      common /porder/ pnum, iorder(NAMESS)
      common /pnames/ pnames                                            ! The names of the processes called
c
      parameter ( NSETPROC = 1)                                         ! The number of sets of procs
      parameter ( MAXSETL = 4)                                          ! The max. length of a set
      parameter ( NFKPROC = 4)
c      character*7 setproc(NSETPROC)
      character*7 holdproc(MAXSETL)
      character*7 fkproc(NFKPROC)                                       ! List of known FK processes
c
      character*7 cpproc
      integer     holdnum
      logical     found
c
      save  setproc, fkproc
c
c      data  setproc/'FK'/
      data  fkproc / 'FK2TX ','FKMIGR','FKFILT','TX2FK'/
c
      cpproc = ' '
      cpproc = inproc
      nchars = lenstr(cpproc)
      call upcase(cpproc,nchars)
c
c.. Check to see if it is a process set
      if (cpproc.eq.'FK') then
        do 5 i = 1,NFKPROC
   5      holdproc(i) = fkproc(i)
        holdnum = NFKPROC
      else
        holdproc(1) = cpproc
        holdnum     = 1
      endif
c
c Check the procs list for the given process or set of processes
      i = 1
      found = .FALSE.
 10   if ( (i.le.pnum) .AND. .not.found) then
        j = 1
 15     if ( (j.le.holdnum) .AND. .not.found) then
          if (pnames(i).eq.holdproc(j) ) found = .TRUE.
          j = j + 1
          GO TO 15
        endif
        i = i + 1
        GO TO 10
      endif
c
      chkpra = found
      if (chkpra) then
c        print *,'CHKPRA found  ', inproc
      else
c        print *,'CHKPRA did not find ',inproc
      endif
      return
      end

max_stars_repo_name: apthorpe/sioseis

max_issues_repo_name: henkart/sioseis

max_forks_repo_name: henkart/sioseis


text:
! tl_ad_integrator.f90
!
!> Tangent Linear (TL) and Adjoint (AD) model versions of MAOOAM.
!> Integrators module.
!
!> @copyright
!> 2016 Lesley De Cruz & Jonathan Demaeyer.
!> See LICENSE.txt for license information.
!
!---------------------------------------------------------------------------!
!
!>  @remark
!>  This module actually contains the Heun algorithm routines.
!>  The user can modify it according to its preferred integration scheme.
!>  For higher-order schemes, additional buffers will probably have to be defined.
!
!---------------------------------------------------------------------------

MODULE tl_ad_integrator

  USE params, only: ndim
  USE tl_ad_tensor, only: ad,tl
  IMPLICIT NONE

  PRIVATE

  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_y1 !< Buffer to hold the intermediate position (Heun algorithm) of the tangent linear model
  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_f0 !< Buffer to hold tendencies at the initial position of the tangent linear model
  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: buf_f1 !< Buffer to hold tendencies at the intermediate position of the tangent linear model

  PUBLIC :: init_tl_ad_integrator, ad_step, tl_step

CONTAINS

  !> Routine to initialise the integration buffers.
  SUBROUTINE init_tl_ad_integrator
    INTEGER :: AllocStat
    ALLOCATE(buf_y1(0:ndim),buf_f0(0:ndim),buf_f1(0:ndim),STAT=AllocStat)
    IF (AllocStat /= 0) then
      print *, "init_tl_ad_integrator:: AllocStat = ", AllocStat
      STOP "*** Not enough memory ! ***"
    ENDIF
  END SUBROUTINE init_tl_ad_integrator

  !-----------------------------------------------------!
  !                                                     !
  ! Adjoint model integrator                            !
  !                                                     !
  !-----------------------------------------------------!

  !> Routine to perform an integration step (Heun algorithm) of the adjoint model. The incremented time is returned.
  !> @param y Initial point.
  !> @param ystar Adjoint model at the point ystar.
  !> @param t Actual integration time
  !> @param dt Integration timestep.
  !> @param res Final point after the step.
  SUBROUTINE ad_step(y,ystar,t,dt,res)
    REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y,ystar
    REAL(KIND=8), INTENT(INOUT) :: t
    REAL(KIND=8), INTENT(IN) :: dt
    REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res

    CALL ad(t,ystar,y,buf_f0)
    buf_y1 = y+dt*buf_f0
    CALL ad(t+dt,ystar,buf_y1,buf_f1)
    res=y+0.5*(buf_f0+buf_f1)*dt
    t=t+dt
  END SUBROUTINE ad_step

  !-----------------------------------------------------!
  !                                                     !
  ! Tangent linear model integrator                     !
  !                                                     !
  !-----------------------------------------------------!

  !> Routine to perform an integration step (Heun algorithm) of the tangent linear model. The incremented time is returned.
  !> @param y Initial point.
  !> @param ystar Adjoint model at the point ystar.
  !> @param t Actual integration time
  !> @param dt Integration timestep.
  !> @param res Final point after the step.
  SUBROUTINE tl_step(y,ystar,t,dt,res)
    REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: y,ystar
    REAL(KIND=8), INTENT(INOUT) :: t
    REAL(KIND=8), INTENT(IN) :: dt
    REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: res

    CALL tl(t,ystar,y,buf_f0)
    buf_y1 = y+dt*buf_f0
    CALL tl(t+dt,ystar,buf_y1,buf_f1)
    res=y+0.5*(buf_f0+buf_f1)*dt
    t=t+dt
  END SUBROUTINE tl_step

END MODULE tl_ad_integrator

max_stars_repo_name: JCSDA/MAOOAM

max_issues_repo_name: JCSDA/MAOOAM

max_forks_repo_name: JCSDA/MAOOAM


submodule (random) randinit

implicit none (type, external)

contains

module procedure rand_init

call random_init(repeatable, image_distinct)

end procedure rand_init

end submodule randinit


!-------------------------------------------------------------------------------
!> module administrator for restart
!!
!! @par Description
!!          Restart administrator
!!
!! @author Team SCALE
!<
!-------------------------------------------------------------------------------
module mod_admin_restart
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: ADMIN_restart_setup

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  logical,               public :: RESTART_RUN          = .false.   !< is this run restart?
  logical,               public :: RESTART_OUTPUT       = .false.   !< output restart file?

  character(len=H_LONG), public :: RESTART_IN_BASENAME  = ''        !< basename of the restart file
  character(len=H_LONG), public :: RESTART_OUT_BASENAME = ''        !< basename of the output file
  character(len=H_MID),  public :: RESTART_OUT_TITLE    = ''        !< title    of the output file
  character(len=H_MID),  public :: RESTART_OUT_DTYPE    = 'DEFAULT' !< REAL4 or REAL8

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  !> Setup
  subroutine ADMIN_restart_setup
    use scale_process, only: &
       PRC_MPIstop
    use mod_atmos_vars, only: &
       ATMOS_RESTART_OUTPUT,       &
       ATMOS_RESTART_IN_BASENAME,  &
       ATMOS_RESTART_OUT_BASENAME, &
       ATMOS_RESTART_OUT_TITLE,    &
       ATMOS_RESTART_OUT_DTYPE
    use mod_atmos_dyn_vars, only: &
       ATMOS_DYN_RESTART_OUTPUT,       &
       ATMOS_DYN_RESTART_IN_BASENAME,  &
       ATMOS_DYN_RESTART_OUT_BASENAME, &
       ATMOS_DYN_RESTART_OUT_TITLE,    &
       ATMOS_DYN_RESTART_OUT_DTYPE
    use mod_atmos_phy_cp_vars, only: &
       ATMOS_PHY_CP_RESTART_OUTPUT,       &
       ATMOS_PHY_CP_RESTART_IN_BASENAME,  &
       ATMOS_PHY_CP_RESTART_OUT_BASENAME, &
       ATMOS_PHY_CP_RESTART_OUT_TITLE,    &
       ATMOS_PHY_CP_RESTART_OUT_DTYPE
    use mod_atmos_phy_mp_vars, only: &
       ATMOS_PHY_MP_RESTART_OUTPUT,       &
       ATMOS_PHY_MP_RESTART_IN_BASENAME,  &
       ATMOS_PHY_MP_RESTART_OUT_BASENAME, &
       ATMOS_PHY_MP_RESTART_OUT_TITLE,    &
       ATMOS_PHY_MP_RESTART_OUT_DTYPE
    use mod_atmos_phy_rd_vars, only: &
       ATMOS_PHY_RD_RESTART_OUTPUT,       &
       ATMOS_PHY_RD_RESTART_IN_BASENAME,  &
       ATMOS_PHY_RD_RESTART_OUT_BASENAME, &
       ATMOS_PHY_RD_RESTART_OUT_TITLE,    &
       ATMOS_PHY_RD_RESTART_OUT_DTYPE
    use mod_atmos_phy_sf_vars, only: &
       ATMOS_PHY_SF_RESTART_OUTPUT,       &
       ATMOS_PHY_SF_RESTART_IN_BASENAME,  &
       ATMOS_PHY_SF_RESTART_OUT_BASENAME, &
       ATMOS_PHY_SF_RESTART_OUT_TITLE,    &
       ATMOS_PHY_SF_RESTART_OUT_DTYPE
    use mod_atmos_phy_tb_vars, only: &
       ATMOS_PHY_TB_RESTART_OUTPUT,       &
       ATMOS_PHY_TB_RESTART_IN_BASENAME,  &
       ATMOS_PHY_TB_RESTART_OUT_BASENAME, &
       ATMOS_PHY_TB_RESTART_OUT_TITLE,    &
       ATMOS_PHY_TB_RESTART_OUT_DTYPE
    use mod_atmos_phy_ch_vars, only: &
       ATMOS_PHY_CH_RESTART_OUTPUT,       &
       ATMOS_PHY_CH_RESTART_IN_BASENAME,  &
       ATMOS_PHY_CH_RESTART_OUT_BASENAME, &
       ATMOS_PHY_CH_RESTART_OUT_TITLE,    &
       ATMOS_PHY_CH_RESTART_OUT_DTYPE
    use mod_atmos_phy_ae_vars, only: &
       ATMOS_PHY_AE_RESTART_OUTPUT,       &
       ATMOS_PHY_AE_RESTART_IN_BASENAME,  &
       ATMOS_PHY_AE_RESTART_OUT_BASENAME, &
       ATMOS_PHY_AE_RESTART_OUT_TITLE,    &
       ATMOS_PHY_AE_RESTART_OUT_DTYPE
    use mod_ocean_vars, only: &
       OCEAN_RESTART_OUTPUT,       &
       OCEAN_RESTART_IN_BASENAME,  &
       OCEAN_RESTART_OUT_BASENAME, &
       OCEAN_RESTART_OUT_TITLE,    &
       OCEAN_RESTART_OUT_DTYPE
    use mod_land_vars, only: &
       LAND_RESTART_OUTPUT,       &
       LAND_RESTART_IN_BASENAME,  &
       LAND_RESTART_OUT_BASENAME, &
       LAND_RESTART_OUT_TITLE,    &
       LAND_RESTART_OUT_DTYPE
    use mod_urban_vars, only: &
       URBAN_RESTART_OUTPUT,       &
       URBAN_RESTART_IN_BASENAME,  &
       URBAN_RESTART_OUT_BASENAME, &
       URBAN_RESTART_OUT_TITLE,    &
       URBAN_RESTART_OUT_DTYPE

    implicit none

    NAMELIST / PARAM_RESTART / &
       RESTART_RUN,          &
       RESTART_OUTPUT,       &
       RESTART_IN_BASENAME,  &
       RESTART_OUT_BASENAME, &
       RESTART_OUT_TITLE,    &
       RESTART_OUT_DTYPE

    integer :: ierr
    !---------------------------------------------------------------------------

    if( IO_L ) write(IO_FID_LOG,*)
    if( IO_L ) write(IO_FID_LOG,*) '++++++ Module[ADMIN] / Categ[RESTART] / Origin[SCALE-LES]'

    !--- read namelist
    rewind(IO_FID_CONF)
    read(IO_FID_CONF,nml=PARAM_RESTART,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( IO_L ) write(IO_FID_LOG,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_RESTART. Check!'
       call PRC_MPIstop
    endif
    if( IO_LNML ) write(IO_FID_LOG,nml=PARAM_RESTART)

    !--- set default output switch
    ATMOS_RESTART_OUTPUT        = RESTART_OUTPUT
    ATMOS_DYN_RESTART_OUTPUT    = RESTART_OUTPUT
    ATMOS_PHY_CP_RESTART_OUTPUT = RESTART_OUTPUT
    ATMOS_PHY_MP_RESTART_OUTPUT = RESTART_OUTPUT
    ATMOS_PHY_RD_RESTART_OUTPUT = RESTART_OUTPUT
    ATMOS_PHY_SF_RESTART_OUTPUT = RESTART_OUTPUT
    ATMOS_PHY_TB_RESTART_OUTPUT = RESTART_OUTPUT
    ATMOS_PHY_CH_RESTART_OUTPUT = RESTART_OUTPUT
    ATMOS_PHY_AE_RESTART_OUTPUT = RESTART_OUTPUT
    OCEAN_RESTART_OUTPUT        = RESTART_OUTPUT
    LAND_RESTART_OUTPUT         = RESTART_OUTPUT
    URBAN_RESTART_OUTPUT        = RESTART_OUTPUT

    !--- set default input filename
    if ( RESTART_IN_BASENAME /= '' ) then
       ATMOS_RESTART_IN_BASENAME         = RESTART_IN_BASENAME
       ATMOS_DYN_RESTART_IN_BASENAME     = RESTART_IN_BASENAME
       ATMOS_PHY_AE_RESTART_IN_BASENAME  = RESTART_IN_BASENAME
       ATMOS_PHY_CH_RESTART_IN_BASENAME  = RESTART_IN_BASENAME
       ATMOS_PHY_CP_RESTART_IN_BASENAME  = RESTART_IN_BASENAME
       ATMOS_PHY_MP_RESTART_IN_BASENAME  = RESTART_IN_BASENAME
       ATMOS_PHY_RD_RESTART_IN_BASENAME  = RESTART_IN_BASENAME
       ATMOS_PHY_SF_RESTART_IN_BASENAME  = RESTART_IN_BASENAME
       ATMOS_PHY_TB_RESTART_IN_BASENAME  = RESTART_IN_BASENAME
       OCEAN_RESTART_IN_BASENAME         = RESTART_IN_BASENAME
       LAND_RESTART_IN_BASENAME          = RESTART_IN_BASENAME
       URBAN_RESTART_IN_BASENAME         = RESTART_IN_BASENAME
    endif

    !--- set default output filename
    if ( RESTART_OUT_BASENAME /= '' ) then
       ATMOS_RESTART_OUT_BASENAME        = RESTART_OUT_BASENAME
       ATMOS_DYN_RESTART_OUT_BASENAME    = RESTART_OUT_BASENAME
       ATMOS_PHY_CP_RESTART_OUT_BASENAME = RESTART_OUT_BASENAME
       ATMOS_PHY_MP_RESTART_OUT_BASENAME = RESTART_OUT_BASENAME
       ATMOS_PHY_RD_RESTART_OUT_BASENAME = RESTART_OUT_BASENAME
       ATMOS_PHY_SF_RESTART_OUT_BASENAME = RESTART_OUT_BASENAME
       ATMOS_PHY_TB_RESTART_OUT_BASENAME = RESTART_OUT_BASENAME
       ATMOS_PHY_CH_RESTART_OUT_BASENAME = RESTART_OUT_BASENAME
       ATMOS_PHY_AE_RESTART_OUT_BASENAME = RESTART_OUT_BASENAME
       OCEAN_RESTART_OUT_BASENAME        = RESTART_OUT_BASENAME
       LAND_RESTART_OUT_BASENAME         = RESTART_OUT_BASENAME
       URBAN_RESTART_OUT_BASENAME        = RESTART_OUT_BASENAME
    endif

    !--- set default output title
    if ( RESTART_OUT_TITLE /= '' ) then
       ATMOS_RESTART_OUT_TITLE           = RESTART_OUT_TITLE
       ATMOS_DYN_RESTART_OUT_TITLE       = RESTART_OUT_TITLE
       ATMOS_PHY_CP_RESTART_OUT_TITLE    = RESTART_OUT_TITLE
       ATMOS_PHY_MP_RESTART_OUT_TITLE    = RESTART_OUT_TITLE
       ATMOS_PHY_RD_RESTART_OUT_TITLE    = RESTART_OUT_TITLE
       ATMOS_PHY_SF_RESTART_OUT_TITLE    = RESTART_OUT_TITLE
       ATMOS_PHY_TB_RESTART_OUT_TITLE    = RESTART_OUT_TITLE
       ATMOS_PHY_CH_RESTART_OUT_TITLE    = RESTART_OUT_TITLE
       ATMOS_PHY_AE_RESTART_OUT_TITLE    = RESTART_OUT_TITLE
       OCEAN_RESTART_OUT_TITLE           = RESTART_OUT_TITLE
       LAND_RESTART_OUT_TITLE            = RESTART_OUT_TITLE
       URBAN_RESTART_OUT_TITLE           = RESTART_OUT_TITLE
    endif

    !--- set default output data type
    if ( RESTART_OUT_DTYPE /= '' ) then
       ATMOS_RESTART_OUT_DTYPE           = RESTART_OUT_DTYPE
       ATMOS_DYN_RESTART_OUT_DTYPE       = RESTART_OUT_DTYPE
       ATMOS_PHY_AE_RESTART_OUT_DTYPE    = RESTART_OUT_DTYPE
       ATMOS_PHY_CH_RESTART_OUT_DTYPE    = RESTART_OUT_DTYPE
       ATMOS_PHY_CP_RESTART_OUT_DTYPE    = RESTART_OUT_DTYPE
       ATMOS_PHY_MP_RESTART_OUT_DTYPE    = RESTART_OUT_DTYPE
       ATMOS_PHY_RD_RESTART_OUT_DTYPE    = RESTART_OUT_DTYPE
       ATMOS_PHY_SF_RESTART_OUT_DTYPE    = RESTART_OUT_DTYPE
       ATMOS_PHY_TB_RESTART_OUT_DTYPE    = RESTART_OUT_DTYPE
       OCEAN_RESTART_OUT_DTYPE           = RESTART_OUT_DTYPE
       LAND_RESTART_OUT_DTYPE            = RESTART_OUT_DTYPE
       URBAN_RESTART_OUT_DTYPE           = RESTART_OUT_DTYPE
    endif

    return
  end subroutine ADMIN_restart_setup

end module mod_admin_restart


!*==dlasd3.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b DLASD3 finds all square roots of the roots of the secular equation, as defined by the values in D and Z, and then updates the singular vectors by matrix multiplication. Used by sbdsdc.
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download DLASD3 + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE DLASD3( NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2,
!                          LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z,
!                          INFO )
!
!       .. Scalar Arguments ..
!       INTEGER            INFO, K, LDQ, LDU, LDU2, LDVT, LDVT2, NL, NR,
!      $                   SQRE
!       ..
!       .. Array Arguments ..
!       INTEGER            CTOT( * ), IDXC( * )
!       DOUBLE PRECISION   D( * ), DSIGMA( * ), Q( LDQ, * ), U( LDU, * ),
!      $                   U2( LDU2, * ), VT( LDVT, * ), VT2( LDVT2, * ),
!      $                   Z( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> DLASD3 finds all the square roots of the roots of the secular
!> equation, as defined by the values in D and Z.  It makes the
!> appropriate calls to DLASD4 and then updates the singular
!> vectors by matrix multiplication.
!>
!> This code makes very mild assumptions about floating point
!> arithmetic. It will work on machines with a guard digit in
!> add/subtract, or on those binary machines without guard digits
!> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
!> It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
!> DLASD3 is called from DLASD1.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] NL
!> \verbatim
!>          NL is INTEGER
!>         The row dimension of the upper block.  NL >= 1.
!> \endverbatim
!>
!> \param[in] NR
!> \verbatim
!>          NR is INTEGER
!>         The row dimension of the lower block.  NR >= 1.
!> \endverbatim
!>
!> \param[in] SQRE
!> \verbatim
!>          SQRE is INTEGER
!>         = 0: the lower block is an NR-by-NR square matrix.
!>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
!>
!>         The bidiagonal matrix has N = NL + NR + 1 rows and
!>         M = N + SQRE >= N columns.
!> \endverbatim
!>
!> \param[in] K
!> \verbatim
!>          K is INTEGER
!>         The size of the secular equation, 1 =< K = < N.
!> \endverbatim
!>
!> \param[out] D
!> \verbatim
!>          D is DOUBLE PRECISION array, dimension(K)
!>         On exit the square roots of the roots of the secular equation,
!>         in ascending order.
!> \endverbatim
!>
!> \param[out] Q
!> \verbatim
!>          Q is DOUBLE PRECISION array, dimension (LDQ,K)
!> \endverbatim
!>
!> \param[in] LDQ
!> \verbatim
!>          LDQ is INTEGER
!>         The leading dimension of the array Q.  LDQ >= K.
!> \endverbatim
!>
!> \param[in,out] DSIGMA
!> \verbatim
!>          DSIGMA is DOUBLE PRECISION array, dimension(K)
!>         The first K elements of this array contain the old roots
!>         of the deflated updating problem.  These are the poles
!>         of the secular equation.
!> \endverbatim
!>
!> \param[out] U
!> \verbatim
!>          U is DOUBLE PRECISION array, dimension (LDU, N)
!>         The last N - K columns of this matrix contain the deflated
!>         left singular vectors.
!> \endverbatim
!>
!> \param[in] LDU
!> \verbatim
!>          LDU is INTEGER
!>         The leading dimension of the array U.  LDU >= N.
!> \endverbatim
!>
!> \param[in] U2
!> \verbatim
!>          U2 is DOUBLE PRECISION array, dimension (LDU2, N)
!>         The first K columns of this matrix contain the non-deflated
!>         left singular vectors for the split problem.
!> \endverbatim
!>
!> \param[in] LDU2
!> \verbatim
!>          LDU2 is INTEGER
!>         The leading dimension of the array U2.  LDU2 >= N.
!> \endverbatim
!>
!> \param[out] VT
!> \verbatim
!>          VT is DOUBLE PRECISION array, dimension (LDVT, M)
!>         The last M - K columns of VT**T contain the deflated
!>         right singular vectors.
!> \endverbatim
!>
!> \param[in] LDVT
!> \verbatim
!>          LDVT is INTEGER
!>         The leading dimension of the array VT.  LDVT >= N.
!> \endverbatim
!>
!> \param[in,out] VT2
!> \verbatim
!>          VT2 is DOUBLE PRECISION array, dimension (LDVT2, N)
!>         The first K columns of VT2**T contain the non-deflated
!>         right singular vectors for the split problem.
!> \endverbatim
!>
!> \param[in] LDVT2
!> \verbatim
!>          LDVT2 is INTEGER
!>         The leading dimension of the array VT2.  LDVT2 >= N.
!> \endverbatim
!>
!> \param[in] IDXC
!> \verbatim
!>          IDXC is INTEGER array, dimension ( N )
!>         The permutation used to arrange the columns of U (and rows of
!>         VT) into three groups:  the first group contains non-zero
!>         entries only at and above (or before) NL +1; the second
!>         contains non-zero entries only at and below (or after) NL+2;
!>         and the third is dense. The first column of U and the row of
!>         VT are treated separately, however.
!>
!>         The rows of the singular vectors found by DLASD4
!>         must be likewise permuted before the matrix multiplies can
!>         take place.
!> \endverbatim
!>
!> \param[in] CTOT
!> \verbatim
!>          CTOT is INTEGER array, dimension ( 4 )
!>         A count of the total number of the various types of columns
!>         in U (or rows in VT), as described in IDXC. The fourth column
!>         type is any column which has been deflated.
!> \endverbatim
!>
!> \param[in,out] Z
!> \verbatim
!>          Z is DOUBLE PRECISION array, dimension (K)
!>         The first K elements of this array contain the components
!>         of the deflation-adjusted updating row vector.
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>         = 0:  successful exit.
!>         < 0:  if INFO = -i, the i-th argument had an illegal value.
!>         > 0:  if INFO = 1, a singular value did not converge
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date June 2017
!
!> \ingroup OTHERauxiliary
!
!> \par Contributors:
!  ==================
!>
!>     Ming Gu and Huan Ren, Computer Science Division, University of
!>     California at Berkeley, USA
!>
!  =====================================================================
      SUBROUTINE DLASD3(Nl,Nr,Sqre,K,D,Q,Ldq,Dsigma,U,Ldu,U2,Ldu2,Vt,   &
     &                  Ldvt,Vt2,Ldvt2,Idxc,Ctot,Z,Info)
      IMPLICIT NONE
!*--DLASD3227
!
!  -- LAPACK auxiliary routine (version 3.7.1) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     June 2017
!
!     .. Scalar Arguments ..
      INTEGER Info , K , Ldq , Ldu , Ldu2 , Ldvt , Ldvt2 , Nl , Nr ,    &
     &        Sqre
!     ..
!     .. Array Arguments ..
      INTEGER Ctot(*) , Idxc(*)
      DOUBLE PRECISION D(*) , Dsigma(*) , Q(Ldq,*) , U(Ldu,*) ,         &
     &                 U2(Ldu2,*) , Vt(Ldvt,*) , Vt2(Ldvt2,*) , Z(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      DOUBLE PRECISION ONE , ZERO , NEGONE
      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0,NEGONE=-1.0D+0)
!     ..
!     .. Local Scalars ..
      INTEGER ctemp , i , j , jc , ktemp , m , n , nlp1 , nlp2 , nrp1
      DOUBLE PRECISION rho , temp
!     ..
!     .. External Functions ..
      DOUBLE PRECISION DLAMC3 , DNRM2
      EXTERNAL DLAMC3 , DNRM2
!     ..
!     .. External Subroutines ..
      EXTERNAL DCOPY , DGEMM , DLACPY , DLASCL , DLASD4 , XERBLA
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , SIGN , SQRT
!     ..
!     .. Executable Statements ..
!
!     Test the input parameters.
!
      Info = 0
!
      IF ( Nl<1 ) THEN
         Info = -1
      ELSEIF ( Nr<1 ) THEN
         Info = -2
      ELSEIF ( (Sqre/=1) .AND. (Sqre/=0) ) THEN
         Info = -3
      ENDIF
!
      n = Nl + Nr + 1
      m = n + Sqre
      nlp1 = Nl + 1
      nlp2 = Nl + 2
!
      IF ( (K<1) .OR. (K>n) ) THEN
         Info = -4
      ELSEIF ( Ldq<K ) THEN
         Info = -7
      ELSEIF ( Ldu<n ) THEN
         Info = -10
      ELSEIF ( Ldu2<n ) THEN
         Info = -12
      ELSEIF ( Ldvt<m ) THEN
         Info = -14
      ELSEIF ( Ldvt2<m ) THEN
         Info = -16
      ENDIF
      IF ( Info/=0 ) THEN
         CALL XERBLA('DLASD3',-Info)
         RETURN
      ENDIF
!
!     Quick return if possible
!
      IF ( K==1 ) THEN
         D(1) = ABS(Z(1))
         CALL DCOPY(m,Vt2(1,1),Ldvt2,Vt(1,1),Ldvt)
         IF ( Z(1)>ZERO ) THEN
            CALL DCOPY(n,U2(1,1),1,U(1,1),1)
         ELSE
            DO i = 1 , n
               U(i,1) = -U2(i,1)
            ENDDO
         ENDIF
         RETURN
      ENDIF
!
!     Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
!     be computed with high relative accuracy (barring over/underflow).
!     This is a problem on machines without a guard digit in
!     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
!     The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
!     which on any of these machines zeros out the bottommost
!     bit of DSIGMA(I) if it is 1; this makes the subsequent
!     subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
!     occurs. On binary machines with a guard digit (almost all
!     machines) it does not change DSIGMA(I) at all. On hexadecimal
!     and decimal machines with a guard digit, it slightly
!     changes the bottommost bits of DSIGMA(I). It does not account
!     for hexadecimal or decimal machines without guard digits
!     (we know of none). We use a subroutine call to compute
!     2*DSIGMA(I) to prevent optimizing compilers from eliminating
!     this code.
!
      DO i = 1 , K
         Dsigma(i) = DLAMC3(Dsigma(i),Dsigma(i)) - Dsigma(i)
      ENDDO
!
!     Keep a copy of Z.
!
      CALL DCOPY(K,Z,1,Q,1)
!
!     Normalize Z.
!
      rho = DNRM2(K,Z,1)
      CALL DLASCL('G',0,0,rho,ONE,K,1,Z,K,Info)
      rho = rho*rho
!
!     Find the new singular values.
!
      DO j = 1 , K
         CALL DLASD4(K,j,Dsigma,Z,U(1,j),rho,D(j),Vt(1,j),Info)
!
!        If the zero finder fails, report the convergence failure.
!
         IF ( Info/=0 ) RETURN
      ENDDO
!
!     Compute updated Z.
!
      DO i = 1 , K
         Z(i) = U(i,K)*Vt(i,K)
         DO j = 1 , i - 1
            Z(i) = Z(i)                                                 &
     &             *(U(i,j)*Vt(i,j)/(Dsigma(i)-Dsigma(j))/(Dsigma(i)+   &
     &             Dsigma(j)))
         ENDDO
         DO j = i , K - 1
            Z(i) = Z(i)                                                 &
     &             *(U(i,j)*Vt(i,j)/(Dsigma(i)-Dsigma(j+1))/(Dsigma(i)  &
     &             +Dsigma(j+1)))
         ENDDO
         Z(i) = SIGN(SQRT(ABS(Z(i))),Q(i,1))
      ENDDO
!
!     Compute left singular vectors of the modified diagonal matrix,
!     and store related information for the right singular vectors.
!
      DO i = 1 , K
         Vt(1,i) = Z(1)/U(1,i)/Vt(1,i)
         U(1,i) = NEGONE
         DO j = 2 , K
            Vt(j,i) = Z(j)/U(j,i)/Vt(j,i)
            U(j,i) = Dsigma(j)*Vt(j,i)
         ENDDO
         temp = DNRM2(K,U(1,i),1)
         Q(1,i) = U(1,i)/temp
         DO j = 2 , K
            jc = Idxc(j)
            Q(j,i) = U(jc,i)/temp
         ENDDO
      ENDDO
!
!     Update the left singular vector matrix.
!
      IF ( K==2 ) THEN
         CALL DGEMM('N','N',n,K,K,ONE,U2,Ldu2,Q,Ldq,ZERO,U,Ldu)
         GOTO 100
      ENDIF
      IF ( Ctot(1)>0 ) THEN
         CALL DGEMM('N','N',Nl,K,Ctot(1),ONE,U2(1,2),Ldu2,Q(2,1),Ldq,   &
     &              ZERO,U(1,1),Ldu)
         IF ( Ctot(3)>0 ) THEN
            ktemp = 2 + Ctot(1) + Ctot(2)
            CALL DGEMM('N','N',Nl,K,Ctot(3),ONE,U2(1,ktemp),Ldu2,       &
     &                 Q(ktemp,1),Ldq,ONE,U(1,1),Ldu)
         ENDIF
      ELSEIF ( Ctot(3)>0 ) THEN
         ktemp = 2 + Ctot(1) + Ctot(2)
         CALL DGEMM('N','N',Nl,K,Ctot(3),ONE,U2(1,ktemp),Ldu2,Q(ktemp,1)&
     &              ,Ldq,ZERO,U(1,1),Ldu)
      ELSE
         CALL DLACPY('F',Nl,K,U2,Ldu2,U,Ldu)
      ENDIF
      CALL DCOPY(K,Q(1,1),Ldq,U(nlp1,1),Ldu)
      ktemp = 2 + Ctot(1)
      ctemp = Ctot(2) + Ctot(3)
      CALL DGEMM('N','N',Nr,K,ctemp,ONE,U2(nlp2,ktemp),Ldu2,Q(ktemp,1), &
     &           Ldq,ZERO,U(nlp2,1),Ldu)
!
!     Generate the right singular vectors.
!
 100  DO i = 1 , K
         temp = DNRM2(K,Vt(1,i),1)
         Q(i,1) = Vt(1,i)/temp
         DO j = 2 , K
            jc = Idxc(j)
            Q(i,j) = Vt(jc,i)/temp
         ENDDO
      ENDDO
!
!     Update the right singular vector matrix.
!
      IF ( K==2 ) THEN
         CALL DGEMM('N','N',K,m,K,ONE,Q,Ldq,Vt2,Ldvt2,ZERO,Vt,Ldvt)
         RETURN
      ENDIF
      ktemp = 1 + Ctot(1)
      CALL DGEMM('N','N',K,nlp1,ktemp,ONE,Q(1,1),Ldq,Vt2(1,1),Ldvt2,    &
     &           ZERO,Vt(1,1),Ldvt)
      ktemp = 2 + Ctot(1) + Ctot(2)
      IF ( ktemp<=Ldvt2 ) CALL DGEMM('N','N',K,nlp1,Ctot(3),ONE,        &
     &                               Q(1,ktemp),Ldq,Vt2(ktemp,1),Ldvt2, &
     &                               ONE,Vt(1,1),Ldvt)
!
      ktemp = Ctot(1) + 1
      nrp1 = Nr + Sqre
      IF ( ktemp>1 ) THEN
         DO i = 1 , K
            Q(i,ktemp) = Q(i,1)
         ENDDO
         DO i = nlp2 , m
            Vt2(ktemp,i) = Vt2(1,i)
         ENDDO
      ENDIF
      ctemp = 1 + Ctot(2) + Ctot(3)
      CALL DGEMM('N','N',K,nrp1,ctemp,ONE,Q(1,ktemp),Ldq,Vt2(ktemp,nlp2)&
     &           ,Ldvt2,ZERO,Vt(1,nlp2),Ldvt)
!
!
!     End of DLASD3
!
      END SUBROUTINE DLASD3


text:
c     # ----------------------------------------------------------------------
c     # This handles the boundary conditions at the block
c     # corners for the pillow sphere.
c     # ----------------------------------------------------------------------

      subroutine fclaw2d_pillow46_copy_block_corner
     &      (mx,my,mbc,meqn,
     &      qthis, qneighbor, icorner, iblock)
      implicit none

      integer mx, my, mbc, meqn, icorner, iblock
      double precision qthis(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)
      double precision qneighbor(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)

      integer mq, ibc, jbc

      do mq = 1,meqn
         if (icorner .eq. 0) then
            do ibc = 1,mbc
               do jbc = 1,mbc
                  qthis(1-ibc,1-jbc,mq) = qneighbor(ibc,jbc,mq)
                  qneighbor(1-ibc,1-jbc,mq) = qthis(ibc,jbc,mq)
               enddo
            enddo
         elseif (icorner .eq. 1) then
            do ibc = 1,mbc
               do jbc = 1,mbc
                  qthis(mx+ibc,1-jbc,mq) = qneighbor(mx+1-ibc,jbc,mq)
                  qneighbor(mx+ibc,1-jbc,mq) = qthis(mx+1-ibc,jbc,mq)
               enddo
            enddo
         elseif (icorner .eq. 2) then
            do ibc = 1,mbc
               do jbc = 1,mbc
                  qthis(1-ibc,my+jbc,mq) = qneighbor(ibc,my+1-jbc,mq)
                  qneighbor(1-ibc,my+jbc,mq) = qthis(ibc,my+1-jbc,mq)
               enddo
            enddo
         elseif (icorner .eq. 3) then
            do ibc = 1,mbc
               do jbc = 1,mbc
                  qthis(mx+ibc,my+jbc,mq) =
     &                  qneighbor(mx+1-ibc,my+1-jbc,mq)
                  qneighbor(mx+ibc,my+jbc,mq) =
     &                  qthis(mx+1-ibc,my+1-jbc,mq)
               enddo
            enddo
         endif
      enddo
      end

      subroutine fclaw2d_pillow46_average_block_corner
     &      (mx,my,mbc,meqn,
     &      refratio, qcoarse, qfine, areacoarse, areafine,
     &      icorner,iblock)
      implicit none

      integer mx, my, mbc, meqn, icorner, iblock, refratio
      double precision areacoarse(-mbc:mx+mbc+1,-mbc:my+mbc+1)
      double precision areafine(-mbc:mx+mbc+1,-mbc:my+mbc+1)

      double precision qcoarse(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)
      double precision qfine(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)

      integer mq, ibc, jbc, ii, jj, ifine, jfine
      double precision sum,qf,kf,kc

      do mq = 1,meqn
         if (icorner .eq. 0) then
            do ibc = 1,mbc
               do jbc = 1,mbc
                  sum = 0
                  do ii = 1,refratio
                     do jj = 1,refratio
                        ifine = (ibc-1)*refratio + ii
                        jfine = (jbc-1)*refratio + jj
                        qf = qfine(ifine,jfine,mq)
                        kf = areafine(ifine,jfine)
                        sum = sum + kf*qf
                     enddo
                  enddo
                  kc = areacoarse(1-ibc,1-ibc)
                  qcoarse(1-ibc,1-jbc,mq) = sum/kc
               enddo
            enddo
         elseif (icorner .eq. 1) then
            do ibc = 1,mbc
               do jbc = 1,mbc
                  sum = 0
                  do ii = 1,refratio
                     do jj = 1,refratio
                        ifine = (ibc-1)*refratio + ii
                        jfine = (jbc-1)*refratio + jj
                        qf = qfine(mx+1-ifine,jfine,mq)
                        kf = areafine(mx+1-ifine,jfine)
                        sum = sum + kf*qf
                     enddo
                  enddo
                  kc = areacoarse(mx+ibc,1-jbc)
                  qcoarse(mx+ibc,1-jbc,mq) = sum/kc
               enddo
            enddo
         elseif (icorner .eq. 2) then
            do ibc = 1,mbc
               do jbc = 1,mbc
                  sum = 0
                  do ii = 1,refratio
                     do jj = 1,refratio
                        ifine = (ibc-1)*refratio + ii
                        jfine = (jbc-1)*refratio + jj
                        qf = qfine(ifine,my+1-jfine,mq)
                        kf = areafine(ifine,my+1-jfine)
                        sum = sum + kf*qf
                     enddo
                  enddo
                  kc = areacoarse(1-ibc,my+jbc)
                  qcoarse(1-ibc,my+jbc,mq) = sum/kc
               enddo
            enddo
         elseif (icorner .eq. 3) then
            do ibc = 1,mbc
               do jbc = 1,mbc
                  sum = 0
                  do ii = 1,refratio
                     do jj = 1,refratio
                        ifine = (ibc-1)*refratio + ii
                        jfine = (jbc-1)*refratio + jj
                        qf = qfine(mx+1-ifine,my+1-jfine,mq)
                        kf = areafine(mx+1-ifine,my+1-jfine)
                        sum = sum + kf*qf
                     enddo
                  enddo
                  kc = areacoarse(mx+ibc,my+jbc)
                  qcoarse(mx+ibc,my+jbc,mq) = sum/kc
               enddo
            enddo
         endif
      enddo

      end

      subroutine fclaw2d_pillow46_interpolate_block_corner
     &      (mx,my,mbc,meqn,
     &      refratio, qcoarse, qfine, icorner_coarse, iblock)
      implicit none
      integer mx, my, mbc, meqn, icorner_coarse, iblock, refratio
      double precision qcoarse(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)
      double precision qfine(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)

      integer mq, ibc, jbc
      integer ic, jc, mth
      double precision gradx, grady, shiftx, shifty
      double precision sl, sr, qc, value
      double precision fclaw2d_clawpatch_compute_slopes

      mth = 5

      if (icorner_coarse .eq. 0) then
         ic = 1
         jc = 1
      elseif (icorner_coarse .eq. 1) then
         ic = mx
         jc = 1
      elseif (icorner_coarse .eq. 2) then
         ic = 1
         jc = my
      elseif (icorner_coarse .eq. 3) then
         ic = mx
         jc = my
      else
         write(6,*) 'pillow : icorner_coarse has unexpected value'
         write(6,*) 'icorner_coarse : ', icorner_coarse
         stop
      endif

c     # This may not even matter
      do mq = 1,meqn
         qc = qcoarse(ic,jc,mq)
         sl = (qc - qcoarse(ic-1,jc,mq))
         sr = (qcoarse(ic+1,jc,mq) - qc)
         gradx = fclaw2d_clawpatch_compute_slopes(sl,sr,mth)

         sl = (qc - qcoarse(ic,jc-1,mq))
         sr = (qcoarse(ic,jc+1,mq) - qc)
         grady = fclaw2d_clawpatch_compute_slopes(sl,sr,mth)

c        # Loop over fine grid ghost cells
         do ibc = 1,mbc
            do jbc = 1,mbc
c              # Fill in interpolated values on fine grid cell
               shiftx = (ibc - refratio/2.d0 - 0.5d0)/refratio
               shifty = (jbc - refratio/2.d0 - 0.5d0)/refratio

               value = qc + shiftx*gradx + shifty*grady
               if (icorner_coarse .eq. 0) then
                  qfine(1-ibc,1-jbc,mq) = value
               elseif (icorner_coarse .eq. 1) then
                  qfine(mx+ibc,1-jbc,mq) = value
               elseif (icorner_coarse .eq. 2) then
                  qfine(1-ibc,my+jbc,mq) = value
               elseif (icorner_coarse .eq. 3) then
                  qfine(mx+ibc,my+jbc,mq) = value
               endif
            enddo
         enddo
      enddo
      end

max_stars_repo_name: mattzett/forestclaw

max_issues_repo_name: mattzett/forestclaw

max_forks_repo_name: mattzett/forestclaw


text:
C
C $Id: mdqtra.f,v 1.6 2008-09-18 00:42:18 kennison Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE MDQTRA (RLAT,RLON,UVAL,VVAL)
C
        DOUBLE PRECISION RLAT,RLON,UVAL,VVAL
C
C Declare a special common block, containing only those variables that
C MDQINI needs to set to make MDQTRA, MDQTRI, and MDQTRN carry out the
C transformation in effect at the time MDQINI was called.
C
        COMMON /MAQCMN/  ALFA,COSO,COSR,CSLS,CSLT,DCSA,DCSB,DSNA,DSNB,
     +                   DTOR,DTRH,OOPI,PLNO,  PI,PIOT,ROTA,RTDD,RTOD,
     +                   SALT,SINO,SINR,SRSS,SSMO,TOPI,UCNM,UMNM,UMXM,
     +                   UOFF,URNM,VCNM,VMNM,VMXM,VOFF,VRNM,UTPA,IPRF,
     +                   IPRJ,IROD,ELPM
        DOUBLE PRECISION ALFA,COSO,COSR,CSLS,CSLT,DCSA,DCSB,DSNA,DSNB,
     +                   DTOR,DTRH,OOPI,PLNO,  PI,PIOT,ROTA,RTDD,RTOD,
     +                   SALT,SINO,SINR,SRSS,SSMO,TOPI,UCNM,UMNM,UMXM,
     +                   UOFF,URNM,VCNM,VMNM,VMXM,VOFF,VRNM,UTPA(15)
C
        INTEGER IPRF,IPRJ,IROD
C
        LOGICAL ELPM
C
        SAVE   /MAQCMN/
C
C Check for an uncleared prior error.
C
        IF (ICFELL('MDQTRA - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
C
C The call to MDQTRA is simply passed on to MDQTRN, but the values
C returned are checked to see if the point lies outside the perimeter;
C if so, the value 1.D12 is substituted for UVAL.
C
        CALL MDQTRN (RLAT,RLON,UVAL,VVAL)
        IF (ICFELL('MDQTRA',2).NE.0) RETURN
C
        IF (ELPM) THEN
          IF (((UVAL-UCNM)/URNM)**2+
     +        ((VVAL-VCNM)/VRNM)**2.GT.1.000002D0) THEN
            UVAL=1.D12
            VVAL=1.D12
          END IF
        ELSE
          IF (UVAL.LT.UMNM.OR.UVAL.GT.UMXM.OR.
     +        VVAL.LT.VMNM.OR.VVAL.GT.VMXM) THEN
            UVAL=1.D12
            VVAL=1.D12
          END IF
        END IF
C
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      subroutine cmo_addatt_def(imsgin,xmsgin,cmsgin,msgtype,nwds,
     *                          ierror_return)
C
C
C#######################################################################
C
C      PURPOSE -
C
C         This Routine Adds Attributes to the Default Table..
C
C      INPUT ARGUMENTS -
C
C         imsgin()  - Integer array of command input tokens
C         xmsgin()  - Real array of command input tokens
C         cmsgin()  - Character array of command input tokens
C         msgtype() - Integer array of command input token types
C         nwds      - Number of command input tokens
C
C      OUTPUT ARGUMENTS -
C
C         ierror_return - Error Return Code (==0 ==> OK, <>0 ==> Error)
C
C      CHANGE HISTORY -
C
C         $Log: cmo_addatt_def.f,v $
C         Revision 2.00  2007/11/05 19:45:47  spchu
C         Import to CVS
C
CPVCS
CPVCS       Rev 1.3   07 Jul 2006 08:48:36   gable
CPVCS    Modified screen output to reduce blank space.
CPVCS
CPVCS       Rev 1.2   10 Apr 2001 11:04:10   dcg
CPVCS    shorten too long name
CPVCS
CPVCS       Rev 1.1   21 Jan 2000 17:02:38   dcg
CPVCS
CPVCS
CPVCS       Rev 1.8   Tue Nov 02 19:07:48 1999   jtg
CPVCS    fixed hardwired character*32 for cmsgin
CPVCS
CPVCS       Rev 1.7   Mon Apr 14 16:39:52 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.6   Mon Nov 18 11:00:50 1996   dcg
CPVCS    remove character literals from argument lists
CPVCS
CPVCS       Rev 1.5   09/11/95 14:43:14   het
CPVCS    Change to the storage block based CMO stuff.
CPVCS
CPVCS       Rev 1.4   03/15/95 15:22:16   ejl
CPVCS    Finished installing the defaults.
CPVCS
CPVCS       Rev 1.3   02/16/95 10:22:58   ejl
CPVCS    Put return afer end statement.
CPVCS
CPVCS       Rev 1.2   02/16/95 09:55:40   ejl
CPVCS    Fixed bugs, fixed hole in the Create command.
CPVCS    Added commands MODATT, LENGTH, MEMORY, & COMPRESS.
CPVCS
CPVCS       Rev 1.1   02/10/95 14:06:28   ejl
CPVCS    Fix bugs left from last update.
C
C#######################################################################
C
      implicit none
C
C#######################################################################
C
      include 'cmo_lg.h'
C
C#######################################################################
C
      integer nwds, imsgin(nwds), msgtype(nwds)
      REAL*8 xmsgin(nwds)
      character*(*) cmsgin(nwds)
C
      integer ierror_return
C
C#######################################################################
C
C     LOCAL VARIABLE DEFINITION
C
      character*32 att_name, partname,cmo_name
      integer ierr, iatt_index, nodnx, icscode,i,natts,len,
     *  posname,postype,posrank,poslen,posint,posio,pospers,
     *  lenattlist
C
      character*132 logmess
C
C#######################################################################
C
      integer icharlnf
C
C#######################################################################
C
      partname='default_cmo_lg'
      att_name=cmsgin(4)
      call mmfindbk('cmo_attlist',partname,ipcmo_attlist,
     *                   lenattlist,icscode)
      do i=1,number_of_default_attributes
         if(cmo_attlist(number_of_params_per_att*(i-1)+1)
     *             .eq.att_name) then
C
C....    This is a existing attribute.
C
            ierror_return=-1
C
            write(logmess,'(a,a)')
     *         'CMO_ADDATT error: attribute already exist: ',
     *          att_name(1:icharlnf(att_name))
            call writloga('default',0,logmess,0,ierr)
            go to 9999
         endif
C
      enddo
C
C....    Change the Default Attribute Table.
C
      ierror_return=0
      natts=number_of_default_attributes
      number_of_default_attributes=natts+1
C
C....    See if there is enough space for new attribute
C
      if(lenattlist.lt.(natts+1)*number_of_params_per_att*4) then
        call mmincblk('cmo_attlist',partname,ipcmo_attlist,2000,icscode)
      endif
C
c.... Find postions of name, type, rank and length
      partname='define_cmo_lg'
      call mmfindbk( 'defcmo_attparam_names',partname,
     *      ipdefcmo_attparam_names,len,icscode)
      do i=1,number_of_default_attparam_name
            if(defcmo_attparam_names(i).eq.'name') posname=i
            if(defcmo_attparam_names(i).eq.'type') postype=i
            if(defcmo_attparam_names(i).eq.'rank') posrank=i
            if(defcmo_attparam_names(i).eq.'length') poslen=i
            if(defcmo_attparam_names(i).eq.'interpolation') posint=i
            if(defcmo_attparam_names(i).eq.'persistence') pospers=i
            if(defcmo_attparam_names(i).eq.'ioflag') posio=i
      enddo
      cmo_attlist(number_of_params_per_att*(natts)+posname)=att_name
      cmo_attlist(number_of_params_per_att*(natts)+postype)=cmsgin(5)
      cmo_attlist(number_of_params_per_att*(natts)+posrank)=cmsgin(6)
      cmo_attlist(number_of_params_per_att*(natts)+poslen)=cmsgin(7)
      cmo_attlist(number_of_params_per_att*(natts)+posint)=cmsgin(8)
      cmo_attlist(number_of_params_per_att*(natts)+pospers)=
     *    'temporary'
      cmo_attlist(number_of_params_per_att*(natts)+posio)=cmsgin(10)
C

9999  return
      end

max_forks_repo_name: daniellivingston/LaGriT


text:
      SUBROUTINE ENDBUF
C
C          ------------------------------------------------
C          ROUTINE NO. (  13)   VERSION (A7.3)    11:FEB:85
C          ------------------------------------------------
C
C
C          THIS RESTORES GRAPHICAL OUTPUT TO THE MAIN STREAM.
C
C
      REAL    RDATA(1)
      INTEGER IDATA(1)
C
      COMMON /T0BUFN/ KBUFR0
      COMMON /T0TRAC/ IPRINT
C
      DATA RDATA /0.0/, IDATA /0/
C
C
      IF (IPRINT.EQ.1) CALL G0MESG(115,0)
C
      KBUFR0= 0
      CALL G3LINK(4,2,-1,IDATA,RDATA)
C
      RETURN
      END

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST


text: C<8, 4, 7>(n, k, l) = A<8, 16, 23>(n, c, l*3+r*2) * W<4, 16, 3>(k, c, r);

max_stars_repo_name: KnowingNothing/ExprAutodiff

max_issues_repo_name: KnowingNothing/ExprAutodiff

max_forks_repo_name: KnowingNothing/ExprAutodiff


text:
	SUBROUTINE SNDSPL ( nparms, np, nlvin, datain,
     +			    nlvout, hghts, dataout, iret )
C************************************************************************
C* SNDSPL								*
C* 									*
C* This subroutine fits a curve to a set of input points and then	*
C* evaluates the curve for values of x.  The input points must be	*
C* strictly monotonic in x:						*
C* 									*
C*     datain (1) < datain (2) < ... < datain (nlvin)			*
C*			or						*
C*     datain (1) > datain (2) > ... > datain (nlvin)			*
C* 									*
C* SNDSPL ( NPARMS, NP, NLVIN, DATAIN, NLVOUT, HGHTS, DATAOUT, IRET )	*
C* 									*
C* Input parameters:							*
C*	NPARMS		INTEGER		Number of parameters		*
C*	NP		INTEGER		Current parameter		*
C*	NLVIN		INTEGER		Number of levels in		*
C*	DATAIN (LLMXLV) REAL		Input sounding data		*
C*	NLVOUT		INTEGER		Number of levels out		*
C*	HGHTS (NLVOUT) 	REAL		Height levels			*
C* 									*
C* Output parameters:							*
C*	DATAOUT(LLMXLV)	REAL		Output sounding data		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal			*
C**									*
C* Log:									*
C* S. Jacobs/SSAI	 4/92		Copied from GCURVE		*
C* J. Whistler/SSAI	 4/93		Cleaned up code			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'ERROR.PRM'
	INCLUDE		'CURVE.CMN'
C*
	REAL		datain(*), dataout(*), hghts(*)
C*
	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	iret = 0
C
C*	Check for enough points.
C
	IF  ( ( nlvin .lt. 3 ) .or. ( nlvout .lt. 1 ) )
     +		iret = NOPNTS
	IF  ((nlvin .gt. MAXPTS) .or. (nlvout .gt. MAXPTS))
     +		iret = NIPNTS
	IF  ( iret .ne. NORMAL ) RETURN
C
C*	Move data into internal arrays sorted in increasing order.
C*	Return with error if data is not strictly monotonic.
C
	delx = datain(1*nparms+1) - datain(0*nparms+1)
	IF  ( delx .lt. 0. ) THEN
	    ibegin = nlvin
	    iend   = 1
	    inc    = -1
	ELSE
	    ibegin = 1
	    iend   = nlvin
	    inc    = 1
	END IF
C*
	x(1) = datain( (ibegin-1)*nparms +  1 )
	y(1) = datain( (ibegin-1)*nparms + np )
	nknt  = 1
	xx    = x (1)
	DO  i = ibegin+inc, iend, inc
	    tx = datain( (i-1)*nparms +  1 )
	    ty = datain( (i-1)*nparms + np )
	    IF  ( ( .not. ERMISS (tx) ) .and.
     +            ( .not. ERMISS (ty) ) ) THEN
		IF  ( tx .le. xx ) THEN
		    iret = NIMONO
		    RETURN
		ELSE
		    nknt     = nknt + 1
		    x (nknt) = tx
		    y (nknt) = ty
		    xx       = tx
		END IF
	    END IF
	END DO
C
C*	Check to see if there are still at least two points.
C
	IF  ( nknt .lt. 2 ) THEN
	    iret = NOPNTS
	    RETURN
	END IF
C
C*	Generate coefficients for evaluating curve.
C
	jend = 1
	CALL SNDFSP ( x, y, nknt, jend, coeffs, iret )
C
	IF  ( iret .ne. NORMAL ) THEN
	    iret = NORMAL
	    DO  i = 1, nlvout
		dataout( (i-1)*nparms + np ) = RMISSD
	    END DO
	    RETURN
	END IF
C
C*	Transform evaluating points to normalized coordinates.
C
	DO  i = 1, nlvout
	    xout(i) = hghts(i)
	END DO
C
C*	Evaluate curve.
C
	CALL SNDESP ( x, nknt, coeffs, xout, nlvout, yout, iret )
C
C*	Move points into output array.
C
	DO  i = 1, nlvout
	    dataout( (i-1)*nparms + np ) = yout(i)
	END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


! { dg-do run }
!
! Test the fix for PR84115 comment #1.
!
! Contributed by G Steinmetz  <gscfq@t-online.de>
!
  character(:), allocatable :: dum
  dum = "s1"
  call s1 (dum)
  dum = "s2"
  call s2 (dum)
  dum = "s3"
  call s3 (dum)
contains
  subroutine s1(x)
    character(:), allocatable :: x
    associate (y => x//x)
      if (y .ne. x//x) stop 1
    end associate
  end

  subroutine s2(x)
    character(:), allocatable :: x
    associate (y => [x])
      if (any(y .ne. [x])) stop 2
    end associate
  end

  subroutine s3(x)
    character(:), allocatable :: x
    associate (y => [x,x])
      if (any(y .ne. [x,x])) stop 3
    end associate
  end
end


text:
      SUBROUTINE CALBAL1
C
C CHANGE RECORD
C **  SUBROUTINES CALBAL CALCULATE GLOBAL VOLUME, MASS, MOMENTUM,
C **  AND ENERGY BALANCES
C
      USE GLOBAL
	IMPLICIT NONE
	INTEGER::L,LN,K
      IF(NBAL.GT.1) RETURN
C
C **  INITIALIZE VOLUME, SALT MASS, DYE MASS, MOMENTUM, KINETIC ENERGY
C **  AND POTENTIAL ENERGY, AND ASSOCIATED FLUXES
C
      VOLBEG=0.
      SALBEG=0.
      DYEBEG=0.
      UMOBEG=0.
      VMOBEG=0.
      UUEBEG=0.
      VVEBEG=0.
      PPEBEG=0.
      BBEBEG=0.
      VOLOUT=0.
      SALOUT=0.
      DYEOUT=0.
      UMOOUT=0.
      VMOOUT=0.
      UUEOUT=0.
      VVEOUT=0.
      PPEOUT=0.
      BBEOUT=0.
      DO L=2,LA
        LN=LNC(L)
        VOLBEG=VOLBEG+SPB(L)*DXYP(L)*HP(L)
        UMOBEG=UMOBEG+SPB(L)*0.5*DXYP(L)*HP(L)*(DYIU(L)*HUI(L)*UHDYE(L)
     &      +DYIU(LEAST(L))*HUI(LEAST(L))*UHDYE(LEAST(L)))
        VMOBEG=VMOBEG+SPB(L)*0.5*DXYP(L)*HP(L)*(DXIV(L)*HVI(L)*VHDXE(L)
     &      +DXIV(LN)*HVI(LN)*VHDXE(LN))
        PPEBEG=PPEBEG+SPB(L)*0.5*DXYP(L)
     &      *(GI*P(L)*P(L)-G*BELV(L)*BELV(L))
      ENDDO
      AMOBEG=SQRT(UMOBEG*UMOBEG+VMOBEG*VMOBEG)
      DO K=1,KC
        DO L=2,LA
          LN=LNC(L)
          SALBEG=SALBEG+SCB(L)*DXYP(L)*HP(L)*SAL(L,K)*DZC(K)
          DYEBEG=DYEBEG+SCB(L)*DXYP(L)*HP(L)*DYE(L,K)*DZC(K)
          UUEBEG=UUEBEG+SPB(L)*0.125*DXYP(L)*HP(L)*DZC(K)
     &        *( (U(L,K)+U(LEAST(L),K))*(U(L,K)+U(LEAST(L),K)) )
          VVEBEG=VVEBEG+SPB(L)*0.125*DXYP(L)*HP(L)*DZC(K)
     &        *( (V(L,K)+V(LN,K))*(V(L,K)+V(LN,K)) )
          BBEBEG=BBEBEG+SPB(L)*GP*DXYP(L)*HP(L)*DZC(K)*( BELV(L)
     &        +0.5*HP(L)*(Z(K)+Z(K-1)) )*B(L,K)
        ENDDO
      ENDDO
      RETURN
      END