      SUBROUTINE MINMAX(A,IDIM,JDIM,AMIN,AMAX)
C
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    MINMAX      RETURNS (MIN,MAX) IN ARRAY
C   PRGRMMR: TREADON         ORG: W/NP2      DATE: 93-01-08
C     
C ABSTRACT:
C     THIS ROUTINE FINDS THE MINIMUM AND MAXIMUM VALUES 
C     IN THE PASSED ARRAY "A" WHICH CONTAINS NPTS ELEMENTS.
C   .     
C     
C PROGRAM HISTORY LOG:
C   ??-??-??  DAVID PLUMMER - SUBROUTINE MINMAX IN ETAPACKC.
C   93-01-08  RUSS TREADON - ADDED COMMENTS.
C     
C USAGE:    CALL MINMAX(A,IDIM,JDIM,AMIN,AMAX)
C   INPUT ARGUMENT LIST:
C     A        - INPUT ARRAY TO SCAN.
C     IDIM,JDIM- DIMENSIONS OF ARRAY "A".
C
C   OUTPUT ARGUMENT LIST: 
C     AMIN     - MINIMUM VALUE IN ARRAY "A".
C     AMAX     - MAXIMUM VALUE IN ARRAY "A".
C     
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C     UTILITIES:
C       NONE
C     LIBRARY:
C       NONE
C     
C   ATTRIBUTES:
C     LANGUAGE: FORTRAN
C     MACHINE : CRAY C-90
C$$$  
C     
C     
C     DECLARE VARIABLES.
      REAL AMIN, AMAX
      REAL A(IDIM,JDIM)
C     
C*******************************************************************
C     START MINMAX HERE.
C     
C     INITIALIZE VARIABLES.
      AMIN =  1.E30
      AMAX = -1.E30
C     
C     FIND MIN AND MAX VALUES IN ARRAY.
C
      DO J=1,JDIM
      DO I=1,IDIM
        AMIN = AMIN1(AMIN,A(I,J))
        AMAX = AMAX1(AMAX,A(I,J))
      ENDDO
      ENDDO
C     
C     END OF ROUTINE.
      RETURN      
      END     


   
      SUBROUTINE EFILL(EGRID,IM,JM)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    EFILL       FILLS MISSING VALUES ON E-GRID
C   PRGRMMR: TREADON         ORG: W/NP2 
C     
C ABSTRACT:  
C     THIS ROUTINE REPLACES ARRAY ELEMENTS ON A E-GRID MARKED
C     BY A MISSING DATA FLAG WITH THE FIELD MEAN.  THE FIELD
C     MEAN IS COMPUTED FROM ALL NON-MISSING ARRAY VALUES.
C     THE MISSING DATA FLAG, SPVAL,  IS SET IN INCLUDE FILE
C     (AND COMMON BLOCK) OPTIONS.
C   .     
C     
C PROGRAM HISTORY LOG:
C   93-02-03  RUSS TREADON
C   98-06-01  BLACK - CONVERSION FROM 1-D TO 2-D
C     
C USAGE:    CALL EFILL(EGRID,IM,JM)
C   INPUT ARGUMENT LIST:
C     EGRID    - DATA ON STAGGERED E-GRID
C     IM,JM    - DIMENSIONS OF E-GRID
C
C   OUTPUT ARGUMENT LIST: 
C     EGRID    - SAME A INPUT EXCEPT MISSING VALUES
C                ARE REPLACED BY THE FIELD MEAN
C     
C   OUTPUT FILES:
C     STDOUT  - RUN TIME STANDARD OUT.
C     
C   SUBPROGRAMS CALLED:
C     UTILITIES:
C       NONE
C     LIBRARY:
C       COMMON   - OPTIONS
C                  IOUNIT
C     
C   ATTRIBUTES:
C     LANGUAGE: FORTRAN
C     MACHINE : CRAY C-90
C$$$  
C     
C     SET PARAMETER
      PARAMETER (SMALL=1.E-12)
C
C     DECLARE VARIABLES
      REAL EGRID(IM,JM)
      DOUBLE PRECISION SUM
C     
C     
      SPVAL= 999999.
