ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     this corresponds to the file slgcm.zangs                c
c         in the new version of the 3d model                  c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine uptv(lenx,leny, vin, vou, tx1, tx2, tx3, lenpx,lenpy)
c
      include "slgcm.slgprm"
c
      dimension vin(lenx,leny), vou(lenx,leny) 
     *,         tx1(lenpx,lenpy), tx2(lenpx,lenpy), tx3(lenpx,lenpy)
c
c     compute v at u-points i.e. vupt(i,j)
c
       do j=1,jmm1
       do i=1,imp3
       tx3(i,j+2)=vin(i,j)
       enddo
       enddo
c
      call fillv(tx3)
c
       do j=1,jmp3
       tx3(im  +1,j) = tx3(1,j)
       tx3(imp2,j) = tx3(2,j)
       tx3(imp3,j) = tx3(3,j)
       enddo

       
c
       do j=1,jmp3
       do i=1,imp2
       tx1(i,j) = tx3(i+1,j) + tx3(i,j)
       enddo
       enddo


       do j=1,jmp3
       tx1(imp3,j)=tx1(3,j)
       enddo


      if (fourpt) then

          do j=2,jmp1
          vou(i,j)   = (tx1(i ,j+1)+tx1(i,j+2))* 0.25
          enddo
      else
          do j=1,jmp3
          do i=1,im
          tx2(i+1,j) = tx3(i+3,j) + tx3(i,j)
          enddo
          enddo

          do j=1,jmp3
          tx2(1,j)      = tx2(imp1,j)
          tx2(imp2,j)   = tx2(3,j)
          tx2(imp3,j)   = tx2(4,j)
          enddo

    
         do j=1,jmm1
         do i=1,imp3
         vou(i,j) = (tx1(i,j+2) + tx1(i,j+1)) * c81
     *            + (tx2(i,j)    + tx2(i,j+3)) * c01
     *            - (tx2(i,j+2) + tx2(i,j+1) 
     *            + tx1(i,j+3) + tx1(i,j)) * c09
         enddo
         enddo


      endif
c
      return
      end

      subroutine vptu(lenx,leny,uin,uou,tx1,tx2,tx3,lenpx,lenpy)
c
      include "slgcm.slgprm"
c
      dimension uin(lenx,leny), uou(lenx,leny) 
     *,         tx1(lenpx,lenpy), tx2(lenpx,lenpy), tx3(lenpx,lenpy)
c
c      compute u at v-points i.e. uvpt(i,j)
c
c
       do j=1,jm
       do i=1,imp3
       tx3(i,j+2)=uin(i,j)
       enddo
       enddo

      call fillu(tx3)
c
       do j=1,jmp3
       tx3(imp1,j) = tx3(1,j)
       tx3(imp2,j) = tx3(2,j)
       tx3(imp3,j) = tx3(3,j)
       enddo

       do j=2,jmp3
       do i=1,imp2
       tx1(i+1,j) = tx3(i+1,j) + tx3(i,j)
       enddo
       enddo

       do j=2,jmp3
       tx1(1,j)=tx1(imp1,j)
       enddo


      if (fourpt) then
          do j=1,jmm1
          do i=1,imp3
          uou(i ,j )=(tx1(i,j+2) + tx1(i,j+3)) * 0.25
          enddo
          enddo

      else
          do j=2,jmp3
          do i=1,im
        tx2(i+2,j)  = tx3(i+3,j) + tx3(i,j)
          enddo
          enddo
      

           do j=2,jmp3
         tx2(1,j)      = tx2(imp1,j)
         tx2(2,j)    = tx2(imp2,j)
         tx2(imp3,j) = tx2(3,j)
          enddo


 
          do j=1,jmm1
          do i=1,imp3
         uou(i,j)  = (tx1(i,j+3)  + tx1(i,j+2)) * c81
     *           + (tx2(i,j+1)  + tx2(i,j+4)) * c01
     *           - (tx2(i,j+3)  + tx2(i,j+2)
     *           + tx1(i,j+4) + tx1(i,j+1))*c09
          enddo
          enddo



      endif
c
      return
      end
      subroutine anguvr(lenx,leny, cph, sph, 
     *    uin, vin, pho, cbv, sbv, cpa)
