CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     This file corresponds to the slgcm.misc for the 3D model    C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE AVRX3(U,V)
      include "slgcm.slgprm"
C
C*********************************************************************
C*********************** SLGCM   MODEL *******************************
C********************* SUBROUTINE AVRX *******************************
C********************** 08 OCTO    1992 ******************************
C*********************************************************************
C
C
      PARAMETER (NLATS=JM/8)
C
      INTEGER   TOTLATS, TEMPSIZE, FORWARD, BACKWARD
C
      PARAMETER (TOTLATS=4*NLATS)
      PARAMETER (TEMPSIZE=TOTLATS*IMP2, ITRG=3*IMP2+1)
      PARAMETER (FORWARD=-1,BACKWARD=1)
C
C  ARGUMENTS
C
      DIMENSION    U(IMP3,JM), V(IMP3,JM)
     *,            S(IMP2,TOTLATS), TR(ITRG), IX(IM)
C
C   LOCALS
C
      LOGICAL        FIRST
C
C   SCRATCH SPACE ALLOCATION
C
      DIMENSION      A(IMP2,TOTLATS), B(IM*2,TOTLATS)
      DIMENSION      F(IMP2,JMM1),    Z(IMP3,TOTLATS)
C
      DATA           FIRST/.TRUE./
      DATA S/TEMPSIZE*0.0/, TR/ITRG*0.0/, IX/IM*0/
      SAVE
C
C   INITIALIZATION ON FIRST CALL
C
      IF (FIRST) THEN
         FIRST = .FALSE.
C
         DO 10 I=1,IMP2*JMM1
         F(I,1) = 1.
   10    CONTINUE
         CONS   = ONE / COS(-PI/TWO+(NLATS+1)*(PI/JMM1))
C
         DO 100 I=4,IM
         ANGI = SIN((I/2)*PI/IM)
         DO 100 J=2,JMM1
         TEM      = CONS * COS(-PI/TWO + (J-1)*(PI/JMM1))
         F(I+1,J) = MIN(TEM*TEM/ANGI,ONE)
100      CONTINUE
C
         DO 200 I=1,IMP2
         DO 200 J=1,NLATS
         S(I,J      ) = F(I,         1 +J)
         S(I,J+NLATS) = F(I,JM-NLATS-1 +J)
         S(I,J+2*NLATS) = S(I,J)
         S(I,J+3*NLATS) = S(I,J+NLATS)
200      CONTINUE
C
         CALL FFTFAX(IM, IX, TR)
C
      ENDIF
C
C       END INITIALIZATION
C**********************************************************************
C       BEGIN STANDARD CASE
C
C   COPY LATITUDES TO BE SMOOTHED INTO CONTIGUOUS BUFFER
C
      IOFF = 1
      DO 1003 I=1,IMP3*NLATS
      Z(I,IOFF)         = U(I,2)
      Z(I,NLATS+IOFF)   = U(I,JM-NLATS)
      Z(I,2*NLATS+IOFF) = V(I,2)
      Z(I,3*NLATS+IOFF) = V(I,JM-NLATS)
1003  CONTINUE
C
C   TRANSFORM AND FILTER
C
      DO 1100 J=1,TOTLATS
      DO 1100 I=1,IMP2
      A(I,J) = Z(I,J)
1100  CONTINUE
C
c     CALL RFFTMLT(A,B,TR,IX,1,IMP2,IM,TOTLATS,FORWARD)
      CALL FFT99(A,B,TR,IX,1,IMP2,IM,TOTLATS,FORWARD)
C
      DO 1110 I=1,IMP2*TOTLATS
      A(I,1) = A(I,1)*S(I,1)
 1110 CONTINUE
C
c     CALL RFFTMLT(A,B,TR,IX,1,IMP2,IM,TOTLATS,BACKWARD)
      CALL FFT99(A,B,TR,IX,1,IMP2,IM,TOTLATS,BACKWARD)
C
      DO 1500 J=1,TOTLATS
      DO 1450 I=1,IM
      Z(I,J) = A(I,J)
1450  CONTINUE
      Z(IMP1,J) = Z(1,J)
      Z(IMP2,J) = Z(2,J)
      Z(IMP3,J) = Z(3,J)
