c
c     @(#)meshob.f	1.1
c
      SUBROUTINE MESHOB(IWAL)
c
C----------------------------------------------------------------------
C  IWAL=1 - READ REF FIELD, IWAL=0 - READ BNDRY FIELD
C  THIS SUBROUTINE COMPUTES STREAMFUNCTION FIELD USING SLA DATA 
C  AND FIRST MODE EOF; 
C  READ VORTICITY FROM UNIT=52
C  COMPUTES TOPDENSITY FIELD BY EXTENTION IT FROM LEVEL 1;  
C  "AL" IS SLA ARRAY, "PROF" IS 1-ST EOF MODE VERTICAL PROFILE 
C  0.23162 - CONVERTION COEEFICIENT FROM SLA(MM) TO STREAMFUNCTION(SM2/SEC)
c     INCLUDE 'qg_par.inc'
      INCLUDE 'mladv.inc'
C
      dimension vor(108),al(mn0),
     +          alv(mn0), ald(mn0),prof(k0),hz1(k0)
      COMMON/PARM/KZ,M,MM1,MM2,MP1,MSQ,TIM,TMAX,TSTART,HY,REA,DT,
     +   XBASIN,YY(MX0),ICALC,KZP1,KZM1,IFTOP,IFBOT,IFTVV,IFREL,THETA,
     +   ifsbl,ifeva
      common/assfield/ch3(mn0),ch4(mn0,k0),ch5(mn0,k0),bnd3(mn0),
     +     bnd4(mn0,k0),bnd5(mn0,k0)
      common/reffield/ref3(mn0),ref4(mn0,k0),ref5(mn0,k0)
      common/pref/irefs
      common/asspar/isd,ibeg,iasday,iastep,ifsdo,imf,icit
      common/philtr/iltord,iltfrq,iltcnt
      common/cnn/n,nm1,nm2,np1
      common/threed/hz(k0),scsig,sigz(k01),eigval(k0),ainv(k0,k0),
     +      amat(k0,k0)
C
       DATA inach/0/
        if(inach.ne.0)go to 3
         inach=1
           do 1 i=2,kz
             hz1(i)=0.5*(hz(i-1)+hz(i))
    1   continue
            read(50,*)(prof(i),i=1,6)

    3   CONTINUE
            read(24,*)(al(j),j=1,msq)
c
c --- read alongbnd vort and compute streamfunction at 6 levels
c
             do 140 k=1,kz
             
           read(52,*)(vor(j),j=1,108)
c ----
         if(iwal.eq.1)then
           do 25 j=1,msq
             ref4(j,k)=0.23162*al(j)*prof(k)/prof(1)
   25   continue
           goto 35
        endif
c ----
           do 30 j=1,msq             
             bnd4(j,k)=0.23162*al(j)*prof(k)/prof(1)
   30   continue
   35 CONTINUE
c
c --- compute topdensity
c
        if(k.ne.1) go to 37

        do 36 j=1,msq
        ald(j)=sigz(1)*(al(j)*(1-prof(2)/prof(1)))/hz1(2)       
   36   continue
   37   continue

         if(iwal.eq.1)then

           do 125 j=1,msq
             ref5(j,k)=ald(j)
  125   continue
           goto 140
        endif
c ---    compute alongboundary vorticity
c
c ---   southern          
           do 126 j=1,m
             bnd5(j,k)=vor(j)*2.3162 
  126   continue
c ---  eastern
           ip=m
           do 127 j=2*m, msq-m, m
            ip=ip+1
            bnd5(j,k)=vor(ip)*2.3162 
  127   continue
c --- northern
           ip=2*m+n-1
           do 128 j=msq-m+1,msq 
           ip=ip-1
           bnd5(j,k)=vor(ip)*2.3162 
  128   continue
c --- western
           ip=2*(m+n)-3
           do 129 j=m+1, msq-2*m+1, m 
           ip=ip-1
           bnd5(j,k)=vor(ip)*2.3162 
  129   continue  
c
  130   continue
  140  CONTINUE
c --- compute topdensity
         if(iwal.eq.1)then
           do 225 j=1,msq
             ref3(j)=0.23162*scsig*ald(j)
  225   continue
           goto 240
        endif
c
c --- compute topdensity
           do 230 j=1,msq
              bnd3(j)=0.23162*scsig*ald(j)             
  230   continue
  240   CONTINUE
        return
        end 

