      SUBROUTINE SLGCST(DT, FIRST)

      include "slgcm.slgprm"
      include "slgcm.slgcon"
      include "slgcm.slgprg"
      include "slgcm.slgclm"
      include "slgcm.slgwrk"

      LOGICAL FIRST

      DIMENSION SLUH(IJM0M0), SLVH(IJM0M0), SPUH(IJM0M0), SPVH(IJM0M0)
     *,         SLAM(IJM0M0), SPHI(IJM0M0)
     *,         CBV(IJM0M0),  SBV(IJM0M0)
     *,         sluh1(imp3,jm),slvh1(imp3,jm), spuh1(imp3,jm)
     *,         spvh1(imp3,jm),slam1(imp3,jm),sphi1(imp3,jm)
      DIMENSION    BLS(IJM0M0),  BTS(IJM0M0)
     *,         UNSS(IJM0M0)

     *,         UST(IJP3M0),   VST(IJP3M0)
     *,         ust1(imp3,jmp3),vst1(imp3,jmp3)  
     *,       YUUPT(IJP3M0), YVVPT(IJP3M0)
     *,          BL(IJP3M0), PPNH1(IJP3M0)

      DIMENSION TX1(IJP3M0),   TX2(IJP3M0),   TX3(IJP3M0)

      DIMENSION UTS(IJM0M0),    VTS(IJM0M0)
     *,         UMH(IJM0M0),    VMH(IJM0M0)
     *,       AUUPT(IJM0M0),  AVVPT(IJM0M0)
     *,      YUUPTS(IJM0M0), YUVPTS(IJM0M0)
     *,      YVUPTS(IJM0M0), YVVPTS(IJM0M0)

      DIMENSION PPN(IJM0M0), PPNH(IJM0M0), PPNO(IJM0M0), WF1(IJM0M0)


      DIMENSION  C1A(IJM0M0), C2A(IJM0M0)
     *,          c1a1(imp3,jm),c2a1(imp3,jm)

      DIMENSION  RHS(IJM0M0,2), PHI(IJM0M0,2), WSV(100*IM), IX(IM)
      DIMENSION  WKTMP(16*IJM0M0)

      LOGICAL INIT1, INIT2
      DATA INIT1/.TRUE./, INIT2/.TRUE./
      DATA NSTEP/0/
      data ia/158/
      SAVE

c     The following is to test the results
c	open (unit=15, FILE='bill1.out',FORM='FORMATTED',
c     +        ACCESS='SEQUENTIAL', STATUS='NEW')
c	open (unit=16, FILE='bill2.out',FORM='FORMATTED',
c     +        ACCESS='SEQUENTIAL', STATUS='NEW')
c	open (unit=17, FILE='bill3.out',FORM='FORMATTED',
c     +        ACCESS='SEQUENTIAL', STATUS='NEW')
      
      NSTEP = NSTEP + 1

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     STEP 1: ----continued                                         C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      CALL WRAP(I1P1,IJM1M0, UN0(1))
      CALL WRAP(1,   IJM1M0, VN0(1))
      CALL WRAP(1,   IJM0M0, GTP(1))

      CALL FILPU(UN0(1))
      CALL FILPV(VN0(1))
C     
C     Calculate the initial divergence     
C
      IF(FIRST) THEN
         DO 21 I=2,IJM2M0
            WK1(I1+I) = UN0(I1+I) - UN0(I1+I-1)
 21      CONTINUE
         DO 25 I=I1P1,IJM1M0,IMP3
            WK1(I) = WK1(IM+I)
 25      CONTINUE
         DO 30 I=1,IJM1M0
            WK5(I) = VN0(I) * CPHV(I)
 30      CONTINUE
         DO 35 I=I1P1,IJM1M0
            DIV(I) = WK1(I) * CDU(I) + (WK5(I) - WK5(I-I1)) * CDV(I)
 35      CONTINUE

C
C     Calculate the dibergence at the poles
C
         DNPOLE = 0.0
         DSPOLE = 0.0
         DO 36 I=1,IM
            DNPOLE = DNPOLE + VN0(IJM2M0+I)
            DSPOLE = DSPOLE + VN0(I)
 36      CONTINUE
         TEM    = 4.0 / (DY*IM)
         DNPOLE = DNPOLE * TEM
         DSPOLE = DSPOLE * TEM
         
         DO 40 I=1,IMP3
            DIV(I+IJM1M0) = -DNPOLE
            DIV(I)        =  DSPOLE
 40      CONTINUE

         CALL WRAP(1,IJM0M0,DIV(1))
