c
c     @(#)infld.f	1.1
c
      SUBROUTINE INFLD(R,K,T,ISORV,IFBND)
C----------------------------------------------------------------------
C
C  THIS ROUTINE DETERMINES ROSSBY WAVE PARAMETRS, AND STREAMFUNCTION
C  AND VORTICITY FOR INITIALIZATION, BOUNDARY CONDITIONS, AND 
C  PERTURBATION FIELDS.
C
C
C  ARGUMENTS:
C
C  R      THIS IS THE ARRAY WHICH IS TO BE FIELLED WITH THE FIELD.
C         IT HAS DIMENSION (M,M) WHERE M IS IN COMMON/PARM/ AND IS
C         LESS THAN OR EQUAL TO 65
C
C  K      THE LEVEL FOR WHICH THE FIELD IS TO BE COMPUTED.
C
C  T      THE TIME FOR WHICH THE FIELD IS TO BE COMPUTED (NON-
C         DIMENSIONAL).
C
C  ISORV  PARAMETER WHICH INDICATES:
C
C              1        STREAMFUNCTION  
C              2        VORTICITY
C              3        PERTURBATION STREAMFUNCTION (EQUIV. TO 1, HERE)
C              4          "      "   VORTICITY (EQUIV. TO 2, HERE)
c              5        top vertical velocity
c              7        top density
C
C  IFBND  PARAMETER WHICH INDICATES:
C      
C              0        COMPUTED FULL FIFLD.
C         OTHERWISE     COMPUTE BOUNDARY ONLY.
C
C
C
C  ADDITIONAL PARAMETERS PASSED THROUGH COMMON/USER/PRAM
C
C  PRAM(1)  MODE OF THE VERTICAL OPERATOR TO USE.
C  
C      (2)  SCALED VELOCITY AMPLITUDE AT THE UPPER LEVEL.
C           (NONDIMENSIONAL)
C
C      (3)  WAVELENGTH. (NONDIMENSIONAL)
C
C      (4)  PROGATION ANGLE OF THE WAVE.
C
C      (5)  ADVECTION, IN DIMENSIONLESS UNITS.
C
C
C---------------------------------------------------------------
C
C
C
c     INCLUDE 'qg_par.inc'
      INCLUDE 'mladv.inc'
C
      parameter (maxmode = 4)
      parameter (pi = 3.141592654)
      real ch,zn
      dimension r(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/SOLVER/VAM(M0),VBM(M0),VAN(N0),VBN(N0)
      COMMON/EDGES/ICORN(4),NCORN(4),MAST(4),NAST(4)
      COMMON/THREED/HZ(K0),SCSIG,SIGZ(K01),EIGVAL(K0),AINV(K0,K0),
     +      AMAT(K0,K0)
      COMMON/USER/PRAM(16)
      COMMON/DYNAM/ALPHA,BETA,RF
      COMMON/UNIT/IIN,IOUT,IDFF,IDIF,IORR,IOWR,IPLT,IFNDD
      COMMON/CNN/N,NM1,NM2,NP1
      COMMON/EXTPAR/RLAT0,RLNG0,T0,V0,DHOR,HT,TIME0,
     +                R0,F0STAR
      common/control/ddt(4),adt(4),nalev(4),ialev(4,k02),iprnt(k02),
     +      lprnt(k0)
      common/fields/vort(mn0,k0),stream(mn0,k0),prvjac(mn0,k0),
     +     topden(mn0),botden(mn0),relief(mn0),prvtop(mn0),
     +     prvbot(mn0),topfct(k0),botfct(k0),pvort(mn0)
      common/asspar/isd,ibeg,iasday,iastep,ifsdo,imf,icit
      common/assfield/ch3(mn0),ch4(mn0,k0),ch5(mn0,k0),
     +     bnd3(mn0),bnd4(mn0,k0),bnd5(mn0,k0)
      common/unit1/isdo,iref,istc,istr,ivoc,ivor,isbc,isbr,
     +              ivbc,ivbr
c

   99 CONTINUE

        zn=iasday
        ch=isd
  105   if(icalc.eq.0)ch=0

      IF(ISORV .NE. 1 .AND. ISORV .NE. 3) GOTO 199
  100 CONTINUE
C
C --- COMPUTE STREAMFUNCTION.
C
      IF(IFBND.NE.0) GO TO 149
C
C --- FULL FIELD COMPUTATION.
C
         IP = 0
         DO 120 J = 1, N
           DO 110 I = 1, M
             IP = IP + 1
             r(ip) = ch4(ip,k)
  110       continue
  120     continue
  130  continue

      return

  149 continue
C
C --- FOR BOUNDARY COMPUTATION

         IP = 0
         DO 170 J = 1, N, NM1
            DO 160 I = 1, M
               IP = IP + 1
               r(ip) = ch4(ip,k)+ch*(bnd4(ip,k)-ch4(ip,k))/zn
  160        CONTINUE
           IP = MSQ - M
  170     CONTINUE
C
         IP = 1
         DO 190 I = 1, M, MM1
            DO 180 J = 2, NM1
               IP = IP + M  
               r(ip) = ch4(ip,k)+ch*(bnd4(ip,k)-ch4(ip,k))/zn

  180        CONTINUE
            IP = M
  190     CONTINUE

  195  CONTINUE

      RETURN

  199 CONTINUE

      IF(ISORV .NE. 2 .AND. ISORV .NE. 4) GO TO 299

  200 CONTINUE
C
C --- VORTICITY COMPUTATION.
C
      IF(IFBND .NE. 0) GO TO 249
C
C --- FULL FIELD COMPUTATION.
C
         IP = 0
         DO 220 J = 1, N
            DO 210 I = 1, M
               IP = IP + 1
               r(ip) = ch5(ip,k)
  210        CONTINUE
  220     CONTINUE
  230  CONTINUE
       RETURN

  249 CONTINUE
C
C --- FOR BOUNDARY COMPUTATION.
C
c          call proect(r,k,1)
c       return

         IP = 0
         DO 270 J = 1, N, NM1
            DO 260 I = 1, M
               IP = IP + 1
               r(ip) = ch5(ip,k)+ch*(bnd5(ip,k)-ch5(ip,k))/zn
  260        CONTINUE
            IP  = MSQ - M
  270     CONTINUE

         IP = 1
         DO 290 I = 1,M,MM1
            DO 280 J = 2, NM1
               IP = IP + M
               r(ip) = ch5(ip,k)+ch*(bnd5(ip,k)-ch5(ip,k))/zn
  280        CONTINUE
            IP = M
  290     CONTINUE
  295  CONTINUE

       RETURN

  299  CONTINUE

      if (isorv.ne.5) goto 599
c
c     top vertical velocity
c
      call curl

      fact = dhor / (v0 * ht)

      do 500 ip=1,msq
        r(ip) = fact * wsurf(ip)
  500  continue

       return

  599 continue
      if (isorv.ne.7) go to 899
c
c     top density
        if (ifbnd.ne.0) go to 740
c
c
         ip = 0
         do 720 j = 1, n
           do 710 i = 1, m
             ip = ip + 1
             r(ip) = ch3(ip)
  710      continue
  720    continue
  730  continue

       return
C
C --- FOR BOUNDARY COMPUTATION.
C
  740    continue
c          call proect(r,1,0)
c         return

         IP=0
         DO 770 J = 1, N, NM1
            DO 760 I = 1, M
               IP = IP + 1
               r(ip) = ch3(ip)+ch*(bnd3(ip)-ch3(ip))/zn
  760        CONTINUE
            IP  = MSQ - M
  770     CONTINUE

         IP = 1
         DO 790 I = 1,M,MM1
            DO 780 J = 2, NM1
               IP = IP + M
               r(ip) = ch3(ip)+ch*(bnd3(ip)-ch3(ip))/zn
  780        CONTINUE
            IP = M
  790     CONTINUE
        CONTINUE

        RETURN


  899 continue
      if(isorv.ne.9) go to 998
      return

  998 if(isorv.ne.10) goto 999
c --------      call mladv
      return
c
c
c
  999 continue
c
c --- if here, then error.
c
      write(iout,1001)isorv
 1001 format('  infld called with invalid isorv:',i5)
      stop
      end
