c %Z% %M% %I% (IDS) %G%

      subroutine step(dt, first)
 
      include "slgcm.slgprm"
      include "slgcm.slgcon"
      include "slgcm.slgprg"
      include "slgcm.slgclm"

c
      dimension
     *          wk1(imp3,jm), wk2(imp3,jm), wk3(imp3,jm)
     *,         wk4(imp3,jm), wk5(imp3,jm), wk6(imp3,jm)
     *,         wk7(imp3,jm), wk8(imp3,jm), wk9(imp3,jm)
     *,         wk0(imp3,jm), wka(imp3,jm)

      logical first
 
      dimension sluh(imp3,jm), slvh(imp3,jm)
     *, spuh(imp3,jm), spvh(imp3,jm)
     *,         slam(imp3,jm), sphi(imp3,jm)
     *,         cbv(imp3,jm),  sbv(imp3,jm)
 
      dimension    bls(imp3,jm),  bts(imp3,jm)
     *,         unss(imp3,jm)
 
     *,         ust(imp3,jmp3),   vst(imp3,jmp3)
     *,       yuupt(imp3,jmp3), yvvpt(imp3,jmp3)
     *,          bl(imp3,jmp3), ppnh1(imp3,jmp3)
 
      dimension tx1(imp3,jmp3),   tx2(imp3,jmp3),   tx3(imp3,jmp3)
 
      dimension uts(imp3,jm),    vts(imp3,jm)
     *,         umh(imp3,jm),    vmh(imp3,jm)
     *,       auupt(imp3,jm),  avvpt(imp3,jm)
     *,      yuupts(imp3,jm), yuvpts(imp3,jm)
     *,      yvupts(imp3,jm), yvvpts(imp3,jm)
 
      dimension ppn(imp3,jm), ppnh(imp3,jm), ppno(imp3,jm), wf1(imp3,jm)
 
 
      dimension  c1a(imp3,jm), c2a(imp3,jm)
 
 
      dimension  rhs(imp3,jm,2), phi(imp3,jm,2), wsv(100*im), ix(im)
      dimension  wktmp(16*imp3,jm)
 
      logical init1, init2
      data init1/.true./, init2/.true./
      data nstep/0/
      data ia/158/
      save

c 
c+     Debugging Purposes
c   
       integer count, dj, dk
c-    

c 
c+     Debugging Purposes
c   
c        count = 0
c        open (unit = 10, file = 'win.out', status = 'new' )
c        write (10, '(// 1X, T5)') 
c     &   'This is the O/P for Serial Program'
        print *, 'This is the O/P for Serial Program'
c- 
      nstep = nstep + 1
 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     step 1: ----continued                                         c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      call wrap(2,jmm1, un0)
      call wrap(1,   jmm1, vn0)
      call wrap(1,   jm, gtp)
 
      call filpu(un0(1,1))
      call filpv(vn0(1,1))
c     
c     calculate the initial divergence     
c
      if(first) then

         do j=2,jmm1
         do i=2,imp3
           wk1(i,j) = un0(i,j) - un0(i-1,j)
         enddo
            wk1(1,j) = wk1(imp1,j)
         enddo
c
         do i=1,imp3
         do j=1,jmm1
          wk5(i,j) = vn0(i,j) * cphv(i,j)
         enddo
         enddo
c
         do i=1,imp3
         do j=2,jmm1
            div(i,j) = wk1(i,j) * cdu(i,j) + 
     .                (wk5(i,j) - wk5(i,j-1)) * cdv(i,j)
         enddo
         enddo


 
c
c     calculate the divergence at the poles
c
         dnpole = 0.0
         dspole = 0.0
         do 36 i=1,im
            dnpole = dnpole + vn0(i,jmm1)
            dspole = dspole + vn0(i,1)
 36      continue
         tem    = 4.0 / (dy*im)
         dnpole = dnpole * tem
         dspole = dspole * tem
         
         do 40 i=1,imp3
            div(i,jm) = -dnpole
            div(i,1)    =  dspole
 40      continue
 
         call wrap(1,jm,div)
c
c    specify the initial minus step model variables
c     
c      do j1x = 1, ijm0m0
         do i=1,imp3
         do j=1,jm

         un1(i,j) = un0(i,j)
         vn1(i,j) = vn0(i,j)
         gtpo(i,j) = gtp(i,j)
         divo(i,j) = div(i,j)
      enddo
      enddo
         first = .false.        
      endif

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c          start the main loop here                             c
c                                                               c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                      step 2                                   c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
      do 205 j=1,jmm1
         do 205 i=1,imp3
            wk7(i,j) = fcor(j,2)
            wk8(i,j) = fcor(j,1)
 205     continue
 