1500  CONTINUE
C
C   COPY BACK FROM CONTIGUOUS BUFFER
C
      IOFF = 1
      DO 1645 I=1,IMP3*NLATS
      U(I,2)        = Z(I,IOFF)
      U(I,JM-NLATS) = Z(I,NLATS+IOFF)
      V(I,2)        = Z(I,2*NLATS+IOFF)
      V(I,JM-NLATS) = Z(I,3*NLATS+IOFF)
1645  CONTINUE
C
      RETURN
      END
      SUBROUTINE CONSL

      include "slgcm.slgprm"
      include "slgcm.slgclm"
      include "slgcm.slgcon"
      real wk1(ijm0m0)
C
      DO 220 J=1,JM
      PH(J,2) = -PIO2  + (J-1)*DPH
      PH(J,1) = PH(J,2) + 0.5*DPH
  220 CONTINUE
C
      DO 225 K=1,2
      DO 225 J=1,JM
      SPH(J,K) = SIN(PH(J,K))
      CPH(J,K) = COS(PH(J,K))
      DX(J,K) = AE*DLM*CPH(J,K)
  225 CONTINUE
      CPH(1,2)  = 0.0
      CPH(JM,2) = 0.0
      DX(1,2)   = 0.0
      DX(JM,2)  = 0.0
C
      DO 226 J=1,JMM1
      RH1JK(J,1) = ONE / (AE*CPH(J,1))
  226 CONTINUE
      DO 227 J=2,JMM1
      RH1JK(J,2) = ONE / (AE*CPH(J,2))
  227 CONTINUE
      RH1JK(1,2)  = 0.0
      RH1JK(JM,1) = 0.0
      RH1JK(JM,2) = 0.0
C
      DO 5010 J=1,JM
      II = (J-1)*IMP3
      DO 5009 I=1,IMP3
      PHV(II+I) =  PH(J,1)
      PHU(II+I) =  PH(J,2)
      SPHV(II+I) = SPH(J,1)
      SPHU(II+I) = SPH(J,2)
      CPHV(II+I) = CPH(J,1)
      CPHU(II+I) = CPH(J,2)
      RH1V(II+I) = RH1JK(J,1)
      RH1U(II+I) = RH1JK(J,2)
      WK1(II+I) =  DX(J,2)
 5009 CONTINUE
 5010 CONTINUE
C
      DO 215 I=1,IJM2M0
      TEM       = 1.0 / ((DY*DY)  * CPHU(I1+I))
      CD1(I1+I) = 1.0 / (WK1(I1+I)*WK1(I1+I))
      CD2(I1+I) = TEM * CPHV(I1+I)
      CD3(I1+I) = TEM * CPHV(I)
  215 CONTINUE
C
C  K=1 GIVES THE ANGLE LAMDA AT THE U-POINTS
C  K=2 GIVES THE ANGLE LAMDA AT THE H-POINTS AND THE V-POINTS
C
      DO 221 I=1,IM
      BLMIK(I,2) = (I-1)*DLM
      BLMIK(I,1) = BLMIK(I,2) +  0.5*DLM
  221 CONTINUE
      DO 222 I=1,3
      BLMIK(IM+I,1) = BLMIK(I,1)
      BLMIK(IM+I,2) = BLMIK(I,2)
  222 CONTINUE
      DO 223 K=1,2
      DO 223 I=1,IM
      CLM(I,K) = COS(BLMIK(I,K))
      SLM(I,K) = SIN(BLMIK(I,K))

      CLM(I,K+2) = COS(2.0*BLMIK(I,K))
      SLM(I,K+2) = SIN(2.0*BLMIK(I,K))

      CLM(I,K+4) = COS(3.0*BLMIK(I,K))
      SLM(I,K+4) = SIN(3.0*BLMIK(I,K))
  223 CONTINUE
C
C
      DO 230 J=2,JMM1
      CDIVU(J)  = ONE  / DX(J,2)
      CDIVV(J)  = ONE  / (DY*CPH(J,2))
  230 CONTINUE
