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 j=1,jmm1
         do 10 i=1,imp2
         f(i,j) = 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 j=1,nlats
      do 1003 i=1,imp3
      z(i,j)         = u(i,j+1)
      z(i,nlats+j)   = u(i,jm-nlats-1+j)
      z(i,2*nlats+j) = v(i,j+1)
      z(i,3*nlats+j) = v(i,jm-nlats-1+j)
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 j=1,totlats
      do 1110 i=1,imp2
      a(i,j) = a(i,j)*s(i,j)
 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 j=1,nlats
      do 1645 i=1,imp3
      u(i,j+1)        = z(i,j)
      u(i,jm-nlats-1+j) = z(i,nlats+j)
      v(i,j+1)        = z(i,2*nlats+j)
      v(i,jm-nlats-1+j) = z(i,3*nlats+j)
1645  continue

c
      return
      end
      subroutine consl

      include "slgcm.slgprm"
      include "slgcm.slgclm"
      include "slgcm.slgcon"
      real wk1(imp3,jm)
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
      do 5009 i=1,imp3
      phv(i,j) =  ph(j,1)
      phu(i,j) =  ph(j,2)
      sphv(i,j) = sph(j,1)
      sphu(i,j) = sph(j,2)
      cphv(i,j) = cph(j,1)
      cphu(i,j) = cph(j,2)
      rh1v(i,j) = rh1jk(j,1)
      rh1u(i,j) = rh1jk(j,2)
      wk1(i,j) =  dx(j,2)
 5009 continue
 5010 continue
c
       do j=1,jmm2
       do i=1,imp3
       tem        = 1.0 / ((dy*dy)  * cphu(i,j+1))
       cd1(i,j+1) = 1.0 / (wk1(i,j+1)*wk1(i,j+1))
       cd2(i,j+1) = tem * cphv(i,j+1)
       cd3(i,j+1) = tem * cphv(i,j)
       enddo
       enddo


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
      do 2007 i=1,imp3
      cdu(i,j) = cdivu(j)
      cdv(i,j) = cdivv(j)
 2007 continue
c
      do 8888 j=1,jm
      do 8888 i=1,imp3
      xcr(i,j,1) = (i - 1.0) * dlm
      xcr(i,j,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,jm,fsl)
c
      return
      end
      subroutine wrap(ib,ie,xx)
      include "slgcm.slgprm"
      dimension xx(imp3,jm)
       do j=ib,ie
       xx(imp1,j)   = xx(1,j)
       xx(imp2,j)   = xx(2,j)
       xx(imp3,j)   = xx(3,j)
       enddo

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,1) = tem1 + tem2
      cde(ii+i,1) = tem1 - tem2
      cdn(ii+i,1) = tem3
      cds(ii+i,1) = 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
            tkeu = tkeu + 0.5*(un0(ii,jj)**2)*au(jj)*fi0
            pol = pol + 0.5*(gtp(ii,jj)**2)*ag(jj)/9.8
            gmass = gmass + gtp(ii,jj)*ag(jj)
         enddo
      enddo

      tkeu = tkeu + 0.5*(un0(1,1)**2)*au(1)*fi0
     #     + 0.5*(un0(imp3,jm)**2)*au(jm)*fi0
      pol = pol + 0.5*(gtp(1,1)**2)*ag(1)/9.8
     #          + 0.5*(gtp(imp3,jm)**2)*ag(jm)/9.8
      gmass = gmass + gtp(1,1)*ag(1) + gtp(imp3,jm)*ag(jm)

      gmass = gmass/(aaa0*9.81)

      do jj = 1, jmm1
         do ii = 1, im

            tkev = tkev + 0.5*(vn0(ii,jj)**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
      end