c
c ****************  at u-grid points  *********************
c
c     calculate v at u-pts  (saved in wk3)
c     
      call uptv(imp3,jm, vn0(1,1), wk3, tx1, tx2, tx3, imp3,jmp3)
c
c     calculate yu (in yuupt) at u-pts
         do j=1,jm
         do i=1,imp3
            wk5(i,j) = gtp(i,j) - div(i,j)*cdb
         enddo
         enddo
c
         do j=2,jmm1
         do i=1,imp3
         yuupt(i,j+2) = un0(i,j) + (wk7(i,j)*wk3(i,j)
     *        - (wk5(i+1,j)-wk5(i,j))*cdu(i,j))*dl2dt
         enddo
         enddo
 
      call wrap(2,jmm1,yuupt(i2+1,1))
      call filpu(yuupt(i2+1,1))
      call fillu(yuupt(1,1))
c
c ****************  at v-grid points  *********************
c
c     calculate u at v-pts (saved in wk3)
c
      call vptu(imp3,jm, un0(1,1), wk3, tx1, tx2, tx3, imp3,jmp3)
c
c     calculate yv (in yvvpt) at v-pts
c
         do j=1,jmm1
         do i=1,imp3
         yvvpt(i,j+2) = vn0(i,j) - (wk8(i,j) * wk3(i,j)
     *        + (wk5(i,j+1) - wk5(i,j)) * (1.0/dy) ) * dl2dt
         enddo
         enddo

 
      call wrap(1,jmm1,yvvpt(i2+1,1))
      call filpv(yvvpt(i2+1,1))
      call fillv(yvvpt(1,1))
c
c ****************  at phi-grid points  *********************
c
c  calculate n(t=n), n(t=n-1) and b1 at phi-grid points
c

c      do i=1, ijm0m0
c         ppn(i,1) = dt*gtp(i,1)*div(i,1)
c         ppno(i,1) = dt*gtpo(i,1)*divo(i,1)
c        bl(i+i2,1) = gtp(i,1)+fi0 - fi0*dl2dt*div(i,1)
c      enddo
c
c
         do j=1,jm
         do i=1,imp3
         ppn(i,j) = dt*gtp(i,j)*div(i,j)
         ppno(i,j) = dt*gtpo(i,j)*divo(i,j)
         bl(i,j+2) = gtp(i,j)+fi0 - fi0*dl2dt*div(i,j)
         enddo
         enddo

 
      call wrap(1,jm,ppn(1,1))
      call wrap(1,jm,ppno(1,1))
      call wrap(1,jm,bl(1,3))
      call fillp(bl(1,1))


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                            step 3                              c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c          time extrapolation to obtain model variables 
c                       at half time levels
c
       
       do  js  =   1  ,    jsu-1
       do  i   =   1  ,    imp3
       ust(i,js+2)  = un0(i,js) + texfa * (un0(i,js)-un1(i,js))
       vst(i,js+2)  = vn0(i,js) + texfa * (vn0(i,js)-vn1(i,js))
       ppnh1(i,js+2)= ppn(i,js) + texfa * (ppn(i,js)-ppno(i,js))
       enddo
       enddo


 

         do j=jsu,jnu
         do i=1,imp3
          ust(i,j+2)  = un0(i,j) + texf1 * (un0(i,j) - un1(i,j))
          vst(i,j+2)  = vn0(i,j) + texf1 * (vn0(i,j) - vn1(i,j))
          ppnh1(i+2,j)     = ppn(i,j) + texf1 * (ppn(i,j) - ppno(i,j))
          enddo
          enddo

         do j=jnu+1,jm
         do i=1,imp3
         ust(i,j+2)  = un0(i,j) + texfa * (un0(i,j) - un1(i,j))
         vst(i,j+2)  = vn0(i,j) + texfa * (vn0(i,j) - vn1(i,j)) 
         ppnh1(i+2,j)= ppn(i,j) + texfa * (ppn(i,j) - ppno(i,j))
         enddo
         enddo
      call filpu(ust(1,3))
      call filpv(vst(1,3))
      call fillu(ust(1,1))
      call fillv(vst(1,1))
      call fillp(ppnh1(1,1))