C**********************************************************************
C     START EFILL HERE
C     
C     SEE IF THERE ARE ANY MISSING VALUES IN THE ARRAY.  IF NOT,
C     RETURN TO CALLING PROGRAM.
C
      CALL MINMAX(EGRID,IM,JM,FMIN,FMAX)
      IF (ABS(FMAX-SPVAL).GT.SMALL) RETURN
C     
C     COMPUTE FIELD MEAN, EXCLUDING MISSING POINTS.
C
      NPTS=0
      SUM =0.
C
      DO J=1,JM
      DO I=1,IM
        IF(ABS(EGRID(I,J)-SPVAL).GE.SMALL)THEN
          NPTS=NPTS+1
          SUM =SUM+EGRID(I,J)
        ENDIF
      ENDDO
      ENDDO
      IF(NPTS.EQ.0)THEN
        WRITE(6,*)'EFILL:  ALL DATA ON EGRID EQUAL SPVAL'
        RETURN
      ENDIF
      AVG=SUM/NPTS
C     
C     REPLACE MISSING VALUES WITH FIELD MEAN.
C
!$omp  parallel do
      DO J=1,JM
      DO I=1,IM
        IF(ABS(EGRID(I,J)-SPVAL).LT.SMALL)EGRID(I,J)=AVG
      ENDDO
      ENDDO
CX      CALL MINMAX(EGRID,IM,JM,FMIN,FMAX)
CX      WRITE(STDOUT,*)'EFILL:  EGRID MIN,MAX,AVG:  ',FMIN,FMAX,AVG
C     
C     END OF ROUTINE.
C     
      RETURN
      END


      SUBROUTINE COEFI(ARRX,IDIM,JDIM,I,J,C1,C2,C3,C4)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    COEFI       INTRP WGTS FOR FILLED E-GRID ROWS
C   PRGRMMR: TREADON         ORG: W/NP2 
C     
C ABSTRACT:  
C     THIS ROUTINE RETURNS WEIGHTED E-GRID VALUES FOR CUBIC
C     INTERPOLATION OF ROWS ON FILLED E-GRID.
C   .     
C     
C PROGRAM HISTORY LOG:
C   93-05-12  RUSS TREADON
C     
C USAGE:    CALL COEFI(ARRX,IDIM,JDIM,I,J,C1,C2,C3,C4)
C   INPUT ARGUMENT LIST:
C     ARRX     - ARRAY TO CONTAIN FILLED E-GRID.
C     IDIM     - 2*IM-1, FIRST DIMENSION OF ARRX
C     JDIM     - JM, SECOND DIMENSION OF ARRX
C     I        - TARGET POINT I INDEX
C     J        - TARGET POINT J INDEX
C
C   OUTPUT ARGUMENT LIST: 
C     C1,C2,C3,C4 - WEIGHTED FILLED E-GRID VALUES.
C
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C     UTILITIES:
C       NONE
C     LIBRARY:
C       NONE
C     
C   ATTRIBUTES:
C     LANGUAGE: FORTRAN
C     MACHINE : CRAY C-90
C$$$  
C
      PARAMETER (C11=.166667,C12=-.5,C13=.5,C14=-.166667)
      PARAMETER (C21=-1.,C22=3.5,C23=-4.,C24=1.5)
      PARAMETER (C31=1.833333,C32=-7.,C33=9.5,C34=-4.333333)
      PARAMETER (C41=-1.,C42=4.,C43=-6.,C44=4.)
      DIMENSION ARRX(IDIM,JDIM)
      C1=C11*ARRX(I+6,J)+C12*ARRX(I+4,J)+C13*ARRX(I+2,J)+C14*ARRX(I,J)
      C2=C21*ARRX(I+6,J)+C22*ARRX(I+4,J)+C23*ARRX(I+2,J)+C24*ARRX(I,J)
      C3=C31*ARRX(I+6,J)+C32*ARRX(I+4,J)+C33*ARRX(I+2,J)+C34*ARRX(I,J)
      C4=C41*ARRX(I+6,J)+C42*ARRX(I+4,J)+C43*ARRX(I+2,J)+C44*ARRX(I,J)
      RETURN
      END


      SUBROUTINE COEFJ(ARRX,IDIM,JDIM,I,J,C1,C2,C3,C4)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    COEFJ       INTRP WGTS FOR FILLED E-GRID COLS