c
      dimension  pho(lenx,leny), cph(lenx,leny)
     *, sph(lenx,leny), cbv(lenx,leny), sbv(lenx,leny)
     *, uin(lenx,leny), vin(lenx,leny), cpa(lenx,leny)
c
c
c      do 10 i=1,lenx*leny
       do j=1,leny
       do i=1,lenx

      sinv   = sin(-vin(i,j))
      cosv   = cos(-vin(i,j))
c
      sinu   = cosv * sin(-uin(i,j))
      cosu   = cosv * cos(-uin(i,j))

      tema   = cph(i,j) * sinv + sph(i,j) * cosu
      temb   = cph(i,j) * cosu - sph(i,j) * sinv
c
      pho(i,j) = asin(tema)
c
      uin(i,j) = - atan2(sinu,temb)
c
      tem    = sqrt(1.0 - tema*tema) * cph(i,j)
     *       + sph(i,j) * tema * cos(-uin(i,j))
c
      cpa(i,j) =  1.0  /  cosv
      cbv(i,j) =  tem  *  cpa(i,j)
c
      sbv(i,j) = sph(i,j) * sin(-uin(i,j)) * cpa(i,j)
c
c   10 continue
      enddo
      enddo
c
      return
      end
      subroutine zangur(uin, vin, cphu, sphu, phu, xcr
     *,                 pho, blo, cbv, sbv
     *,                 cp2, wk1, wk2, tx1, tx2, tx3)
c
      include "slgcm.slgprm"
c
c
      dimension  uin(imp3,jm),  vin(imp3,jm)
     *,          pho(imp3,jm),  blo(imp3,jm)
     *,          cphu(imp3,jm), sphu(imp3,jm), phu(imp3,jm)
     *,          cbv(imp3,jm),  sbv(imp3,jm),  xcr(imp3,jm)
      dimension  wk1(imp3,jm),  wk2(imp3,jm),  cp2(imp3,jm)
     *,          tx1(imp3,jmp3),  tx2(imp3,jmp3),  tx3(imp3,jmp3)
c
c      compute v at u-points i.e. vupt(i,j)
c

       call uptv(imp3,jm, vin, wk1, tx1, tx2, tx3, imp3,jmp3)
c
       do j=1,jmm1
       do i=1,imp3
       wk1(i,j) =  wk1(i,j) * 0.5
       wk2(i,j) =  uin(i,j) * 0.5
       enddo
       enddo


      if (rotate) then
c
c     compute for south of thetam
c
c
         jc = i1 + 1
         call anguvr(imp3,jsu-2, cphu(1,2), sphu(1,2), wk2(1,2)
     *,  wk1(1,2), pho(1,2),  cbv(1,2),  sbv(1,2), cp2(1,2))
c
          do j=2,jsu-1
          do i=1,imp3
          cbv(i,j) = cbv(i,j) * cp2(i,j)
          sbv(i,j) = sbv(i,j) * cp2(i,j)
          enddo
          enddo


c
c        points equatorward of thetam
c
c
          do j=jsu,jnu
          do i=1,imp3
          pho(i,j) = phu(i,j) - wk1(i,j)
          enddo
          enddo


c
c        northern high latitudes (poleward of thetam)
c
         jn = jnu*imp3
         jc = jn + 1
c
         call anguvr(imp3,jmm1-jnu, cphu(1,jnu+1), 
     * sphu(1,jnu+1), wk2(1,jnu+1)
     *,  wk1(1,jnu+1), pho(1,jnu+1),  cbv(1,jnu+1)
     *,  sbv(1,jnu+1), cp2(1,jnu+1))
c
          do j=jnu+1,jm
          do i=1,imp3
          cbv(i,j) = cbv(i,j) * cp2(i,j)
          sbv(i,j) = sbv(i,j) * cp2(i,j)
          enddo
          enddo

      else

          do j=2,jmm1
          do i=1,imp3
         pho(i,j) = phu(i,j) - wk1(i,j)
          enddo
          enddo


      endif
c

c
          do j=2,jmm1
          do i=1,imp3
      blo(i,j) = xcr(i,j) - wk2(i,j)
          enddo
          enddo