C
C    Specify the initial minus step model variables
C     
      DO J1X = 1, IJM0M0
         UN1(J1X) = UN0(J1X)
         VN1(J1X) = VN0(J1X)
         GTPO(J1X) = GTP(J1X)
         DIVO(J1X) = DIV(J1X)
      END DO
         FIRST = .FALSE.
         
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          START THE MAIN LOOP HERE                             C
C                                                               C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                      STEP 2                                   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 205 J=1,JMM1
         II = (J-1)*IMP3
         DO 205 I=1,IMP3
            WK7(II+I) = FCOR(J,2)
            WK8(II+I) = FCOR(J,1)
 205     CONTINUE

C
C ****************  AT U-GRID POINTS  *********************
C
C     CALCULATE V AT U-PTS  (SAVED IN WK3)
C     
      CALL UPTV(IJM0M0, VN0(1), WK3, TX1, TX2, TX3, IJP3M0)
C
C     CALCULATE YU (IN YUUPT) AT U-PTS
C
      DO 1226 I=1,IJM0M0
         WK5(I) = GTP(I) - DIV(I)*CDB
 1226 CONTINUE

      DO 217 I=I1P1,IJM1M1
         YUUPT(I+I2) = UN0(I) + (WK7(I)*WK3(I)
     *        - (WK5(I+1)-WK5(I))*CDU(I))*DL2DT
  217 CONTINUE

      CALL WRAP(I1P1,IJM1M0,YUUPT(I2+1))
      CALL FILPU(YUUPT(I2+1))
      CALL FILLU(YUUPT(1))
C
C ****************  AT V-GRID POINTS  *********************
C
C     CALCULATE U AT V-PTS (SAVED IN WK3)
C
      CALL VPTU(IJM0M0, UN0(1), WK3, TX1, TX2, TX3, IJP3M0)
C
C     CALCULATE YV (IN YVVPT) AT V-PTS
C
      DO 216 I=1,IJM1M0
         YVVPT(I+I2) = VN0(I) - (WK8(I) * WK3(I)
     *        + (WK5(I1+I) - WK5(I)) * (1.0/DY) ) * DL2DT
 216  CONTINUE

      CALL WRAP(1,IJM1M0,YVVPT(I2+1))
      CALL FILPV(YVVPT(I2+1))
      CALL FILLV(YVVPT(1))
C
C ****************  AT PHI-GRID POINTS  *********************
C
C  CALCULATE N(t=n), N(t=n-1) AND B1 AT PHI-GRID POINTS
C
      DO I=1, IJM0M0
         PPN(I) = DT*GTP(I)*DIV(I)
         PPNO(I) = DT*GTPO(I)*DIVO(I)
         BL(I+I2) = GTP(I)+FI0 - FI0*DL2DT*DIV(I)
      ENDDO

      CALL WRAP(1,IJM0M0,PPN(1))
      CALL WRAP(1,IJM0M0,PPNO(1))
      CALL WRAP(1,IJM0M0,BL(I2+1))
      CALL FILLP(BL(1))

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                            STEP 3                              C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          TIME EXTRAPOLATION TO OBTAIN MODEL VARIABLES 
C                       AT HALF TIME LEVELS
C
      DO 310 I=1,JSU1M0
         UST(I+I2)  = UN0(I) + TEXFa * (UN0(I) - UN1(I))
         VST(I+I2)  = VN0(I) + TEXFa * (VN0(I) - VN1(I))

         PPNH1(I+I2)     = PPN(I) + TEXFa * (PPN(I) - PPNO(I))
 310  CONTINUE

      DO 311 II=1,JUP0M0
         I = II + JSU1M0
         UST(I+I2)  = UN0(I) + TEXF1 * (UN0(I) - UN1(I))
         VST(I+I2)  = VN0(I) + TEXF1 * (VN0(I) - VN1(I))

         PPNH1(I+2)     = PPN(I) + TEXF1 * (PPN(I) - PPNO(I))
 311  CONTINUE

      DO 312 II=1,JNU1M0
         I = II + JNU*IMP3
         UST(I+I2)  = UN0(I) + TEXFa * (UN0(I) - UN1(I))
         VST(I+I2)  = VN0(I) + TEXFa * (VN0(I) - VN1(I))

         PPNH1(I+2)     = PPN(I) + TEXFa * (PPN(I) - PPNO(I))
 312  CONTINUE

      CALL FILPU(UST(I2+1))
      CALL FILPV(VST(I2+1))
      CALL FILLU(UST(1))
      CALL FILLV(VST(1))
      CALL FILLP(PPNH1(1))
