      subroutine diagout(krs0,iiii,dataname)
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                       #     # ######   #####   #####
c                       #     # #     # #     # #     #
c                       #     # #     # #       #
c                       ####### ######  #       #
c                       #     # #       #       #
c                       #     # #       #     # #     #
c                       #     # #        #####   #####
c                       
c                          #######  #####   #####
c                          #       #     # #     #
c                          #       #       #
c                          #####    #####   #####
c                          #             #       #
c                          #       #     # #     #
c                          #######  #####   #####
c                       
c                          ######     #     #####
c                          #     #   # #   #     #
c                          #     #  #   #  #     #
c                          #     # #     # #     #
c                          #     # ####### #     #
c                          #     # #     # #     #
c                          ######  #     #  #####
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       jose zero nasa gsfc dao -- zero@talisman.gsfc.nasa.gov
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      include "slgcm.slgprm"
      include "slgcm.slgprg"
      logical first0
      integer ilev(1),lmax(50)
      character*3  month(12)
      character*8  qname(50)
      character*40 qdesc(50)
      character*80 xlab
      character*20  dataname

      dimension fiwk(im,jm)
      data xlab/'barotropic semilagrangian model run'/
      data month  /'jan','feb','mar','apr','may','jun',
     +             'jul','aug','sep','oct','nov','dec'/
      data qdesc(1)/'no way'/
      data first0/.true./
      data irec/1/
c
c
       if(first0)then
c
       first0=.false.
c
c open unit for data file
c
        print *,'opening unit ',krs0,' file ',dataname
       open (krs0,file=dataname,
     +      form='unformatted',access='direct',recl=im*jm*4)
c
c miscelaneous constants for grads format table
c
c
              nfreq=1
              ncnt=4
              nymd=990101
              nhms=000000
              iyear  = nymd/10000
              imon   = nymd/100 - iyear*100
              iday   = nymd - iyear*10000 - imon*100
              ihr    = nhms/10000
              imin   = nhms/100 - ihr*100
              qname(1)='geop'
              qname(2)='geopo'
              qname(3)='un0'
              qname(4)='vn0'
              qname(5)='un1'
              qname(6)='vn1'
              qname(7)='div'
              qname(8)='divo'
c
             qdesc(1)='geop'
             qdesc(2)='geopo'
             qdesc(3)='un0'
             qdesc(4)='vn0'
             qdesc(5)='un1'
             qdesc(6)='vn1'
             qdesc(7)='div'
             qdesc(8)='divo'

        undef = -999.9
        dxx    = 360./im
         dyy    = 180./(jm-1)
         km=1
         ilev(1)=1
        nfld=8
        lmax(1)=1
c
c
c write actual grads table (in fort.81)
c
        print *,'writing grads table'
         write(81,5001) dataname
         write(81,5002) xlab
        write(81,5003) undef
        write(81,5004) im,dxx
        write(81,5005) jm,dyy
        write(81,5006) km,(ilev(k),k=1,km)
        write(81,5007) ncnt,ihr,imin,iday,month(imon),iyear,nfreq
        write(81,5008) nfld
        do n=1,nfld
        lmax(n)=1
        write(81,5009) qname(n),lmax(n),nits,qdesc(n)
        enddo
        write(81,5010)   
      endif
c
c unfortunately all grids must be im by jm for grads
c
        do i=1,im
         do j=1,jm
             fiwk(i,j) = gtp(i,j)+fmean
         enddo    
        enddo
      write(krs0,rec=irec)fiwk
      irec =irec+1
c
        do i=1,im
         do j=1,jm
             fiwk(i,j) = gtpo(i,j)+fmean
         enddo    
        enddo
      write(krs0,rec=irec)fiwk
      irec =irec+1
c
        do i=1,im
         do j=1,jm
             fiwk(i,j) = un0(i,j)
         enddo    
        enddo
      write(krs0,rec=irec)fiwk

      irec =irec+1

c
        do i=1,im
         do j=1,jm
             fiwk(i,j) = vn0(i,j)
         enddo    
        enddo
      write(krs0,rec=irec)fiwk
      irec =irec+1

c
        do i=1,im
         do j=1,jm
             fiwk(i,j) = un1(i,j)
         enddo    
        enddo
      write(krs0,rec=irec)fiwk
      irec =irec+1