C   PRGRMMR: TREADON         ORG: W/NP2      DATE: ??-??-??
C     
C ABSTRACT:  
C     THIS ROUTINE RETURNS WEIGHTED E-GRID VALUES FOR CUBIC
C     INTERPOLATION OF COLUNMS ON FILLED E-GRID.
C   .     
C     
C PROGRAM HISTORY LOG:
C   93-05-12  RUSS TREADON - ADDED DOC BLOC
C     
C USAGE:    CALL COEFJ(ARRX,IDIM,JDIM,I,J,C1,C2,C3,C4)
C   INPUT ARGUMENT LIST:
C     ARRX     - ARRAY TO CONTAIN FILLED E-GRID.
C     IDIM     - 2*IM-1, FIRST DIMENSION OF ARRX
C     JDIM     - JM, SECOND DIMENSION OF ARRX
C     I        - TARGET POINT I INDEX
C     J        - TARGET POINT J INDEX
C
C   OUTPUT ARGUMENT LIST: 
C     C1,C2,C3,C4 - WEIGHTED FILLED E-GRID VALUES.
C
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C     UTILITIES:
C       NONE
C     LIBRARY:
C       NONE
C     
C   ATTRIBUTES:
C     LANGUAGE: FORTRAN
C     MACHINE : CRAY C-90
C$$$  
C
      PARAMETER (C11=.166667,C12=-.5,C13=.5,C14=-.166667)
      PARAMETER (C21=-1.,C22=3.5,C23=-4.,C24=1.5)
      PARAMETER (C31=1.833333,C32=-7.,C33=9.5,C34=-4.333333)
      PARAMETER (C41=-1.,C42=4.,C43=-6.,C44=4.)
      DIMENSION ARRX(IDIM,JDIM)
      C1=C11*ARRX(I,J+6)+C12*ARRX(I,J+4)+C13*ARRX(I,J+2)+C14*ARRX(I,J)
      C2=C21*ARRX(I,J+6)+C22*ARRX(I,J+4)+C23*ARRX(I,J+2)+C24*ARRX(I,J)
      C3=C31*ARRX(I,J+6)+C32*ARRX(I,J+4)+C33*ARRX(I,J+2)+C34*ARRX(I,J)
      C4=C41*ARRX(I,J+6)+C42*ARRX(I,J+4)+C43*ARRX(I,J+2)+C44*ARRX(I,J)
      RETURN
      END


      SUBROUTINE FILLH(ARR1,ARR2,IDIM,JDIM)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    FILLH       FILL H POINTS ON E-GRID V FIELD
C   PRGRMMR: BLACK           ORG: W/NP2      DATE: 92-12-23       
C     
C ABSTRACT:
C     THIS ROUTINE FILLS MASS (H) POINTS ON A VELOCITY (V)
C     POINT E-GRID.
C   .     
C     
C PROGRAM HISTORY LOG:
C   ??-??-??  BLACK - ORIGINATOR
C   93-12-23  RUSS TREADON - ADDED COMMENTS AND VARIABLE
C                            GRID DIMENSIONS TO 
C                            SUBROUTINE CALL.
C   98-06-04  BLACK - CONVERSION TO 2-D
C     
C USAGE:    CALL FILLH(ARR1,ARR2,IDIM,JDIM)
C   INPUT ARGUMENT LIST:
C     ARR2     - VELOCITY POINT DATA ON STAGGERED E-GRID.
C     IDIM     - FIRST DIMENSION OF FILLED E-GRID.
C     JDIM     - SECOND DIMENSION OF FILLED E-GRID.
C
C   OUTPUT ARGUMENT LIST: 
C     ARR1     - VELOCITY POINT DATA ON FILLED E-GRID.
C     
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C     UTILITIES:
C       EFILL - FILL MISSING VALUES ON STAGGERED
C               E-GRID WITH FIELD MEAN.
C     LIBRARY:
C       COMMON
C     
C   ATTRIBUTES:
C     LANGUAGE: FORTRAN
C     MACHINE : CRAY C-90
C$$$  
C
C--------------------------------------------------------------------
                        D I M E N S I O N
     & ARR1((IDIM+1)/2,JDIM), ARR2(IDIM,JDIM)
C     
C********************************************************************
C     START FILLH HERE.
C     
C     REPLACE MISSING VALUES ON EGRID WITH FIELD MEAN
      IM=(IDIM+1)/2
      CALL EFILL(ARR1,IM,JDIM)
