      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
            ii = (j-1)*imp3 + i
             fiwk(i,j) = gtp(ii)+fmean
         enddo    
        enddo
      write(krs0,rec=irec)fiwk
      irec =irec+1
c
        do i=1,im
         do j=1,jm
            ii = (j-1)*imp3 + i
             fiwk(i,j) = gtpo(ii)+fmean
         enddo    
        enddo
      write(krs0,rec=irec)fiwk
      irec =irec+1
c
        do i=1,im
         do j=1,jm
            ii = (j-1)*imp3 + i
             fiwk(i,j) = un0(ii)
         enddo    
        enddo
      write(krs0,rec=irec)fiwk

      irec =irec+1

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

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

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

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

        do i=1,im
         do j=1,jm
            ii = (j-1)*imp3 + i
             fiwk(i,j) = divo(ii)
         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
            IIJJ = (JJ-1)*IMP3 + II
            UN0(IIJJ) = AAAU(II,JJ)
            VN0(IIJJ) = AAAV(II,JJ)
            GTP(IIJJ) = 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
            IIJJ = (JJ-1)*IMP3 + II
            UN0(IIJJ) = AU(II,JJ)
            VN0(IIJJ) = AV(II,JJ)
            GTP(IIJJ) = AF(II,JJ)
            IF(FMAX.LT.AF(II,JJ)) FMAX = AF(II,JJ)
         ENDDO
      ENDDO
      endif

      CALL WRAP(1, IJM0M0, GTP(1))
      CALL GLMNP2(IMP3, JM, GTP, DX, DY, ARE2, DXV, FMEAN)

      DO J1X = 1, IJM0M0
         GTP(J1X) = GTP(J1X) - FMEAN
      END DO

      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))

      CLOSE(KRSI)

      PRINT*, ' FMEAN =',FMEAN, ' FI0=', FI0, ' FMAX=', FMAX
                                                      
 711    format(10x,A10,F6.1)
 911    FORMAT(1X,10E13.6)

        RETURN
        END