c
      return
      end

      subroutine zangu(uin, vin, fac, cphu, sphu, phu, xcr
     *,                pho, blo, wk1, wk2)

      include "slgcm.slgprm"

      dimension  uin(imp3,jm),   vin(imp3,jm)
     *,          pho(imp3,jm),   blo(imp3,jm)
     *,          cphu(imp3,jm), sphu(imp3,jm)
     *, phu(imp3,jm), xcr(imp3,jm)
      dimension  wk1(imp3,jm),  wk2(imp3,jm)

          do j=1,jm
          do i=1,imp3
          wk1(i,j) = uin(i,j) * fac
          wk2(i,j) = vin(i,j) * fac
          enddo
          enddo



 

      if (rotate) then

c        compute for south of thetam

         call anguvp(imp3,jsu-2,  cphu(1,2), sphu(1,2)
     *,                       wk1(1,2),  wk2(1,2), pho(1,2))

         do j=jsu,jnu
         do i=1,imp3
         pho(i,j) = phu(i,j) - wk2(i,j)
         enddo
         enddo

c        northern high latitudes (poleward of thetam)

         jc=jnu * imp3 + 1
         call anguvp(imp3,jmm1-jnu,  cphu(1,jnu+1), sphu(1,jnu+1)
     *,                       wk1(1,jnu+1),  wk2(1,jnu+1)
     *, pho(1,jnu+1))
      else
          do j=2,jmm1
          do i=1,imp3
          pho(i,j) = phu(i,j) - wk2(i,j)
          enddo
          enddo
      endif


      do j=2,jmm1
      do i=1,imp3
      blo(i,j) = xcr(i,j) - wk1(i,j)
      enddo
      enddo



      return
      end

      subroutine zangvr(uin, vin, cphv, sphv, phv, xcr
     *,                 pho, blo, cbv, sbv
     *,                 cp1, wk1, wk2, tx1, tx2, tx3)

      include "slgcm.slgprm"
c
      dimension  uin(imp3,jm),  vin(imp3,jm)
     *,          pho(imp3,jm),  blo(imp3,jm)
     *,          cphv(imp3,jm), sphv(imp3,jm), phv(imp3,jm)
     *,          cbv(imp3,jm),  sbv(imp3,jm),  xcr(imp3,jm)
      dimension  wk1(imp3,jm),  wk2(imp3,jm),  cp1(imp3,jm)
     *,          tx1(imp3,jmp3),  tx2(imp3,jmp3),tx3(imp3,jmp3)
c
c      compute u at v-points i.e. uvpt(i,j)
c
      call vptu(imp3,jm, uin, wk1, tx1, tx2, tx3, imp3,jmp3)
       do j=1,jmm1
       do i=1,imp3
       wk1(i,j) =  wk1(i,j) * 0.5
       wk2(i,j) =  vin(i,j) * 0.5
       enddo
       enddo


c
      if (rotate) then
c
c     compute for south of thetam
c

         call anguvr(imp3,jsv-1, cphv, sphv,
     *   wk1, wk2, pho, cbv, sbv, cp1)
c
c     compute for the equatorial points
c

          do j=jsv,jnv
          do i=1,imp3
          pho(i,j) =  phv(i,j) - wk2(i,j)
          enddo
          enddo


c
c        compute for latitude north of thetam
c
c         jc = jnv*imp3 + 1

         call anguvr(imp3,jmm1-jnv, cphv(1,jnv+1)
     *, sphv(1,jnv+1), wk1(1,jnv+1)
     *, wk2(1,jnv+1), pho(1,jnv+1), cbv(1,jnv+1)
     *, sbv(1,jnv+1), cp1(1,jnv+1))

      else
          do j=1,jmm1
          do i=1,imp3
         pho(i,j) =  phv(i,j) - wk2(i,j)
          enddo
          enddo

          
      endif

       do j=1,jmm1
       do i=1,imp3
       blo(i,j) = -wk1(i,j) + xcr(i,j)
       enddo
       enddo

    

      return
      end
      subroutine zangv(uin, vin, fac, cphv, sphv, phv, xcr
     *,                pho, blo, wk1, wk2)

      include "slgcm.slgprm"

      dimension  uin(imp3,jm),  vin(imp3,jm)
     *,          pho(imp3,jm),  blo(imp3,jm)
     *,         cphv(imp3,jm), sphv(imp3,jm)
     *, phv(imp3,jm), xcr(imp3,jm)
      dimension  wk1(imp3,jm),  wk2(imp3,jm)

       do j=1,jmm1
       do i=1,imp3
       wk1(i,j) = uin(i,j) * fac
       wk2(i,j) = vin(i,j) * fac
       enddo
       enddo


   

      if (rotate) then