C
      DO J=1,JDIM
      DO I=1,IDIM
        ARR2(I,J)=0.
      ENDDO
      ENDDO
C***
C***  REASSIGN VALUES OF MODEL GRID TO V POINTS ON IDIM X JDIM GRID
C***  ON ODD ROWS.
C***
      DO J=1,JDIM,2
      DO I=2,IDIM-1,2
        II=I/2
        ARR2(I,J)=ARR1(II,J)
      ENDDO
      ENDDO
C***
C***  REASSIGN VALUES OF MODEL GRID TO V POINTS ON IDIM X JDIM GRID
C***  ON EVEN ROWS.
C***
      DO J=2,JDIM-1,2
      DO I=1,IDIM,2
        II=(I+1)/2
        ARR2(I,J)=ARR1(II,J)
      ENDDO
      ENDDO
C***
C***  INTERPOLATE VALUES TO THOSE CORRESPONDING TO H POINTS ON THE
C***  IDIM X JDIM GRID ON THE BOTTOM AND TOP ROWS.
C***
      DO I=3,IDIM-2,2
        II=I/2
        ARR2(I,1)=(ARR1(II,1)+ARR1(II+1,1))*0.5
        ARR2(I,JDIM)=(ARR1(II,JDIM)+ARR1(II+1,JDIM))*0.5
      ENDDO
C***
C***  INTERPOLATE VALUES TO THOSE CORRESPONDING TO H POINTS ON THE
C***  IDIM X JDIM GRID ON THE WEST AND EAST SIDES.
C***
      DO J=3,JDIM-2,2
        ARR2(1,J)=(ARR1(1,J+1)+ARR1(1,J-1))*0.5
        ARR2(IDIM,J)=(ARR1(IM,J+1)+ARR1(IM,J-1))*0.5
      ENDDO
C***
C***  INTERPOLATE TO CORRESPONDING H POINTS ON ODD ROWS
C***  EXCLUDING THE SIDES.
C***
      DO J=3,JDIM-2,2
      DO I=3,IDIM-2,2
        II=I/2
        ARR2(I,J)=0.25*(ARR1(II,J)+ARR1(II+1,J)+ARR1(II+1,J+1)+
     1                  ARR1(II+1,J-1))
      ENDDO
      ENDDO
C***
C***  INTERPOLATE TO CORRESPONDING H POINTS ON EVEN ROWS
C***  EXCLUDING THE SIDES.
C***
      DO J=2,JDIM-1,2
      DO I=2,IDIM-1,2
        II=I/2
        ARR2(I,J)=0.25*(ARR1(II,J)+ARR1(II+1,J)+ARR1(II,J+1)+
     1                  ARR1(II,J-1))
      ENDDO
      ENDDO
C***
C***  EXTRAPOLATE TO THE FOUR CORNER H POINTS.
C***
      ARR2(1,1)=0.5*(1.5*(ARR1(1,1)+ARR1(1,2))
     1              -0.5*(ARR1(2,1)+ARR1(1,4)))
      ARR2(IDIM,1)=0.5*(1.5*(ARR1(IM-1,1)+ARR1(IM,2))
     1              -0.5*(ARR1(IM-2,1)+ARR1(IM,4)))
      ARR2(1,JDIM)=0.5*(1.5*(ARR1(1,JDIM)+ARR1(1,JDIM-1))
     1              -0.5*(ARR1(2,JDIM)+ARR1(1,JDIM-3)))
      ARR2(IDIM,JDIM)=0.5*(1.5*(ARR1(IM-1,JDIM)+ARR1(IM,JDIM-1))
     1              -0.5*(ARR1(IM-2,JDIM)+ARR1(IM,JDIM-3)))
      RETURN
      END


      SUBROUTINE FILLV(ARR1,ARR2,IFLAG,IDIM,JDIM)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    FILLV       FILL V POINTS ON E-GRID H FIELD