C     
C     Updating the model variables at minus time levels (n-1)
C
      DO J1X = 1, IJM0M0
         UN1(J1X) = UN0(J1X)
         VN1(J1X) = VN0(J1X)
         DIVO(J1X) = DIV(J1X)
         PPNO(J1X) = PPN(J1X)
         GTPO(J1X) = GTP(J1X)
      END DO
C
C    NOW 'CENTER IN SPACE'.
C

      IF (ROTATE) THEN
         DO 320 I=1,IJM0M0
            VTS(I) = VST(I+I2)*DTOA
 320     CONTINUE
         DO 321 I=1,JSU1M0
            UTS(I) = UST(I+I2)*DTOA
 321     CONTINUE
         DO 322 I=1,JNU1M0
            UTS(I+JNU*IMP3) = UST(I+JNU*IMP3+I2)*DTOA
 322     CONTINUE
         DO 323 I=1,JUP0M0
            UTS(JSU1M0+I) = UST(JSU1M0+I+I2)*RH1U(JSU1M0+I)*DT
 323     CONTINUE
      ELSE
         DO 324 I=1,IJM0M0
            VTS(I) = VST(I+I2)*DTOA
            UTS(I) = UST(I+I2)*RH1U(I)*DT
 324     CONTINUE
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                       STEP 4                                    C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     
C    Estimating the mid-trajectory point locations for all the 
C                          U-grid points.
C
      CALL ZANGUR(UTS, VTS, CPHU, SPHU, PHU, XCR(1,2)
     *,              SPUH, SLUH, CBV, SBV
     *,              WK1, WK2, WK3, TX1, TX2, TX3)
C     
C     BILINEAR INTERPOLATION TO OBTAIN U(N+1/2) AND V(N+1/2) AT THE
C                    MID-TRAJECTORY LOCATIONS.
      do j=1,jmp3
      do i=1,imp3
      ust1(i,j)=ust((j-1)*imp3+i)
      vst1(i,j)=vst((j-1)*imp3+i)
      enddo
      enddo
      do j=1,jm
      do i=1,imp3
      spuh1(i,j)=spuh((j-1)*imp3+i)
      sluh1(i,j)=sluh((j-1)*imp3+i)
      c1a1(i,j)=c1a((j-1)*imp3+i)
      c2a1(i,j)=c2a((j-1)*imp3+i)
      enddo
      enddo
      print *, 'The First BILL***************************************'
      CALL BILL(UST1  , C1A1  ,  0.5, 1.0, 
     .          VST1  , C2A1  ,  1.0, 0.5,
     .          SPUH1 , SLUH1 ,    2, JMM1 )
      do j=1,jm
      do i=1,imp3
      c1a((j-1)*imp3+i)=c1a1(i,j)
      c2a((j-1)*imp3+i)=c2a1(i,j)
c      write(UNIT=15, FMT=19) c1a1(i,j), c2a1(i,j)
c      spuh((j-1)*imp3+i)=spuh1(i,j)
c      sluh((j-1)*imp3+i)=sluh1(i,j)
      enddo
      enddo
c      do j=1,jmp3
c      do i=1,imp3
c      ust((j-1)*imp3+i)=ust1(i,j)
c      vst((j-1)*imp3+i)=vst1(i,j)
c      enddo
c      enddo