c
c        compute for south of thetam
c
         call anguvp(imp3,jsv-1, cphv, sphv, wk1, wk2, pho)
c
c        compute for the equatorial points
c
          do j=jsv,jnv
          do i=1,imp3
          pho(i,j) =  phv(i,j) - vin(i,j)
          enddo
          enddo


c
c        compute for latitude north of thetam
c
c        jc=jnv * imp3 + 1
         call anguvp(imp3,jmm1-jnv,  cphv(1,jnv+1), sphv(1,jnv+1)
     *,                       wk1(1,jnv+1),  wk2(1,jnv+1), pho(1,jnv+1))
      else
          do j=1,jmm1
          do i=1,imp3
         pho(i,j) =  phv(i,j) - vin(i,j)
          enddo
          enddo

 
      endif

          do j=1,jmm1
          do i=1,imp3
          blo(i,j) = -wk1(i,j) + xcr(i,j)
          enddo
          enddo
   

      return
      end
      subroutine zangpr(uin,  vin,  fac
     *,                cphu, sphu, phu, xcr
     *,                pho,  blo,  cbv, sbv
     *,                cp2, wk1, wk2, wk3, wk4, wk5, wk6)

      include "slgcm.slgprm"

      dimension  uin(imp3,jm),  vin(imp3,jm)
     *,          pho(imp3,jm),  blo(imp3,jm)
     *,         cphu(imp3,jm), sphu(imp3,jm)
     *, phu(imp3,jm), xcr(imp3,jm)
      dimension  wk1(imp3,jm),  wk2(imp3,jm)
     *, wk3(imp3,jm),  wk4(imp3,jm)
     *,          wk5(imp3,jm),  wk6(imp3,jm), cp2(imp3,jm)
     *,          cbv(imp3,jm),  sbv(imp3,jm)

c
c      first compute u & v at fi - points
c
c at j=2
c
       do j=1,jmm1
       do i=1,imp2
       wk1(i+1,j) = uin(i+1,j) + uin(i,j)
       enddo
       wk1(1,j)   = wk1(imp1,j)
       enddo
   

       

       do j=2,jm
       do i=1,imp3
       wk2(i,j) = vin(i,j-1) + vin(i,j)
       enddo
       enddo

       do i=1,imp3
       wk4(i,2) = wk1(i,2) * half
       wk3(i,2) = wk2(i,2) * half
       enddo

      

c
c at j=3,jmm2
c
       do j=1,jm-4
       do i=1,imp3
       wk5(i,j+2) = vin(i,j+3) + vin(i,j)
       enddo
       enddo



       do j=1,jm-4
       do i=1,im
       wk6(i+2,j+2) = uin(i+3,j+2) + uin(i,j+2)
       enddo
       enddo

       do j=3,jm-2
       wk6(1,j)      = wk6(imp1,j)
       wk6(2,j)      = wk6(imp1+1,j)
       wk6(imp2+1,j) = wk6(3,j)
       enddo

    


        do j=3,jmm2
        do i=1,imp3
        wk3(i,j) =  wk2(i,j) * c9o16 - wk5(i,j) * c1o16
        wk4(i,j) =  wk1(i,j) * c9o16 - wk6(i,j) * c1o16
        enddo
        enddo


        

c
c at j=jmm1
c
       do i=1,imp3
       wk4(i,jmm1) = wk1(i,jmm1) * half
       wk3(i,jmm1) = wk2(i,jmm1) * half
       enddo



        do j=2,jmm1
        do i=1,imp3
        wk4(i,j) = wk4(i,j) * fac
        wk3(i,j) = wk3(i,j) * fac
        enddo
        enddo




      if (rotate) then