C   PRGRMMR: BLACK           ORG: W/NP2      DATE: 92-12-23       
C     
C ABSTRACT:
C     THIS ROUTINE FILLS VELOCITY (V) POINTS ON A MASS
C     POINT E-GRID.  PASSED VARIABLE IFLAG TELLS THE ROUTINE
C     WHAT TO DO WITH NEGATIVE VALUES OBTAINED DURING THE
C     FILLING PROCESS
C        IFLAG=0, TAKE WHATEVER THE INTERPOLATION GIVES.
C        IFLAG=1, TRUNCATE AT ZERO.
C        IFLAG=2, TRUNCATE AT A VERY SMALL POSITIVE VALUE (EPS).
C     THE INTERPOLATION CAN USE A BILINEAR SCHEME OR CUBIC
C     INTERPOLATING POLYMONIALS.  PARAMETERS WGTCUB AND WGTLIN
C     CONTROL THIS OPTION.  SETTING WGTCUB TO 1.0 ACTIVATES
C     INTERPOLATION USING CUBIC INTERPOLATING POLYNOMIALS.
C     SETTING WGTCUB TO 0.0 ACTIVATES BILINEAR INTERPOLATION.
C
C   .     
C     
C PROGRAM HISTORY LOG:
C   ??-??-??  BLACK - ORIGINATOR
C   93-12-23  RUSS TREADON - ADDED COMMENTS AND VARIABLE
C                            GRID DIMENSIONS TO 
C                            SUBROUTINE CALL.
C   98-06-04  BLACK - CONVERSION TO 2-D
C     
C USAGE:    CALL FILLV(ARR1,ARR2,IFLAG,IDIM,JDIM)
C   INPUT ARGUMENT LIST:
C     ARR2     - VELOCITY POINT DATA ON STAGGERED E-GRID.
C     IFLAG    - 
C     IDIM     - FIRST DIMENSION OF FILLED E-GRID.
C     JDIM     - SECOND DIMENSION OF FILLED E-GRID.
C
C   OUTPUT ARGUMENT LIST: 
C     ARR1     - VELOCITY POINT DATA ON FILLED E-GRID.
C     
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C     UTILITIES:
C       EFILL - FILL MISSING VALUES ON STAGGERED
C               E-GRID WITH FIELD MEAN.
C       COEFI - COMPUTES ROW-WISE INTERPOLATION VALUES.
C       COEFJ - COMPUTES COLUMN-WISE INTERPOLATION VALUES.
C     LIBRARY:
C       COMMON
C     
C   ATTRIBUTES:
C     LANGUAGE: FORTRAN
C     MACHINE : CRAY C-90
C$$$  
C     
C----------------------------------------------------------------------
      PARAMETER (WGTCUB=0.0,WGTLIN=1.-WGTCUB,EPS=1.E-5)
      PARAMETER (P5C=.125,P15C=3.375,P25C=15.625
     1,          P35C=42.875,P45C=91.125)
      PARAMETER (P5S=.25,P15S=2.25,P25S=6.25
     1     ,          P35S=12.25,P45S=20.25)
C----------------------------------------------------------------------
                     D I M E N S I O N
     1 ARR1((IDIM+1)/2,JDIM),ARR2(IDIM,JDIM)
C     
C**********************************************************************
C     START FILLV HERE.
C     
C     REPLACE MISSING VALUES ON EGRID WITH FIELD MEAN
      IM=(IDIM+1)/2
      CALL EFILL(ARR1,IM,JDIM)
C
      DO J=1,JDIM
      DO I=1,IDIM
        ARR2(I,J)=0.
      ENDDO
      ENDDO
C***
C***  REASSIGN VALUES OF MODEL GRID TO IDIM X JDIM GRID ON ODD ROWS
C***
      DO J=1,JDIM,2
      DO I=1,IDIM,2
        II=(I+1)/2
        ARR2(I,J)=ARR1(II,J)
      ENDDO
      ENDDO
C***
C***  REASSIGN VALUES OF MODEL GRID TO IDIM X JDIM GRID ON EVEN ROWS
C***
      DO J=2,JDIM-1,2
      DO I=2,IDIM-1,2
        II=I/2
        ARR2(I,J)=ARR1(II,J)
      ENDDO
      ENDDO
C
      IF(WGTCUB.GT.0.99)GO TO 100