c      CALL BILL(UST  , C1A  ,  0.5, 1.0, 
c     .          VST  , C2A  ,  1.0, 0.5,
c     .          SPUH , SLUH ,    2, JMM1 )
      IF (ROTATE) THEN
         DO 610 I=1,JSU2M0
            UMH(I1+I) = C1A(I1+I) * CBV(I1+I) - C2A(I1+I) * SBV(I1+I)
            VMH(I1+I) = C2A(I1+I) * CBV(I1+I) + C1A(I1+I) * SBV(I1+I)
 610     CONTINUE

         DO 612 II=1,JUP0M0
            I = II + JSU1M0
            UMH(I) = C1A(I) / COS(SPUH(I))
            VMH(I) = C2A(I)
 612     CONTINUE

         DO 1614 I=1,JNU2M0
            UMH(I+JNU*IMP3) = C1A(I+JNU*IMP3) * CBV(I+JNU*IMP3)
     *           - C2A(I+JNU*IMP3) * SBV(I+JNU*IMP3)
            VMH(I+JNU*IMP3) = C2A(I+JNU*IMP3) * CBV(I+JNU*IMP3)
     *           + C1A(I+JNU*IMP3) * SBV(I+JNU*IMP3)
 1614    CONTINUE

         DO 614 I=1,IJM2M0
            UMH(I1+I) = UMH(I1+I) * DTOA
            VMH(I1+I) = VMH(I1+I) * DTOA
 614     CONTINUE
      ELSE
         DO 617 I=I1P1,IJM1M0
            UMH(I) = C1A(I) / COS(SPUH(I)) * DTOA
            VMH(I) = C2A(I) * DTOA
 617     CONTINUE
      ENDIF

      CALL WRAP(I1P1,IJM1M0,UMH)
      CALL WRAP(I1P1,IJM1M0,VMH)
      CALL FILPU(UMH)
      CALL FILPU(VMH)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                          STEP 5                                  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C              FIND THE DEPARTURE POINT LOCATIONS FOR ALL 
C                      THE U-GRID POINT LOCATIONS.
C     
      CALL ZANGU(UMH, VMH, ONE, CPHU, SPHU, PHU, XCR(1,2)
     *,          SPHI,  SLAM, WK1, WK2)

C
C     BICUBIC INTERPOLATION TO Y(lambda)* and Y(theta)*
C     FOR ALL THE DEPARTURE POINTS OF U-GRID POINTS.

      do 761 I=1,IJM0M0
         yuupts(I) = -999.0
         yvupts(I) = -999.0
 761  CONTINUE

      CALL BIC1nu(YUUPT, YUUPTS, 0.5, 1.0,
     *            YVVPT, YVUPTS ,1.0, 0.5,
     *                     SPHI, SLAM, 2, JMM1)
C
C CALCULATE AUUPT(I,J) AT U-POINTS USING THE FULLY CENTERED APPROACH
C
      DO 9602 J=2,JMM1
         II = (J-1)*IMP3
         DO 9602 I=1,IMP3
            WK1(II+I) = BLMIK(I,1) - SLUH(II+I)
            WK2(II+I) = BLMIK(I,1) - SLAM(II+I)
 9602    CONTINUE

      DO 765 I=I1P1,IJM1M0
         WK4(I) = COS(SPUH(I))
         WK5(I) = COS(WK1(I))
         WK6(I) = SIN(WK1(I))
         WK7(I) = COS(WK2(I))
         WK8(I) = SIN(WK2(I))
         WK2(I) = SIN(SPHI(I))

         WKA(I) = SPHU(I) * SIN(SPUH(I))
         WK0(I) = CPHU(I) * WK4(I)
         WK0(I) =  1.0    / (WKA(I) + WK0(I) * WK5(I))

         WK9(I) = CPHU(I) * WK8(I) * WK4(I) * WK6(I) * WK0(I) + WK7(I)

         WK2(I) =  WK2 (I) * WK8(I)
     *       + (SPHU(I) * COS(SPHI(I)) - CPHU(I) * WK2(I) * WK7(I))
     *                  * WK6(I) * WK4(I) * WK0(I)

         AUUPT(I) =  WK9(I) * YUUPTS(I)  +  WK2(I) * YVUPTS(I)
 765  CONTINUE

      CALL WRAP(I1P1,IJM1M0,AUUPT)
      CALL FILPU(AUUPT)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                               STEP 6                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     Estimating the mid-trajectory point locations for all 
C                       the V-grid points
C
      CALL ZANGVR(UTS, VTS, CPHV, SPHV, PHV, XCR(1,1)
     *,              SPVH, SLVH, CBV, SBV
     *,              WK1, WK2, WK3, TX1, TX2, TX3)

