c    @(#)genrl.f	1.2
c
c
c
      subroutine genrl(s,zp,klevel)
c
      include 'qg_par.inc'
      common/scrat/r(mn0),r1(mn0),r2(mn0),r3(mn0)
      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/edges/icorn(4),ncorn(4),mast(4),nast(4)
      common/bsolv/ab(mx0),bb(mx0),cb(mx0),workb(mnwork)
      common/dynam/alpha,beta,rf
      common/cnn/n,nm1,nm2,np1
      dimension s(mn0),zp(mn0)
      dimension aac(5),fx(mx0),fy(mx0),inout(mx0)
      do 5 iside=1,4
        ib=icorn(iside)
        ip=ncorn(iside)
        aac(iside)=r(ib)-r(ip)
    5 continue
      aac(5)=aac(1)
      do 100 iside=1,4
        i1=icorn(iside)
        i2=i1+mast(iside)
        i3=i2+mast(iside)
c
c --- determine inflow points
c
	ist=mast(iside)
	if(ist.eq.1.or.ist.eq.-1)mn2=mm2
	if(ist.eq.m.or.ist.eq.-m)mn2=nm2
        do 10 ii=1,mn2
          inout(ii)=0
          if(s(i3).gt.s(i1)) inout(ii)=1
          i1=i2
          i2=i3
          i3=i3+mast(iside)
   10   continue
        ib=icorn(iside)+mast(iside)
        if(inout(1).eq.1) go to 15
        i5=ib+nast(iside)
        i6=i5+mast(iside)
        fy(1)=r(ib)-4.*r(i5)-r(i6)-0.5*aac(iside)
        bb(1)=7.
        cb(1)=2.
        go to 20
   15   continue
        fy(1)=r1(ib)-zp(ib)
        bb(1)=1.
        cb(1)=0.
   20   continue
	if(ist.eq.1.or.ist.eq.-1)mn3=m-3
	if(ist.eq.m.or.ist.eq.-m)mn3=n-3
        do 40 iz=2,mn3
          ib=ib+mast(iside)
          if(inout(iz).eq.1) go to 35
          i5=ib+nast(iside)
          i4=i5-mast(iside)
          i6=i5+mast(iside)
          fy(iz)=r(ib)-4.*r(i5)-r(i4)-r(i6)
          bb(iz)=8.
          ab(iz)=2.
          cb(iz)=2.
          go to 40
   35     continue
          fy(iz)=r1(ib)-zp(ib)
          bb(iz)=1.
          ab(iz)=0.
          cb(iz)=0.
   40   continue
        ib=ib+mast(iside)
        if(inout(mn2).eq.1) go to 45
        i5=ib+nast(iside)
        i4=i5-mast(iside)
        fy(mn2)=r(ib)-4.*r(i5)-r(i4)-0.5*aac(iside+1)
        bb(mn2)=7.
        ab(mn2)=2.
        go to 50
   45   continue
        fy(mn2)=r1(ib)-zp(ib)
        bb(mn2)=1.
        ab(mn2)=0.
   50   continue
c
c --- call the ncar tri-diagonal matrix inverter
c
        call trdi(mn2,ab,bb,cb,fy,fx,0,workb)
c
        ib=icorn(iside)
        do 60 iz=1,mn2
          ib=ib+mast(iside)
          r(ib)=fx(iz)
   60   continue
  100 continue
      do 120 iside=1,4
        ib=icorn(iside)
        i1=ib+nast(iside)
        i2=ib+mast(iside)
        if(s(i2).gt.s(i1)) go to 110
        r(ib)=0.25*aac(iside)-0.5*(r(i1)+r(i2))
        go to 120
  110   continue
        r(ib)=r1(ib)-zp(ib)
  120 continue
      return
      end