C***
C***  INTERPOLATE VALUES TO THOSE CORRESPONDING TO V POINTS NOW ON THE
C***  IDIM X JDIM GRID FOR THE BOTTOM AND TOP ROWS.
C***
      DO I=2,IDIM-1,2
        II=I/2
        ARR2(I,1)=(ARR1(II,1)+ARR1(II+1,1))*0.5
        ARR2(I,JDIM)=(ARR1(II,JDIM)+ARR1(II+1,JDIM))*0.5
      ENDDO
C***
C***  INTERPOLATE VALUES TO THOSE CORRESPONDING TO V POINTS NOW ON THE
C***  IDIM X JDIM GRID FOR WEST AND EAST SIDES.
C***
      DO J=2,JDIM-1,2
        ARR2(1,J)=(ARR1(1,J+1)+ARR1(1,J-1))*0.5
        ARR2(IDIM,J)=(ARR1(IM,J+1)+ARR1(IM,J-1))*0.5
      ENDDO
C***
C***  INTERPOLATE TO CORRESPONDING V POINTS ON ALL INNER POINTS
C***  OF ODD ROWS.
C***
      DO J=3,JDIM-2,2
      DO I=2,IDIM-1,2
        II=I/2
        ARR2(I,J)=0.25*(ARR1(II,J)+ARR1(II+1,J)
     1                 +ARR1(II,J+1)+ARR1(II,J-1))
      ENDDO
      ENDDO
C***
C***  INTERPOLATE TO CORRESPONDING V POINTS ON ALL INNER POINTS
C***  OF EVEN ROWS.
C***
      DO J=2,JDIM-1,2
      DO I=3,IDIM-2,2
        II=I/2
        ARR2(I,J)=0.25*(ARR1(II,J)+ARR1(II+1,J)
     1                 +ARR1(II+1,J+1)+ARR1(II+1,J-1))
      ENDDO
      ENDDO
      IF(WGTCUB.LT.0.01)GO TO 825
C***
C***  FILL IN THE VELOCITY POINTS USING CUBIC INTERPOLATING POLYNOMIALS.
C***  FIRST DO THE ROWS.
C***
  100 WGT1=0.5*WGTCUB+WGTLIN
      DO 200 J=1,JDIM,2
      DO I=3,IDIM-8,2
        CALL COEFI(ARR2,IDIM,JDIM,I,J,C1,C2,C3,C4)
        ARR2(I+3,J)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1               WGTLIN*ARR2(I+3,J))/WGT1
      ENDDO
C
      CALL COEFI(ARR2,IDIM,JDIM,1,J,C1,C2,C3,C4)
      ARR2(2,J)=(0.5*WGTCUB*(C1*P15C+C2*P15S+C3*1.5+C4)+
     1           WGTLIN*ARR2(2,J))/WGT1
      ARR2(4,J)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1           WGTLIN*ARR2(4,J))/WGT1
      CALL COEFI(ARR2,IDIM,JDIM,IDIM-6,J,C1,C2,C3,C4)
      ARR2(IDIM-3,J)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1                WGTLIN*ARR2(IDIM-3,J))/WGT1
      ARR2(IDIM-1,J)=(0.5*WGTCUB*(C1*P35C+C2*P35S+C3*3.5+C4)+
     1                WGTLIN*ARR2(IDIM-1,J))/WGT1
  200 CONTINUE
C
      DO 400 J=2,JDIMM,2
      DO I=4,IDIM-9,2
        CALL COEFI(ARR2,IDIM,JDIM,I,J,C1,C2,C3,C4)
        ARR2(I+3,J)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1               WGTLIN*ARR2(I+3,J))/WGT1
      ENDDO
C
      CALL COEFI(ARR2,IDIM,JDIM,2,J,C1,C2,C3,C4)
      ARR2(1,J)=(0.5*WGTCUB*(C1*P5C+C2*P5S+C3*.5+C4)+
     1           WGTLIN*ARR2(1,J))/WGT1
      ARR2(3,J)=(0.5*WGTCUB*(C1*P15C+C2*P15S+C3*1.5+C4)+
     1           WGTLIN*ARR2(3,J))/WGT1
      ARR2(5,J)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1           WGTLIN*ARR2(5,J))/WGT1
      CALL COEFI(ARR2,IDIM,JDIM,IDIM-7,J,C1,C2,C3,C4)
      ARR2(IDIM-4,J)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1                WGTLIN*ARR2(IDIM-4,J))/WGT1
      ARR2(IDIM-2,J)=(0.5*WGTCUB*(C1*P35C+C2*P35S+C3*3.5+C4)+
     1                WGTLIN*ARR2(IDIM-2,J))/WGT1
      ARR2(IDIM,J)=(0.5*WGTCUB*(C1*P45C+C2*P45S+C3*4.5+C4)+
     1                WGTLIN*ARR2(IDIM,J))/WGT1
  400 CONTINUE