C     
C        Bilinear interpolation to obtain velocities at 
C                  the mid-trajectory points
C
      do j=1,jmp3
      do i=1,imp3
      ust1(i,j)=ust((j-1)*imp3+i)
      vst1(i,j)=vst((j-1)*imp3+i)
      enddo
      enddo
      do j=1,jm
      do i=1,imp3
      spvh1(i,j)=spvh((j-1)*imp3+i)
      slvh1(i,j)=slvh((j-1)*imp3+i)
      enddo
      enddo
      print *, 'Second BILL********************************'
      CALL BILL(UST1, C1A1,  0.5, 1.0,
     .          VST1, C2A1,  1.0, 0.5,
     .	               SPVH1, SLVH1, 1, JMM1)
      do j=1,jm
      do i=1,imp3
      c1a((j-1)*imp3+i)=c1a1(i,j)
      c2a((j-1)*imp3+i)=c2a1(i,j)
c      write(UNIT=16, FMT=19) c1a1(i,j), c2a1(i,j)
c      spvh((j-1)*imp3+i)=spvh1(i,j)
c      slvh((j-1)*imp3+i)=slvh1(i,j)
      enddo
      enddo
c      do j=1,jmp3
c      do i=1,imp3
c      ust((j-1)*imp3+i)=ust1(i,j)
c      vst((j-1)*imp3+i)=vst1(i,j)
c      enddo
c      enddo

c      CALL BILL(UST, C1A,  0.5, 1.0,
c     .          VST, C2A,  1.0, 0.5,
c     .	               SPVH, SLVH, 1, JMM1)
      IF (ROTATE) THEN
         DO 616 I=1,JSV1M0
            UMH(I) =  C1A(I) * CBV(I) - C2A(I) * SBV(I)
            VMH(I) =  C2A(I) * CBV(I) + C1A(I) * SBV(I)
 616     CONTINUE

         DO 1618 I=1,JVP0M0
            UMH(I+JSV1M0) = C1A(I+JSV1M0) / COS(SPVH(I+JSV1M0))
            VMH(I+JSV1M0) = C2A(I+JSV1M0)
 1618    CONTINUE

         DO 620 I=1,JNV2M0
            UMH(I+JNV*IMP3) = C1A(I+JNV*IMP3) * CBV(I+JNV*IMP3)
     *           - C2A(I+JNV*IMP3) * SBV(I+JNV*IMP3)
            VMH(I+JNV*IMP3) = C2A(I+JNV*IMP3) * CBV(I+JNV*IMP3)
     *           + C1A(I+JNV*IMP3) * SBV(I+JNV*IMP3)
 620     CONTINUE

         DO 622 I=1,IJM1M0
            UMH(I) =  UMH(I) * DTOA
            VMH(I) =  VMH(I) * DTOA
 622     CONTINUE
      ELSE
         DO 1619 I=1,IJM1M0
            UMH(I) = C1A(I) / COS(SPVH(I)) * DTOA
            VMH(I) = C2A(I) * DTOA
 1619    CONTINUE
      ENDIF

      CALL WRAP(1,IJM1M0,UMH)
      CALL WRAP(1,IJM1M0,VMH)
      CALL FILPV(UMH)
      CALL FILPV(VMH)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                        STEP 7                                C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C           Estimating the departure point locations for 
C                       all the V-grid points
C
      CALL ZANGV(UMH, VMH, ONE, CPHV, SPHV, PHV, XCR(1,1)
     *,          SPHI, SLAM, WK1, WK2)

C     
C         Bicubic interpolation to obtain (Yu)* and (Yv)*     
C
      CALL BIC1nu(YVVPT, YVVPTS, 1.0, 0.5,
     *            YUUPT, YUVPTS, 0.5, 1.0,
     *                 SPHI, SLAM,1, JMM1)
