c
c     @(#)mesho_orig.f	1.1
c
      SUBROUTINE MESHO_ORIG(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  COMPUTES VORTICITY FIELD AS SUM OF RELATIVE AND THERMAL VORTICITY;
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 COEFFICIENT FROM SLA(MM) TO STREAMFUNCTION(CM^2/SEC)
c     INCLUDE 'qg_par.inc'
      INCLUDE 'mladv.inc'
C
      dimension 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
          if(iwal.eq.1)read(51,*)(al(j),j=1,msq)

          if(iwal.eq.0)read(24,*)(al(j),j=1,msq)
c
c --- calculate streamfunction 
c
             do 140 k=1,kz
c 
         if(iwal.eq.1)then
           do 25 j=1,msq
             ref4(j,k)=0.23162*al(j)*prof(k)/prof(1)
   25   continue
          call filter(ref4(1,k),m,n,iltord)
          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
c
c --- compute  vorticity 
c
c --- compute thermal vorticity at all grid points 
c
   37   continue
         do 39 j=1,msq

             if(k.eq.1)then
          alv(j)=scsig/hz(1)*(ald(j)-sigz(2)*(al(j)*(1-
     +    prof(2)/prof(1)))/hz1(2))
              go to 39
             endif

             if(k.eq.kz)then
          alv(j)=scsig/hz(kz)*(sigz(kz)*(al(j)/prof(1)*(prof(5)-
     +    prof(6)))/hz1(kz))
            go to 39
          endif

        alv(j)=scsig/hz(k)*(sigz(k)*(al(j)/prof(1)*
     +  (prof(k-1)-prof(k)))
     +  /hz1(k)-sigz(k+1)*(al(j)/prof(1)*(prof(k)-
     +   prof(k+1)))/hz1(k+1))
   39   continue
   40   continue
c --- compute total vorticity in internal grid points
              ip=m+1
         do 50 i=2,n-1
           do 45 j=2,m-1
              ip=ip+1
         alv(ip)=alv(ip)+prof(k)/prof(1)*(al(ip+1)+
     +           al(ip-1)-al(ip)*4+al(ip+m)+
     +           al(ip-m))/hy/hy
   45  continue
              ip=ip+2
   50   continue
c
c --- southern wall
c
           do 60 i=2,m-1
        alv(i)=alv(i)+prof(k)/prof(1)*(al(i+1)-
     +         al(i)*2+al(i-1))/hy/hy
   60   continue
c --- northern wall
              ip=msq-m+1
           do 70 j=2,m-1
              ip=ip+1
        alv(ip)=alv(ip)+prof(k)/prof(1)*(al(ip+1)-
     +          al(ip)*2+al(ip-1))/hy/hy
   70   continue
c --- western wall
             ip=1
           do 80 j=2,n-1
             ip=ip+m    
        alv(ip)=alv(ip)+prof(k)/prof(1)*(al(ip+m)-
     +          al(ip)*2+al(ip-m))/hy/hy
c --- eastern wall
   80  continue
             ip=m
            do 90 j=2,n-1
             ip=ip+m
        alv(ip)=alv(ip)+prof(k)/prof(1)*(al(ip+m)-
     +          al(ip)*2+al(ip-m))/hy/hy
   90   continue
c
c --- corner points
c
        CONTINUE
            ip=1
            alv(ip)=0.5*(alv(ip+1)+alv(ip+m))
            ip=m
            alv(ip)=0.5*(alv(ip-1)+alv(ip+m))
            ip=msq-m+1
            alv(ip)=0.5*(alv(ip-m)+alv(ip+1))
            ip=msq
            alv(ip)=0.5*(alv(ip-1)+alv(ip-m))
c ----
         if(iwal.eq.1)then
           do 125 j=1,msq
             ref5(j,k)=0.23162*alv(j)
  125   continue
          call filter(ref5(1,k),m,n,iltord)
           goto 140
        endif
c ----
           do 130 j=1,msq
             bnd5(j,k)=0.23162*alv(j)
  130   continue
c --            call filter(bnd5(1,k),m,n,iltord)
  140  CONTINUE
c ---
         if(iwal.eq.1)then
           do 225 j=1,msq
             ref3(j)=0.23162*scsig*ald(j)
  225   continue
          call filter(ref3,m,n,iltord)
           goto 240
        endif
c ----
           do 230 j=1,msq
              bnd3(j)=0.23162 *scsig*ald(j)
  230   continue
c             call filter(bnd3,m,n,iltord)
  240   CONTINUE
        return
        end 