C
C
      AREA   = 0.0
      DO 2005 J=2,JMM1
      DO 2005 I=1,IM
      AREA   = AREA + (DX(J,2)*DY)
 2005 CONTINUE
      ARE2   = AREA + IM * DY * DX(1,1) * 0.5
      AREA   = AREA + 2.0 * PI * DY*DY*0.25
C
      DO 2007 J=2,JMM1
      II = (J-1)*IMP3
      DO 2007 I=1,IMP3
      CDU(II+I) = CDIVU(J)
      CDV(II+I) = CDIVV(J)
 2007 CONTINUE
C
      DO 8888 J=1,JM
      II = (J-1)*IMP3
      DO 8888 I=1,IMP3
      xcr(II+I,1) = (I - 1.0) * dlm
      xcr(II+I,2) = (I - 0.5) * dlm
 8888 continue

      RETURN
      END

      SUBROUTINE GLMNP2(IMP3, JM, SRP, DX, DY, ARE2, DXV, GLBM)
C
      DIMENSION SRP(IMP3,JM), DX(JM)
      double precision dglbm
C
      IM   = IMP3 - 3
      JMM1 = JM - 1
      DYBA = DY / ARE2
C
      DGLBM = (SRP(1,1) + SRP(1,JM)) * (DYBA*DXV*0.25*IM)
      DO 3005 J=2,JMM1
      DO 3005 I=1,IM
      DGLBM = DGLBM + SRP(I,J) * (DX(J)*DYBA)
 3005 CONTINUE
      GLBM  = DGLBM
C
      RETURN
      END

      SUBROUTINE INTER(IEM, JEM, IMP3, JM, FEC, FSL, WRK)
C
      DIMENSION FEC(IEM,JEM), FSL(IMP3,JM)
      DIMENSION WRK(IMP3,JEM)
C
      IM = IMP3 - 3
      JMM1 = JM - 1
C
      if (iem-1 .eq. im .and. jem .eq. jmm1) then
         do 5 j=2,jmm1
         do 5 i=1,im
         fsl(i,j) = fec(i,j-1)
    5    continue
         do 6 i=1,im
         wrk(i,1)     = fec(i,1)
         wrk(i,jem-1) = fec(i,jem-1)
    6    continue
      else
C
         DLON = 360.0 / FLOAT(IM)
         DLAT = 180.0 / FLOAT(JMM1)
         DLOE = FLOAT(IEM-1) / FLOAT(IM)
         DLAE = FLOAT(JEM)   / FLOAT(JMM1)
      DO J1X = 1, JEM
         FEC(1,J1X) = FEC(IEM,J1X)
      END DO

c
         DO 20 J=1,JEM
         DO 10 I=1,IM
         XLON = (I-1) * DLOE + 1
         ILON = XLON
         XLON = XLON - ILON
         WRK(I,J) = (1.0-XLON)*FEC(ILON,J) + XLON*FEC(ILON+1,J)
   10    CONTINUE
   20    CONTINUE
C
         DO 50 J=2,JMM1
         YLAT =  (J-1) * DLAE
         JLAT = YLAT
         YLAT = YLAT - JLAT
         IF (JLAT .GT. 0 .AND. JLAT .LT. JEM-1) THEN
            DO 40 I=1,IM
            FSL(I,J) = (1.0-YLAT)*WRK(I,JLAT) + YLAT*WRK(I,JLAT+1)
   40       CONTINUE
         ELSEIF (JLAT .LE. 0) THEN
            DO 42 I=1,IM
            FSL(I,J) = WRK(I,1)
   42       CONTINUE
         ELSE
            DO 43 I=1,IM
            FSL(I,J) = WRK(I,JEM-1)
   43       CONTINUE
         ENDIF
C
   50    CONTINUE
      endif
C
      TEMS = 0.0
      TEMN = 0.0
      DO 60 I=1,IM
      TEMS = TEMS + WRK(I,1)
      TEMN = TEMN + WRK(I,JEM-1)
   60 CONTINUE
      TEMS = TEMS / FLOAT(IM)
      TEMN = TEMN / FLOAT(IM)
C
      DO 65 I=1,IM
      FSL(I,1)  = TEMS
      FSL(I,JM) = TEMN
   65 CONTINUE