C
C CALCULATE AVVPT(I,J) AT V-POINTS USING THE FULLY CENTERED APPROACH
C
      DO 9601 J=1,JMM1
         II = (J-1)*IMP3
         DO 9601 I=1,IMP3
            WK1(II+I) = BLMIK(I,2) - SLVH(II+I)
            WK2(II+I) = BLMIK(I,2) - SLAM(II+I)
 9601    CONTINUE

      DO 785 I=1,IJM1M0
         WK3(I) = SIN(SPVH(I))
         WK4(I) = COS(SPVH(I))
         WK5(I) = COS(WK1(I))
         WK6(I) = COS(WK2(I))
         WK7(I) = SIN(WK2(I))
         WK8(I) = SIN(SPHI(I))

         WK2(I) = CPHV(I)*WK4(I)*WK5(I) +   SPHV(I)*WK3(I)

         WK8(I) =  (WK8(I)*WK3(I)*WK6(I) 
     #            + WK4(I)*COS(SPHI(I))*WK5(I))
         WK9(I) = - WK3(I) * WK7(I)

         AVVPT(I) =   (WK9(I) * YUVPTS(I) 
     #              + WK8(I) * YVVPTS(I))/ WK2(I)
 785  CONTINUE

      CALL WRAP(1,IJM1M0,AVVPT)
      CALL FILPV(AVVPT)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                          STEP 8                               C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  CALCULATE A TERM IN FORCING FOR HELMHOLTZ EQUATION
C

      CALL UPTV(IJM0M0, AVVPT, WK2, TX1, TX2, TX3, IJP3M0)

      DO 814 I=1,IMP3
         WK7(I) = FDTHF(1,1)
         WK0(I) = ARR1(1,1)
 814  CONTINUE
      DO 815 J=2,JMM1
         II = (J-1)*IMP3
         DO 815 I=1,IMP3
            WK7(II+I) = FDTHF(J,1) 
            WK8(II+I) = FDTHF(J,2)
            WK9(II+I) = ARR1(J,2)
            WK0(II+I) = ARR1(J,1)
 815     CONTINUE

      DO 817 I=I1P1,IJM1M0
         UN0(I) = (WK2(I) * WK8(I) + AUUPT(I)) * WK9(I)
 817  CONTINUE

      DO 820 I=I1P1,IJM1M1
         WK3(I+1) = UN0(I+1) - UN0(I)
  820 CONTINUE

      DO 821 I=I1P1,IJM1M0,IMP3
         WK3(I)      = WK3(IM+I)
  821 CONTINUE

      CALL VPTU(IJM0M0, AUUPT, WK1, TX1, TX2, TX3, IJP3M0)

      DO 823 I=1,IJM1M0
         VN0(I) = (AVVPT(I) - WK1(I)*WK7(I)) * WK0(I)
  823 CONTINUE

      DO 824 I=1,IJM1M0
         WK2(I) = VN0(I) * CPHV(I)
  824 CONTINUE

      DO 825 I=I1P1,IJM1M0
         UNSS(I) =   WK3(I) * CDU(I) 
     #            + (WK2(I) - WK2(I-I1)) * CDV(I)
  825 CONTINUE

C
C CALCULATE TERM1 AT THE POLES
C
      TEMA = 0.0
      TEMB = 0.0
      DO 830 I=1,IM
         TEMA = TEMA + VN0(I)
         TEMB = TEMB - VN0(IJM2M0+I)
  830 CONTINUE

      TEMN     = 4.0 / (IM*DY)
      TEMS     = TEMN
      TM1NP = TEMN * TEMB
      TM1SP = TEMS * TEMA

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                       STEP  9                                C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     
C     Estimating the mid-trajectory point locations for 
C             all the PHI-grid points     
C
      CALL ZANGPR(UTS, VTS, HALF, CPHU, SPHU, PHU, XCR(1,1), SPHI
     *,           SLAM, CBV, SBV, WK1, WK2, WK3, WK4, WK5, WK6, WK7)
C     
C     Bilinear interpolation to obtain velocities at the 
C                    mid-trajectory points
C
      do j=1,jmp3
      do i=1,imp3
      ust1(i,j)=ust((j-1)*imp3+i)
      vst1(i,j)=vst((j-1)*imp3+i)
      enddo
      enddo
      do j=1,jm
      do i=1,imp3
      sphi1(i,j)=sphi((j-1)*imp3+i)
      slam1(i,j)=slam((j-1)*imp3+i)
c      c1a1(i,j)=c1a((j-1)*imp3+i)
c      c2a1(i,j)=c2a((j-1)*imp3+i)
      enddo
      enddo
      print *, 'Third BILL**********************************'
      CALL BILL(UST1, C1A1,  0.5, 1.0,
     .          VST1, C2A1,  1.0, 0.5,
     .	               SPHI1, SLAM1, 2, JMM1)
      do j=1,jm
      do i=1,imp3
      c1a((j-1)*imp3+i)=c1a1(i,j)
      c2a((j-1)*imp3+i)=c2a1(i,j)
c      WRITE(UNIT=17, FMT=19) c1a1(i,j), c2a1(i,j)
 19   FORMAT(1X,f12.10,2X,f12.10)
c      SPHI((j-1)*imp3+i)=SPHI1(i,j)
c      SLAM((j-1)*imp3+i)=SLAM1(i,j)
      enddo
      enddo
c      do j=1,jmp3
c      do i=1,imp3
c      ust((j-1)*imp3+i)=ust1(i,j)
c      vst((j-1)*imp3+i)=vst1(i,j)
c      enddo
c      enddo


c      CALL BILL(UST,  C1A, 0.5, 1.0,  
c     .          VST,  C2A, 1.0, 0.5,
c     .               SPHI, SLAM,2, JMM1)
      IF (ROTATE) THEN
         DO 2610 I=1,JSU2M0
            UMH(I1+I) = C1A(I1+I) * CBV(I1+I) - C2A(I1+I) * SBV(I1+I)
            VMH(I1+I) = C2A(I1+I) * CBV(I1+I) + C1A(I1+I) * SBV(I1+I)
 2610    CONTINUE

         DO 2612 II=1,JUP0M0
            I = II + JSU1M0
            UMH(I) = C1A(I) / COS(SPHI(I))
            VMH(I) = C2A(I)
 2612    CONTINUE

         DO 2614 I=1,JNU2M0
            UMH(I+JNU*IMP3) = C1A(I+JNU*IMP3) * CBV(I+JNU*IMP3)
     *           - C2A(I+JNU*IMP3) * SBV(I+JNU*IMP3)
            VMH(I+JNU*IMP3) = C2A(I+JNU*IMP3) * CBV(I+JNU*IMP3)
     *           + C1A(I+JNU*IMP3) * SBV(I+JNU*IMP3)
 2614    CONTINUE

         DO 2619 I=1,IJM2M0
            UMH(I1+I) = UMH(I1+I) * DTOA
            VMH(I1+I) = VMH(I1+I) * DTOA
 2619    CONTINUE
      ELSE

         DO 2615 I=I1P1,IJM1M0
            UMH(I) = C1A(I) / COS(SPHI(I)) * DTOA
            VMH(I) = C2A(I) * DTOA
 2615    CONTINUE
      ENDIF

      CALL WRAP(1,IJM0M0,UMH)
      CALL WRAP(1,IJM0M0,VMH)

C
C   Bicubic interpolation to obtain Nm(n+1/2) at PHI-points
C     
      CALL BICPnu(VST(I2+1), PPNH1, PPNH, SPHI, SLAM,
     * DTOA, 1.0, 1.0, 2, JMM1)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                   STEP  10                                    C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     
C Estimating the departure point locations for the PHI-Grid points     
C
      CALL ZANGP(UMH, VMH, ONE, CPHU, SPHU, PHU, XCR(1,1)
     *,          SPHI, SLAM, WK1, WK2)
C
C   Bicubic interpolation to obtain B1*(n). 
C
      CALL BICPnu(VST(I2+1), BL,  BLS, SPHI, SLAM, 
     *       DTOA,1.0, 1.0, 2, JMM1)

      DO 712 I=1,IJM0M0
         BTS(I) = BLS(I) - PPNH(I)    
  712 CONTINUE

      CALL WRAP(1,IJM0M0,BTS(1))

      CALL GLMNP2(IMP3,JM,BTS,DX(1,2),DY,
     *                        ARE2,DX(1,1),BTSM)

      DO I= 1, IJM0M0
         BTS(I) = BTS(I) - BTSM
      ENDDO

      DO 751 I=1,IJM0M0
         WK1(I) = 0.0
         WK2(I) = 0.0
  751 CONTINUE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                             STEP 11                           C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C   Calculating the forcing term for the elliptic equation
C
      DO 752 I=I1P1,IJM1M0
         WF1(I)   = SQMU * ( (FI0-BTSM) + DL1DT*FI0*UNSS(I)
     *                         +(CDB*DL1DT*SQMU - 1.0)* BTS(I))  
 752  CONTINUE

      DO I=1, IMP3
         WF1(I) = SQMU*( (FI0-BTSM) + DL1DT*FI0*TM1SP
     *                  +(CDB*DL1DT*SQMU - 1.0)* BTS(I))
         WF1(I+IJM1M0) = SQMU*( (FI0-BTSM) + DL1DT*FI0*TM1NP
     *                  +(CDB*DL1DT*SQMU - 1.0)* BTS(I+IJM1M0))
      ENDDO              

      CALL WRAP(1,IJM0M0,WF1)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                           STEP 12                              C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C    solve the elliptic equation
C
      DO 1120 I=1,IJM0M0
         RHS(I,1) = WF1(I)
         RHS(I,2) = 0.0
 1120 CONTINUE

      IF (INIT1) THEN
         ITRG = 3*IMP2 + 1
         DO 1112 J=1,JM
            DSC1(J) = ARR1(J,2)
            DSC2(J) = 0.0
            DSC3(J) = ARR1(J,1)
            DSC4(J) = ARR1(J,2)*ARR1(J,2)*(1.0-FUSQ(J))*(DL1DT*TWOM)
     *                                   *CPH(J,2) / AE
            DSC5(J) = 0.0
 1112    CONTINUE
      ENDIF
      DO 1116 J=1,JM
         DSC6(J) = - SQMU
 1116 CONTINUE

      CALL ELLSOL(IMP3, IM, JM, AE, INIT1, INIT2, .TRUE.
     *,           DSC1, DSC2, DSC3, DSC4, DSC5, DSC6
     *,           RHS(1,1), PHI(1,1), WSV, WKTMP, IX)

      INIT1 = .FALSE.
      INIT2 = .FALSE.

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                          STEP 13                              C
C          Updating the model variables, time level: (n+1)      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 1125 I=1,IJM0M0
         GTP(I) = PHI(I,1) + CDB*DL1DT*SQMU*BTS(I)
 1125 CONTINUE

      CALL WRAP(1, IJM0M0, GTP(1))
      
      DO 1520 I=1,IJM0M0
         DIV(I) = (1.0/(DL1DT*FI0)) * (BTS(I)+BTSM-GTP(I)-FI0)
 1520 CONTINUE

      CALL WRAP(1, IJM0M0, DIV(1))
C
C     OBTAIN U(T+DT)
C
      DO 1535 I=1,IJM0M0
         WK0(I) =   GTP(I) -   DIV(I)*CDA
 1535 CONTINUE

      DO 1545 I=1,IJM0M1
         WK3(I) =   WK0(I+1) +   WK0(I)
 1545 CONTINUE

      DO 1550 I=I1P1,IJM1M1
         WK4(I) = WK3(I1+I) - WK3(I-I1)
         WK5(I) = WK0(I+1)  - WK0(I)
 1550 CONTINUE

      DO 1565 J=1,JM
         II = (J-1) * IMP3
         DO 1560 I=1,IMP3
            WK6(II+I) = ARR3(J)
            WK7(II+I) = ARR4(J)
 1560    CONTINUE
 1565 CONTINUE

      DO 1570 I=I1P1,IJM1M0
         UN0(I) =  UN0(I)
     *          -  WK5(I) * WK6(I) * RDLM
     *          -  WK4(I) * WK7(I) * (RDPH*0.25)
 1570 CONTINUE

      CALL WRAP(I1P1,IJM1M0,UN0(1))
      CALL FILPU(UN0(1))
C
C     OBTAIN V(T+DT)
C
      DO 1590 I=1,IJM0M2
         WK3(I+1) =   WK0(I+2) -   WK0(I)
 1590 CONTINUE

      DO 1580 I=1,IJM0M0,IMP3
         WK3(IMP2+I) =   WK3(I+2)
         WK3(I)      =   WK3(IM+I)
 1580 CONTINUE

      DO 1598 J=1,JM
         II = (J-1)*IMP3
         DO 1597 I=1,IMP3
            WK6(II+I) = ARR5(J)
            WK7(II+I) = ARR6(J)
 1597    CONTINUE
 1598 CONTINUE

      DO 1600 I=1,IJM1M0
         VN0(I) =  VN0(I)
     *        + (WK3(I)    + WK3(I1+I)) * WK6(I) * (RDLM*0.25)
     *        - (WK0(I1+I) - WK0(I))    * WK7(I) * RDPH
 1600 CONTINUE

      CALL WRAP(1,IJM1M0,VN0(1))
      CALL FILPV(VN0(1))

      RETURN
      END