c
c     compute for south of thetam
c
c         jc = i1 + 1
         call anguvr(imp3,jsu-2, cphu(1,2), sphu(1,2), wk4(1,2)
     *,              wk3(1,2), pho(1,2)
     *,  cbv(1,2),  sbv(1,2), cp2(1,2))

          do j=2,jsu-1
          do i=1,imp3
          cbv(i,j) = cbv(i,j) * cp2(i,j)
          sbv(i,j) = sbv(i,j) * cp2(i,j)
          enddo
          enddo




c
c        points equatorward of thetam
c
          do j=jsu,jnu
          do i=1,imp3
          pho(i,j) = phu(i,j) - wk3(i,j)
          enddo
          enddo
c
c        northern high latitudes (poleward of thetam)
c

         call anguvr(imp3,jmm1-jnu, cphu(1,jnu+1)
     *, sphu(1,jnu+1), wk4(1,jnu+1)
     *, wk3(1,jnu+1), pho(1,jnu+1),  cbv(1,jnu+1)
     *,  sbv(1,jnu+1), cp2(1,jnu+1))
c
          do j=jnu+1,jmm1
          do i=1,imp3
          cbv(i,j) = cbv(i,j) * cp2(i,j)
          sbv(i,j) = sbv(i,j) * cp2(i,j)
          enddo
          enddo


          
      else

          do j=2,jmm1
          do i=1,imp3
          pho(i,j) = phu(i,j) - wk3(i,j) 
          enddo
          enddo       
      endif

       do j=2,jmm1
       do i=1,imp3
       blo(i,j) = xcr(i,j) - wk4(i,j)  
       enddo
       enddo  

      return
      end
      subroutine zangp(uin,  vin,  fac
     *,                cphu, sphu, phu, xcr
     *,                pho,  blo,  wk1, wk2)

      include "slgcm.slgprm"
c
c
      dimension  uin(imp3,jm),  vin(imp3,jm)
     *,          pho(imp3,jm),  blo(imp3,jm)
     *,         cphu(imp3,jm), sphu(imp3,jm)
     *, phu(imp3,jm), xcr(imp3,jm)
      dimension  wk1(imp3,jm),  wk2(imp3,jm)
c
       do j=2,jmm1
       do i=1,imp3
       wk1(i,j) = uin(i,j) * fac
       wk2(i,j) = vin(i,j) * fac
       enddo
       enddo



c
      if (rotate) then
c
c        compute for south of thetam
c
c         jc = i1+1
         call anguvp(imp3,jsu-2, cphu(1,2), sphu(1,2)
     *,                      wk1(1,2),  wk2(1,2), pho(1,2))
c
c         points equatorward of thetam
          do j=jsu,jnu
          do i=1,imp3
          pho(i,j) = phu(i,j) - wk2(i,j)
          enddo
          enddo
c
c    northern high latitudes (poleward of thetam)
c
c         jc    = jnu*imp3 + 1
         call anguvp(imp3,jmm1-jnu, cphu(1,jnu+1), sphu(1,jnu+1)
     *,                      wk1(1,jnu+1),  wk2(1,jnu+1), pho(1,jnu+1))
      else
          do j=2,jmm1
          do i=1,imp3
          pho(i,j) = phu(i,j) - wk2(i,j)
          enddo
          enddo



      endif
       do j=2,jmm1
       do i=1,imp3
       blo(i,j) = xcr(i,j) - wk1(i,j)
       enddo
       enddo


      return
      end
      subroutine anguvp(lenx,leny, cph, sph, uin, vin, pho)
c
      dimension  pho(lenx,leny), cph(lenx,leny),
     *  sph(lenx,leny), uin(lenx,leny), vin(lenx,leny)
c
      do j=1,leny
      do i=1,lenx
      sinv   = sin(-vin(i,j))
      cosv   = cos(-vin(i,j))
c
      sinu   = cosv * sin(-uin(i,j))
      cosu   = cosv * cos(-uin(i,j))

      tema   = cph(i,j) * sinv + sph(i,j) * cosu
      temb   = cph(i,j) * cosu - sph(i,j) * sinv
c
      pho(i,j) = asin(tema)
c
      uin(i,j) = - atan2(sinu,temb)

      enddo
      enddo
c
      return
      end