C
      CALL WRAP (1,IMP3*JM,FSL(1,1))
C
      RETURN
      END
      SUBROUTINE WRAP(IB,IE,XX)
      include "slgcm.slgprm"
      DIMENSION XX(IJM0M0)
C
      DO 10 I=IB,IE,IMP3
      XX(IM+I)   = XX(I)
      XX(IMP1+I) = XX(I+1)
      XX(IMP2+I) = XX(I+2)
   10 CONTINUE
C
      RETURN
      END
      SUBROUTINE CONDT(DT)
C
      include "slgcm.slgprm"
      include "slgcm.slgclm"
      include "slgcm.slgcon"
C
      INTEGER DT
C
C     SAVE
C
      CDA = CDF * CDAMP
      CDB = (2.0 - CDF) * CDAMP
C
      PRINT *,' CDAMP=',CDAMP, ' CDA=',CDA,' CDB=',CDB
C
      DTOA   = DT/AE
      DTHF   = DT*0.5
      DTHFSQ = DTHF*DTHF
      DL1DT  = DL1 * DT
      DL2DT  = DL2 * DT
      HFCDDT = CDA*DL1DT
C
C PRECOMPUTE ARRAYS AT V POINTS
C
        DO 1226 J=1,JM
C
        FCOR(J,1)  = TWOM * SPH(J,1)
        FDTHF(J,1) = FCOR(J,1) * DL1DT
        FVSQ(J)    = FDTHF(J,1) * FDTHF(J,1)
        ARR1(J,1)  = ONE / (ONE + FVSQ(J))
        ARR2(J,1)  = FDTHF(J,1) / (ONE + FVSQ(J))
        ARR5(J)    = DL1DT * ARR2(J,1) * RH1JK(J,1)
        ARR6(J)    = DL1DT * ARR1(J,1) * RA
        F1M(J,1)   =(ONE - FVSQ(J))   * ARR1(J,1)
C
C PRECOMPUTE ARRAYS AT U POINTS
C
        FCOR(J,2)  = TWOM*SPH(J,2)
        FDTHF(J,2) = FCOR(J,2)  * DL1DT
        FUSQ(J)    = FDTHF(J,2) * FDTHF(J,2)
        ARR1(J,2)  = ONE / (ONE + FUSQ(J))
        ARR2(J,2)  = FDTHF(J,2) / (ONE + FUSQ(J))
        ARR3(J)    = DL1DT * ARR1(J,2) * RH1JK(J,2)
        ARR4(J)    = DL1DT * ARR2(J,2) * RA
        F1M(J,2)   = (ONE - FUSQ(J))  * ARR1(J,2)
C
 1226 CONTINUE
C
C
      TEM  = DPH * DPH * DL1DT * TWOM * 0.5 / (DLM*DY*DY)
      AESQ = AE * AE
      DO 240 J=2,JMM1
      TEM1 = 1.0 / (DX(J,2) * DX(J,2))
      TEM2 = TEM * (1.0 - FUSQ(J)) * ARR1(J,2)
      TEM5 = ONE / (ARR1(J,2) * DY * DY *CPH(J,2))
      TEM3 = ARR1(J,1)   * CPH(J,1)   * TEM5
      TEM4 = ARR1(J-1,1) * CPH(J-1,1) * TEM5
      II = (J-1) * IMP3
      DO 235 I=1,IMP3
      CDW(II+I) = TEM1 + TEM2
      CDE(II+I) = TEM1 - TEM2
      CDN(II+I) = TEM3
      CDS(II+I) = TEM4
C
  235 CONTINUE
C
      CDN2(J) = ARR1(J,1) * CPH(J,1) / (CPH(J,2) *  DY*DY)
      CDS2(J) = ARR1(J-1,1) * CPH(J-1,1) / (CPH(J,2) * DY*DY)
      CDC1(J) = ARR1(J,2) / (CPH(J,2) * CPH(J,2) *AESQ)
      CDC2(J) = ARR1(J,2)*ARR1(J,2) *(1.0-FUSQ(J))*(DL1DT*TWOM) / AESQ