c
        do i=1,im
         do j=1,jm
             fiwk(i,j) = vn1(i,j)
         enddo    
        enddo
      write(krs0,rec=irec)fiwk
      irec =irec+1

c
        do i=1,im
         do j=1,jm
             fiwk(i,j) = div(i,j)
         enddo    
        enddo
      write(krs0,rec=irec)fiwk
      irec =irec+1

        do i=1,im
         do j=1,jm
             fiwk(i,j) = divo(i,j)
         enddo    
        enddo
      write(krs0,rec=irec)fiwk
      irec =irec+1
      return

      
 1000 format(1x,'resolution (',i3,',',i3,')')  
 2000 format(1x,'current date: ',i6,1x,i6,2x,
     .          'beginning date: ',i6,1x,i6,2x,'nmax: ',i4)

 5001 format('dset  ',a20)
 5002 format('title  ',a80)
 5003 format('undef  ',f6.1)
 5004 format('xdef  ',i3,' linear  -180.0 ',f6.3)
 5005 format('ydef  ',i3,' linear   -90.0 ',f6.3)
 5006 format('zdef  ',i3,' levels  ',50(i4,1x))
 5007 format('tdef  ',i3,' linear  ',i2.2,':',i2.2,'z',i2.2,a3,i2,
     .       2x,i2,'hr')
 5008 format('vars  ',i3)
 5009 format(a8,2x,i3,2x,i2,2x,a40)
 5010 format('endvars')


 711  format(10x,a10,f6.1)
 811  format(10x,'time step=',f6.1)
       end

      subroutine init(krsi)

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



      character*10 abc
      dimension aaau(im,jm), aaav(im,jm), aaaf(im,jm)
c
c
c in order to give the model the "same" initial conditions
c they will be read in 144x91 dimensions and interpolated
c (if necessary) to the current model resolution. j0
c 
c
      dimension au(144,91), av(144,91), af(144,91)
      common /mean/ fmean

      read(krsi) abc
      read(krsi) au
      read(krsi) av
      read(krsi) af
      if(im.ne.144.or.jm.ne.91)then
      print *,'interpolating to ',im,jm
      rjm=jm
      rim=im
      zp = 1.
      call sph2sph(au,144,91,zp,46.,aaau,im ,jm,zp,(rjm+1)/2.,zp)
      call sph2sph(av,144,91,zp,46.,aaav,im ,jm,zp,(rjm+1)/2.,zp)
      call sph2sph(af,144,91,zp,46.,aaaf,im ,jm,zp,(rjm+1)/2.,zp)
      fmax = aaaf(1,1)
      do ii=1, im
         do jj=1, jm
            un0(ii,jj) = aaau(ii,jj)
            vn0(ii,jj) = aaav(ii,jj)
            gtp(ii,jj) = aaaf(ii,jj)
            if(fmax.lt.aaaf(ii,jj)) fmax = aaaf(ii,jj)
         enddo
      enddo
      print *,'fmax',fmax

      else

      fmax = af(1,1)
      do ii=1, im
         do jj=1, jm
            un0(ii,jj) = au(ii,jj)
            vn0(ii,jj) = av(ii,jj)
            gtp(ii,jj) = af(ii,jj)
            if(fmax.lt.af(ii,jj)) fmax = af(ii,jj)
         enddo
      enddo
      endif

      call wrap(1, jm, gtp)
      call glmnp2(imp3, jm, gtp, dx, dy, are2, dxv, fmean)

c      do j1x = 1, ijm0m0
c         gtp(j1x,1) = gtp(j1x,1) - fmean
c      end do

       do j=1,jm
       do i=1,imp3
         gtp(i,j) = gtp(i,j) - fmean
      end do
      end do

      call wrap(2,jmm1, un0)
      call wrap(1,jmm1, vn0)
      call wrap(1,jm  , gtp)
      call filpu(un0)
      call filpv(vn0)

      close(krsi)

      print*, ' fmean =',fmean, ' fi0=', fi0, ' fmax=', fmax
                                                      
 711    format(10x,a10,f6.1)
 911    format(1x,10e13.6)

        return
        end