C***
C***  NOW DO THE COLUMNS.
C***
      WGT2=0.5*WGTCUB+WGT1
      DO 600 I=1,IDIM,2
      DO J=3,JDIM-8,2
        CALL COEFJ(ARR2,IDIM,JDIM,I,J,C1,C2,C3,C4)
        ARR2(I,J+3)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1               WGT1*ARR2(I,J+3))/WGT2
      ENDDO
C
      CALL COEFJ(ARR2,IDIM,JDIM,I,1,C1,C2,C3,C4)
      ARR2(I,2)=(0.5*WGTCUB*(C1*P15C+C2*P15S+C3*1.5+C4)+
     1           WGT1*ARR2(I,2))/WGT2
      ARR2(I,4)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1           WGT1*ARR2(I,4))/WGT2
      CALL COEFJ(ARR2,IDIM,JDIM,I,JDIM-6,C1,C2,C3,C4)
      ARR2(I,JDIM-3)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1                WGT1*ARR2(I,JDIM-3))/WGT2
      ARR2(I,JDIM-1)=(0.5*WGTCUB*(C1*P35C+C2*P35S+C3*3.5+C4)+
     1                WGT1*ARR2(I,JDIM-1))/WGT2
  600 CONTINUE
C
      DO 800 I=2,IDIM-1,2
      DO J=4,JDIM-9,2
        CALL COEFJ(ARR2,IDIM,JDIM,I,J,C1,C2,C3,C4)
        ARR2(I,J+3)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1               WGT1*ARR2(I,J+3))/WGT2
      ENDDO
C
      CALL COEFJ(ARR2,IDIM,JDIM,I,2,C1,C2,C3,C4)
      ARR2(I,1)=(0.5*WGTCUB*(C1*P5C+C2*P5S+C3*.5+C4)+
     1           WGT1*ARR2(I,1))/WGT2
      ARR2(I,3)=(0.5*WGTCUB*(C1*P15C+C2*P15S+C3*1.5+C4)+
     1           WGT1*ARR2(I,3))/WGT2
      ARR2(I,5)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1           WGT1*ARR2(I,5))/WGT2
      CALL COEFJ(ARR2,IDIM,JDIM,I,JDIM-7,C1,C2,C3,C4)
      ARR2(I,JDIM-4)=(0.5*WGTCUB*(C1*P25C+C2*P25S+C3*2.5+C4)+
     1                WGT1*ARR2(I,JDIM-4))/WGT2
      ARR2(I,JDIM-2)=(0.5*WGTCUB*(C1*P35C+C2*P35S+C3*3.5+C4)+
     1                WGT1*ARR2(I,JDIM-2))/WGT2
      ARR2(I,JDIM)=(0.5*WGTCUB*(C1*P35C+C2*P35S+C3*3.5+C4)+
     1                WGT1*ARR2(I,JDIM))/WGT2
  800 CONTINUE
C
  825 IF(IFLAG.EQ.1)THEN
        DO J=1,JDIM,2
        DO I=2,IDIM-1,2
          IF(ARR2(I,J).LT.0.)ARR2(I,J)=0.
        ENDDO
        ENDDO
C
        DO J=2,JDIM-1,2
        DO I=1,IDIM,2
          IF(ARR2(I,J).LT.0.)ARR2(I,J)=0.
        ENDDO
        ENDDO
      ELSEIF(IFLAG.EQ.2)THEN
        DO J=1,JDIM,2
        DO I=2,IDIM-1,2
          IF(ARR2(I,J).LT.0.)ARR2(I,J)=EPS
        ENDDO
        ENDDO
C
        DO J=2,JDIM-1,2
        DO I=1,IDIM,2
          IF(ARR2(I,J).LT.0.)ARR2(I,J)=EPS
        ENDDO
        ENDDO
      ENDIF
C     
C     END OF ROUTINE.
C
      RETURN
      END