c     updating the model variables at minus time levels (n-1)
c
         do j=1,jm
         do i=1,imp3
         un1(i,j) = un0(i,j)
         vn1(i,j) = vn0(i,j)
         divo(i,j) = div(i,j)
         ppno(i,j) = ppn(i,j)
         gtpo(i,j) = gtp(i,j)
         enddo
         enddo
c
c    now 'center in space'.
c
 
      if (rotate) then
         do j=1,jm
         do i=1,imp3
            vts(i,j) = vst(i,j+2)*dtoa
         enddo
         enddo

         do j=1,jsu-1
         do i=1,imp3
            uts(i,j) = ust(i,j+2)*dtoa
          enddo
          enddo


         do j=jnu+1,jm
         do i=1,imp3
            uts(i,j) = ust(i,j+2)*dtoa
          enddo
          enddo

c         do 323 i=1,jup0m0
c            uts(jsu1m0+i,1) = ust(jsu1m0+i+i2,1)*rh1u(jsu1m0+i,1)*dt
c 323     continue
         do j=jsu,jnu
         do i=1,imp3
            uts(i,j) = ust(i,j+2)*rh1u(i,j)*dt
          enddo
          enddo


      else
         do j=1,jm
         do i=1,imp3
            vts(i,j) = vst(i,j+2)*dtoa
            uts(i,j) = ust(i,j+2)*rh1u(i,j)*dt
          enddo
          enddo

      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,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.
 
      call bill(ust  , c1a  ,  0.5, 1.0, 
     .          vst  , c2a  ,  1.0, 0.5,
     .          spuh , sluh ,    2, jmm1 )

c
c+     Debugging Purposes
c
      count = count + 1