C
  240 CONTINUE
      CDC1(1)  = 4.0 * ARR1(1,1) / (DY*DY)
      CDC1(JM) = 4.0 * ARR1(JMM1,1) / (DY*DY)
C
c     PRINT *,' ARR1=',ARR1
C     PRINT *,' CDW=',CDW
C     PRINT *,' CDE=',CDE
C     PRINT *,' CDN=',CDN
C     PRINT *,' CDS=',CDS
c     PRINT *,' CDN2=',CDN2
c     PRINT *,' CDS2=',CDS2
c     PRINT *,' CDC1=',CDC1
c     PRINT *,' CDC2=',CDC2
C
C     
      SQMU = 1.0/(DL1DT*DL1DT*FI0+CDA*DL1DT)
C
      PRINT *,' SQMU=',SQMU
C
      RETURN
      END
      SUBROUTINE ENGCH(TKEU,TKEV,TKE,POL,EG,GMASS)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     This subroutine calculates the global integral of energy   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      include "slgcm.slgprm"
      include "slgcm.slgprg"
      include "slgcm.slgcon"

      DIMENSION AU(JM), AV(JM), AG(JM)

      DO J = 2, JMM1
         AU(J) = (AE**2)*DPH*DLM*CPH(J,2)
      ENDDO
      
      AU(1) = 0.25*(AE**2)*(DPH**2)*PI
      AU(JM) = AU(1)

      DO J1X = 1, JM
         AG(J1X) = AU(J1X)
      END DO

      DO J = 1, JMM1
         AV(J) = (AE**2)*DPH*DLM*CPH(J,1)
      ENDDO

      TKEU = 0.0
      TKEV = 0.0
      POL = 0.0
      EG = 0.0
      GMASS = 0.0

      AAA0 = 6.0*PI*(AE**2)

      DO JJ = 2, JMM1
         DO II = 1, IM
            I = (JJ-1)*IMP3 + II
            TKEU = TKEU + 0.5*(UN0(I)**2)*AU(JJ)*FI0
            POL = POL + 0.5*(GTP(I)**2)*AG(JJ)/9.8
            GMASS = GMASS + GTP(I)*AG(JJ)
         ENDDO
      ENDDO

      TKEU = TKEU + 0.5*(UN0(1)**2)*AU(1)*FI0
     #     + 0.5*(UN0(IJM0M0)**2)*AU(JM)*FI0
      POL = POL + 0.5*(GTP(1)**2)*AG(1)/9.8
     #          + 0.5*(GTP(IJM0M0)**2)*AG(JM)/9.8
      GMASS = GMASS + GTP(1)*AG(1) + GTP(IJM0M0)*AG(JM)

      GMASS = GMASS/(AAA0*9.81)

      DO JJ = 1, JMM1
         DO II = 1, IM
            I = (JJ-1)*IMP3 + II
            TKEV = TKEV + 0.5*(VN0(I)**2)*AV(JJ)*FI0
         ENDDO
      ENDDO

      TKE = (TKEU + TKEV)/9.8
      EG = TKE + POL

      RETURN
      END