c      write( 10, 99 )  'Count = ', count
 99   format( // 1X, T10, A10, I4)
c      write( 10, 100 ) 'WINA'
 100  format( // 1X, T10, A8 // )
c1 101  format( 1X, 8('+', 4F16.8 //), 3F16.8 //)
 101  format( 1X, 8('+', 4F16.8 //))
 102  format( 1X, 3F16.8 )

c      write (10, 101) ((ust(dk, dj), dj = 1, jmp3), dk = 1, imp3)
c      write (10, 101) ((vst(dk, dj), dj = 1, jmp3), dk = 1, imp3)

C1       print *, 'Count = ', count
C1       print *, '+WOUA'
C1       print 101, ((c1a(dk, dj), dj = 1, jm), dk = 1, imp3)
C1       print *, '+WOUA OVER'
C1       print *, '+WOUB'
C1       print 101, ((c2a(dk, dj), dj = 1, jm), dk = 1, imp3)
C1       print *, '+WOUB OVER'
      
c      do dk = 1, imp3
c         do dj = 1, jm, 4
c            write( 10, 101 ) ust(dk, dj), ust(dk, dj+1)
c     &        ust(dk, dj+2), ust(dk, dj+3)
c         end do
c        
c         write( 10, 102 ) ust(dk, jm+1), ust(dk, jm+2)
c     &        ust(dk, jm+3)
c      end do
c      close(10)
c-
      if (rotate) then
          do j=2,jsu
          do i=1,imp3
            umh(i,j) =    c1a(i,j) * cbv(i,j) 
     .                  - c2a(i,j) * sbv(i,j)
            vmh(i,j) =    c2a(i,j) * cbv(i,j) 
     .                  + c1a(i,j) * sbv(i,j)
          enddo
          enddo


 
c         do 612 ii=1,jup0m0
c            i = ii + jsu1m0
c            umh(i,1) = c1a(i,1) / cos(spuh(i,1))
c            vmh(i,1) = c2a(i,1)
c 612     continue
          do j=jsu,jnu
          do i=1,imp3
            umh(i,j) = c1a(i,j) / cos(spuh(i,j))
            vmh(i,j) = c2a(i,j)
          enddo
          enddo

 
          do j=jnu+1,jmm1
          do i=1,imp3
            umh(i,j) = c1a(i,j) * cbv(i,j)
     *           - c2a(i,j) * sbv(i,j)
            vmh(i,j) = c2a(i,j) * cbv(i,j)
     *           + c1a(i,j) * sbv(i,j)
          enddo
          enddo


 
          do j=2,jmm1
          do i=1,imp3
            umh(i,j) = umh(i,j) * dtoa
            vmh(i,j) = vmh(i,j) * dtoa
          enddo
          enddo

      else
          do j=2,jmm1
          do i=1,imp3
            umh(i,j) = c1a(i,j) / cos(spuh(i,j)) * dtoa
            vmh(i,j) = c2a(i,j) * dtoa
          enddo
          enddo
         
      endif
 
      call wrap(2,jmm1,umh)
      call wrap(2,jmm1,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,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.
 
c     do 761 i=1,ijm0m0
c         yuupts(i,1) = -999.0
c         yvupts(i,1) = -999.0
c 761  continue
       do j=1,jm
       do i=1,imp3
         yuupts(i,j) = -999.0
         yvupts(i,j) = -999.0
       enddo
       enddo

      
      print *, "Calling bic1nu from step"
      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  j=2,jmm1
         do  i=1,imp3
            wk1(i,j) = blmik(i,1) - sluh(i,j)
            wk2(i,j) = blmik(i,1) - slam(i,j)
         enddo
         enddo
 

          do j=2,jmm1
          do i=1,imp3
         wk4(i,j) = cos(spuh(i,j))
         wk5(i,j) = cos(wk1(i,j))
         wk6(i,j) = sin(wk1(i,j))
         wk7(i,j) = cos(wk2(i,j))
         wk8(i,j) = sin(wk2(i,j))
         wk2(i,j) = sin(sphi(i,j))
 
         wka(i,j) = sphu(i,j) * sin(spuh(i,j))
         wk0(i,j) = cphu(i,j) * wk4(i,j)
         wk0(i,j) =  1.0    / (wka(i,j) + wk0(i,j) * wk5(i,j))
 
         wk9(i,j) = cphu(i,j) * wk8(i,j) 
     .              * wk4(i,j) * wk6(i,j) * wk0(i,j) 
     .              + wk7(i,j)
 
         wk2(i,j) =  wk2 (i,j) * wk8(i,j)
     *            + (sphu(i,j ) * cos(sphi(i,j)) -
     *               cphu(i,j  ) * wk2(i,j) * wk7(i,j))
     *                  * wk6(i,j) * wk4(i,j) * wk0(i,j)
 
         auupt(i,j) =  wk9(i,j) * yuupts(i,j)  +  wk2(i,j) * yvupts(i,j)
         enddo
         enddo

 
      call wrap(2,jmm1,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,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
      call bill(ust, c1a,  0.5, 1.0,
     .          vst, c2a,  1.0, 0.5,
     .                 spvh, slvh, 1, jmm1)
      if (rotate) then
          do j=1,jsv-1
          do i=1,imp3
            umh(i,j) =  c1a(i,j) * cbv(i,j) - c2a(i,j) * sbv(i,j)
            vmh(i,j) =  c2a(i,j) * cbv(i,j) + c1a(i,j) * sbv(i,j)
          enddo
          enddo
 
         do j=jsv,jnv
         do i=1,imp3
            umh(i,j) = c1a(i,j) / cos(spvh(i,j))
            vmh(i,j) = c2a(i,j)
          enddo
          enddo
 
          do j=jnv+1,jmm1
          do i=1,imp3
            umh(i,j) = c1a(i,j) * cbv(i,j)
     *           - c2a(i,j) * sbv(i,j)
            vmh(i,j) = c2a(i,j) * cbv(i,j)
     *           + c1a(i,j) * sbv(i,j)
          enddo
          enddo

          do j=1,jmm1
          do i=1,imp3
            umh(i,j) =  umh(i,j) * dtoa
            vmh(i,j) =  vmh(i,j) * dtoa
          enddo
          enddo
      else
          do j=1,jmm1
          do i=1,imp3
            umh(i,j) = c1a(i,j) / cos(spvh(i,j)) * dtoa
            vmh(i,j) = c2a(i,j) * dtoa

          enddo
          enddo


      endif
 
      call wrap(1,jmm1,umh)
      call wrap(1,jmm1,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,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 j=1,jmm1
          do i=1,imp3
            wk1(i,j) = blmik(i,2) - slvh(i,j)
            wk2(i,j) = blmik(i,2) - slam(i,j)
          enddo
          enddo


         do j=1,jmm1
         do i=1,imp3
         wk3(i,j) = sin(spvh(i,j))
         wk4(i,j) = cos(spvh(i,j))
         wk5(i,j) = cos(wk1(i,j))
         wk6(i,j) = cos(wk2(i,j))
         wk7(i,j) = sin(wk2(i,j))
         wk8(i,j) = sin(sphi(i,j))
 
         wk2(i,j) = cphv(i,j)*wk4(i,j)*wk5(i,j) +   sphv(i,j)*wk3(i,j)
 
         wk8(i,j) =  (wk8(i,j)*wk3(i,j)*wk6(i,j) 
     #            + wk4(i,j)*cos(sphi(i,j))*wk5(i,j))
         wk9(i,j) = - wk3(i,j) * wk7(i,j)
 
         avvpt(i,j) =   (wk9(i,j) * yuvpts(i,j) 
     #              + wk8(i,j) * yvvpts(i,j))/ wk2(i,j)
         enddo
         enddo


 
      call wrap(1,jmm1,avvpt)
      call filpv(avvpt)
 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                          step 8                               c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  calculate a term in forcing for helmholtz equation
c
 
      call uptv(imp3,jm, avvpt, wk2, tx1, tx2, tx3, imp3,jmp3)
      do  i=1,imp3
         wk7(i,1) = fdthf(1,1)
         wk0(i,1) = arr1(1,1)
      enddo

         do js=2,jmm1
         do i=1,imp3
         wk7(i,js) = fdthf(js,1) 
            wk8(i,js) = fdthf(js,2)
            wk9(i,js) = arr1(js,2)
            un0(i,js)=(wk2(i,js) * wk8(i,js) + auupt(i,js))
     .               * wk9(i,js)
            wk0(i,js) = arr1(js,1)
       enddo
       enddo

      do j=2,jmm1
      do i=1,imp2
         wk3(i+1,j) = un0(i+1,j) - un0(i,j)
      enddo
      enddo

       do j=2,jmm1
         wk3(1,j)      = wk3(imp1,j)
       enddo 

       
      call vptu(imp3,jm, auupt, wk1, tx1, tx2, tx3, imp3,jmp3)
 
      do j=1,jmm1
      do i=1,imp3
      vn0(i,j) = (avvpt(i,j) -  wk1(i,j)*wk7(i,j)) * wk0(i,j)
      wk2(i,j) =    vn0(i,j) * cphv(i,j)
      enddo
      enddo

      do j=2,jmm1
      do i=1,imp3
         unss(i,j) =   wk3(i,j) * cdu(i,j) 
     #            + (wk2(i,j) - wk2(i,j-1)) * cdv(i,j) 
      enddo
      enddo

 
c
c calculate term1 at the poles
c
      tema = 0.0
      temb = 0.0
      do 830 i=1,im
         tema = tema + vn0(i,1)
         temb = temb - vn0(i,jmm1)
  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,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
      call bill(ust,  c1a, 0.5, 1.0,  
     .          vst,  c2a, 1.0, 0.5,
     .               sphi, slam,2, jmm1)
      if (rotate) then
          do j=2,jsu-1
          do i=1,imp3 
            umh(i,j) = c1a(i,j) * cbv(i,j) - 
     .                    c2a(i,j) * sbv(i,j)
            vmh(i,j) = c2a(i,j) * cbv(i,j) + 
     .                    c1a(i,j) * sbv(i,j)
          enddo
          enddo
 
          do j=jsu,jnu
          do i=1,imp3
            umh(i,j) = c1a(i,j) / cos(sphi(i,j))
            vmh(i,j) = c2a(i,j)
          enddo
          enddo

         do j=jnu+1,jmm1
         do i=1, imp3

            umh(i,j) = c1a(i,j) * cbv(i,j)
     *                        - c2a(i,j) * sbv(i,j)
            vmh(i,j) = c2a(i,j) * cbv(i,j)
     *                      + c1a(i,j) * sbv(i,j)
         enddo
         enddo
          
 
          do j=2,jmm1
          do i=1,imp3
            umh(i,j) = umh(i,j) * dtoa
            vmh(i,j) = vmh(i,j) * dtoa
          enddo
          enddo


      else

          do j=2,jmm1
          do i=1,imp3
            umh(i,j) = c1a(i,j) / cos(sphi(i,j)) * dtoa
            vmh(i,j) = c2a(i,j) * dtoa
          enddo
          enddo
         
      endif
 
      call wrap(1,jm,umh)
      call wrap(1,jm,vmh)
 
c
c   bicubic interpolation to obtain nm(n+1/2) at phi-points
c     
      call bicpnu(vst(i2+1,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,1)
     *,          sphi, slam, wk1, wk2)
c
c   bicubic interpolation to obtain b1*(n). 
c
      call bicpnu(vst(i2+1,1), bl,  bls, sphi, slam, 
     *       dtoa,1.0, 1.0, 2, jmm1)
 
        do j=1,jm
        do i=1,imp3
         bts(i,j) = bls(i,j) - ppnh(i,j)
        enddo
        enddo

 
      call wrap(1,jm,bts(1,1))
 
      call glmnp2(imp3,jm,bts,dx(1,2),dy,
     *                        are2,dx(1,1),btsm)
 
        do j=1,jm
        do i=1,imp3
         bts(i,j) = bts(i,j) - btsm 
         wk1(i,j) = 0.0
         wk2(i,j) = 0.0
        enddo
        enddo

      
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                             step 11                           c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   calculating the forcing term for the elliptic equation
c
       do j=2,jmm1
       do i=1,imp3
         wf1(i,j)   = sqmu * ( (fi0-btsm) + dl1dt*fi0*unss(i,j)
     *                         +(cdb*dl1dt*sqmu - 1.0)* bts(i,j)) 
       enddo
       enddo
      do i=1, imp3
         wf1(i,1) = sqmu*( (fi0-btsm) + dl1dt*fi0*tm1sp
     *                  +(cdb*dl1dt*sqmu - 1.0)* bts(i,1))

         wf1(i,jm) = sqmu*( (fi0-btsm) + dl1dt*fi0*tm1np
     *                  +(cdb*dl1dt*sqmu - 1.0)* bts(i,jm))
      enddo              
 
      call wrap(1,jm,wf1)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                           step 12                              c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c    solve the elliptic equation
c
       do j=1,jm
       do i=1,imp3
         rhs(i,j,1) = wf1(i,j)
         rhs(i,j,2) = 0.0

       enddo
       enddo

 
      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, phi(1,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 j=1,jm
       do i=1,imp3
        gtp(i,j) = phi(i,j,1) + cdb*dl1dt*sqmu*bts(i,j)
       enddo
       enddo


 
      call wrap(1, jm, gtp(1,1))
      
       do j=1,jm
       do i=1,imp3
         div(i,j) = (1.0/(dl1dt*fi0)) * (bts(i,j)+btsm-gtp(i,j)-fi0)
       enddo
       enddo

 
      call wrap(1, jm, div(1,1))
c
c     obtain u(t+dt)
c
       do j=1,jm
       do i=1,imp3
         wk0(i,j) =   gtp(i,j) -   div(i,j)*cda
       enddo
       enddo

 

       do j=1,jm
       do i=1,imp2
         wk3(i,j) =   wk0(i+1,j) +   wk0(i,j) 
       enddo
       enddo

 


       do j=2,jmm1
       do i=1,imp2
         wk4(i,j) = wk3(i,j+1) - wk3(i,j-1)
         wk5(i,j) = wk0(i+1,j)  - wk0(i,j)
       enddo
       enddo


 
c       do j=1,jm
c       do i=1,imp3
c           wk6(i,j) = arr3(j)
c           wk7(i,j) = arr4(j)
c       enddo
c       enddo


 
       do j=1,jmm1
       do i=1,imp3
         un0(i,j) =  un0(i,j)
     *          -  wk5(i,j) * arr3(j) * rdlm
     *          -  wk4(i,j) * arr4(j) * (rdph*0.25) 
       enddo
       enddo
     
 
      call wrap(2,jm,un0)
      call filpu(un0(1,1))
c
c     obtain v(t+dt)
c
       do j=1,jm
       do i=1,imp1
        wk3(i+1,j) =   wk0(i+2,j) -   wk0(i,j)
       enddo
       enddo

 
       do j=1,jm
         wk3(imp3,j) =   wk3(i+2,j)
         wk3(1,j)    =   wk3(imp1,j)
       enddo
 
       do j=1,jm
       do i=1,imp3
         vn0(i,j) =  vn0(i,j)
     *        + (wk3(i,j)    + wk3(i,j+1)) * arr5(j) * (rdlm*0.25)
     *        - (wk0(i,j+1) - wk0(i,j))    * arr6(j) * rdph

       enddo
       enddo
      call wrap(1,jmm1,vn0)
      call filpv(vn0)
 
      return
      end