C***********************************************************************        
C                                                                               
C  SUBROUTINE SPH2SPH                                                           
C                                                                               
C  PURPOSE                                                                      
C     INTERPOLATE FROM ONE SPHERICAL COORDINATE SYSTEM TO                       
C     ANOTHER SPHERICAL COORDINATE SYSTEM                                       
C                                                                               
C  USAGE                                                                        
C     CALL SPH2SPH(AA,IAM,JAM,ANP,AEQ,BB,IBM,JBM,BNP,BEQ,EQB1)                  
C                                                                               
C  DESCRIPTION OF PARAMETERS                                                    
C     AA      - TWO DIMENSIONAL INPUT ARRAY                                     
C     IAM     - 1ST DIMENSION OF AA (FOR AA DLAM = 2*PI / IAM)                  
C     JAM     - 2ND DIMENSION OF AA                                             
C     ANP     - PSEUDO J-INDEX OF NORTH POLE FOR ARRAY AA                       
C     AEQ     - PSEUDO J-INDEX OF EQUATOR FOR ARRAY AA                          
C     BB      - TWO DIMENSIONAL OUTPUT ARRAY                                    
C     IBM     - 1ST DIMENSION OF BB (FOR BB DLAM = 2*PI / IBM)                  
C     JBM     - 2ND DIMENSION OF BB                                             
C     BNP     - PSEUDO J-INDEX OF NORTH POLE FOR ARRAY BB                       
C     BEQ     - PSEUDO J-INDEX OF EQUATOR FOR ARRAY BB                          
C     EQB1    - IN ARRAY AA A PSEUDO I-INDEX OF EQB1                            
C               CORRESPONDS TO AN I-INDEX OF 1 IN ARRAY BB                      
C                                                                               
C  REMARKS                                                                      
C     (1) THE VALUES AT THE NORTH AND SOUTH EDGES OF THE AA GRID                
C         ARE ASSUMED TO EXTEND TO THE POLES                                    
C     (2) THE INTERPOLATION IS VIA AREA WEIGHTED AVERAGES                       
C     (3) THE INPUT PARAMETERS MUST SATISFY THE INEQUALITIES                    
C         ANP.LT.AEQ,  BNP.LT.BEQ,  AND  0.0.LE.EQB1.LT.IAM                     
C                                                                               
C  SUBPROGRAMS REQUIRED                                                         
C     NONE                                                                      
C                                                                               
C  METHOD                                                                       
C     AREA WEIGHTED MEANS                                                       
C                                                                               
C***********************************************************************        
C                                                                               
      SUBROUTINE SPH2SPH(AA,IAM,JAM,ANP,AEQ,BB,IBM,JBM,BNP,BEQ,EQB1) 
      DIMENSION     AA(IAM,JAM),BB(IBM,JBM)                                     
      DIMENSION     QBW(181)                                                    
      PARAMETER ( AHALF = 0.5 )                                                 
      PARAMETER ( ONE   = 1.0 )                                                 
      PARAMETER ( FOUR  = 4.0 )                                                 
C                                                                               
C***********************************************************************        
                                                                                
      write(6,1000) iam,jam,ibm,jbm                                             
 1000 format(1x,'Using SPH2SPH interpolation from (',i3,',',i3,')',             
     .          ' to (',i3,',',i3,')')                                          
                                                                                
C                                                                               
C  SET CONSTANTS                                                                
C  =============                                                                
      PI         = FOUR * ATAN(ONE)                                             
      BBSTART    = MAX (AHALF,BNP)                                              
      BBSTOP     = MIN (JBM+AHALF,BEQ+BEQ-BNP)                                  
      BDIF       = BEQ - BNP                                                    
      ADIF       = AEQ - ANP                                                    
      QBDA       = BDIF / ADIF                                                  
      QADB       = ADIF / BDIF                                                  
      AASTART    = ANP + QADB * (BBSTART - BNP)                                 
      AASTOP     = ANP + QADB * (BBSTOP  - BNP)                                 
      QBA        = (ONE * IBM) / IAM                                            
      SAA        = IAM - EQB1                                                   
      QPIBDIF    = AHALF * PI / BDIF                                            
      ANGU       = QPIBDIF * (BBSTART - BEQ)                                    
      YU         = BBSTART                                                      
      DO 100 JB  = 1,JBM                                                        
      YO         = YU                                                           
      YU         = MIN (JB+AHALF,BBSTOP)                                        
      ANGO       = ANGU                                                         
      ANGU       = QPIBDIF * (YU - BEQ)                                         
  100 QBW(JB)    = ONE / (SIN(ANGO) - SIN(ANGU))                                
C                                                                               
C  ZERO-OUT ARRAY BB                                                            
C  =================                                                            
      DO 110 JB  = 1,JBM                                                        
      DO 110 IB  = 1,IBM                                                        
  110 BB(IB,JB)  = 0.0                                                          
C                                                                               
C  START THE J-LOOP FOR THE AA DATA RECTANGLES                                  
C  ===========================================                                  
      AAUNDER    = AASTART                                                      
      BBUNDER    = BBSTART                                                      
      JASTART    = MAX (1,INT(AASTART+AHALF))                                   
      JASTART    = MIN (JASTART,JAM)                                            
      DO 300 JA  = JASTART,JAM                                                  
C                                                                               
C  COMPUTE UPPER AND LOWER LIMITS OF AA DATA RECTANGLES                         
C  ====================================================                         
      BBOVER     = BBUNDER                                                      
      AAUNDER    = MIN (AASTOP,JA+AHALF)                                        
      IF(JA.EQ.JAM) AAUNDER = AASTOP                                            
      BBUNDER    = BNP + QBDA * (AAUNDER - ANP)                                 
C                                                                               
C  COMPUTE JB AT TOP OF DATA RECTANGLE                                          
C  ===================================                                          
      JBSTART    = INT(BBOVER+AHALF)                                            
C                                                                               
C  START THE I-LOOP FOR THE AA DATA RECTANGLES                                  
C  ===========================================                                  
      AARIGHT    = AHALF                                                        
      BBRIGHT    = ONE + QBA * (AARIGHT + SAA)                                  
      DO 200 IA  = 1,IAM                                                        
      AAIAJA     = AA(IA,JA)                                                    
C                                                                               
C  COMPUTE LEFT AND RIGHT LIMITS OF AA DATA RECTANGLE                           
C  ==================================================                           
      BBLEFT     = BBRIGHT                                                      
      AARIGHT    = IA + AHALF                                                   
      BBRIGHT    = ONE + QBA * (AARIGHT + SAA)                                  
C                                                                               
C  COMPUTE IB AT LEFT EDGE OF AA DATA RECTANGLE                                 
C  ============================================                                 
      IBSTART    = INT(BBLEFT+AHALF)                                            
C                                                                               
C  START THE J-LOOP FOR BB COVERAGE                                             
C  ================================                                             
      JB         = JBSTART - 1                                                  
      YO         = BBOVER                                                       
   20 JB         = JB + 1                                                       
      JBB        = MIN (JB,JBM)                                                 
      JBB        = MAX (JBB,1)                                                  
      YU         = MIN (JB+AHALF,BBUNDER)                                       
C                                                                               
C  COMPUTE WEIGHT FACTOR FOR A FULL DLAM BB RECTANGLE                           
C  ==================================================                           
      WDX        = SIN(QPIBDIF*(YO-BEQ)) - SIN(QPIBDIF*(YU-BEQ))                
      WDX        = WDX * QBW(JBB)                                               
C                                                                               
C  START THE I-LOOP FOR BB COVERAGE                                             
C  ================================                                             
      IB         = IBSTART - 1                                                  
      XL         = BBLEFT                                                       
   30 IB         = IB + 1                                                       
      IBB        = MOD(IB-1,IBM) + 1                                            
      XR         = MIN(IB+AHALF,BBRIGHT)                                        
C                                                                               
C  COMPUTE WEIGHT FACTOR                                                        
C  =====================                                                        
      WEIGHT     = WDX * (XR - XL)                                              
C                                                                               
C  ADD CONTRIBUTION TO BB(IBB,JBB)                                              
C  ===============================                                              
      BB(IBB,JBB)= BB(IBB,JBB) + WEIGHT * AAIAJA                                
C                                                                               
C  CONTINUE THE I-LOOP FOR BB COVERAGE                                          
C  ===================================                                          
      IF(XR.GE.BBRIGHT) GO TO 31                                                
      XL         = XR                                                           
      GO TO 30                                                                  
C                                                                               
C  CONTINUE THE J-LOOP FOR BB COVERAGE                                          
C  ===================================                                          
   31 IF(YU.GE.BBUNDER) GO TO 21                                                
      YO         = YU                                                           
      GO TO 20                                                                  
   21 CONTINUE                                                                  
C                                                                               
C  END THE I-LOOP FOR THE AA RECTANGLES                                         
C  ====================================                                         
  200 CONTINUE                                                                  
C                                                                               
C  END THE J-LOOP FOR THE AA RECTANGLES                                         
C  ====================================                                         
      IF(AAUNDER.EQ.AASTOP) R E T U R N                                         
  300 CONTINUE                                                                  
C                                                                               
      E N D                                                                     


