      PROGRAM MSTRFCN
C     PROGRAM MSTRFCN ( VRES, MPSI, OUTPUT  )
C    2                  TAPE1=VRES, TAPE2=MPSI, TAPE6=OUTPUT)                                                                
C     --------------------------------------------------------------
C
C OCT 6/15 - C. MCLANDRESS
C 
C AUTHOR  - C. MCLANDRESS
C
C PURPOSE - COMPUTES MASS STREAMFUNCTION FROM RESIDUAL MERIDIONAL VELOCITY
C           USING EQN (2) OF MCLANDRESS AND SHEPHERD (JCLI 2009, PP 1516-1540)
C
C INPUT FILE:
C         VRES - RESIDUAL MERIDIONAL VELOCITY (M/S) 
C
C OUTPUT FILE:
C         MPSI - MASS STREAMFUNCTION (KG/M/S)
C----------------------------------------------------------------------------
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      PARAMETER (MAXLAT=100, MAXLEV=100)
C
      LOGICAL OK
      REAL*8 SINY(MAXLAT),COSY(MAXLAT),RAD(MAXLAT)
      REAL*8 WL(MAXLAT),WOSSL(MAXLAT),YFAC(MAXLAT)
      INTEGER LEV(MAXLEV)                                                              
      REAL PR(MAXLEV), DPR(MAXLEV), ZLOGP(MAXLEV)
      REAL VRES(MAXLAT,MAXLEV),MPSI(MAXLAT,MAXLEV)
      COMMON/ICOM/IBUF(8),IDAT(MAXLAT)
      PARAMETER (PR0=1013.25E2, HSCALE=7.E3, GRAVITY=9.81)
C---------------------------------------------------------------------
      NF=2 
      CALL JCLPNT(NF,1,2,6)
C
C     * REWIND INPUT FILE.
C
      REWIND 1
C
C     * DETERMINE THE PRESSURE LEVELS.
C     
      CALL FILEV(LEV,NLEV,IBUF,1)
C
      NLAT=IBUF(5)
C
C     * ABORT IF ARRAY SIZES OR TYPE INCORRECT.
C
      IF((NLEV.LT.3).OR.(NLEV.GT.MAXLEV)) CALL   XIT('MSTRFCN',-1) 
      IF(NLAT.GT.MAXLAT) CALL                    XIT('MSTRFCN',-2) 
      IF(IBUF(1).NE.4HZONL) CALL                 XIT('MSTRFCN',-3) 
C
C     * LATITUDE TRIGONOMETRIC ARRAYS
C
      ILATH=NLAT/2
      CALL GAUSSG(ILATH,SINY,WL,COSY,RAD,WOSSL)
      CALL  TRIGL(ILATH,SINY,WL,COSY,RAD,WOSSL)
      DO J=1,NLAT
        YFAC(J)=-COSY(J)/GRAVITY
      END DO
C
C     * COMPUTE PRESSURE IN PASCALS,
C     * LOG-PRESSURE HEIGHT FOR PRINT OUT.
C
      CALL LVDCODE(PR,LEV,NLEV)
      DO L=1,NLEV
        PR(L)=PR(L)*100.E0 
        ZLOGP(L)=-HSCALE*ALOG(PR(L)/PR0)
      END DO

      DPR(1)=-PR(1)
      DO L=2,NLEV
        DPR(L)=PR(L-1)-PR(L)
      END DO

c ----------------------------------------------------------------
c      write(6,*) 'nlev,nlat=',nlev,nlat
c      write(6,*) 
c      do l=1,nlev
c        write(6,6100) l,zlogp(l)/1.e3,pr(l),dpr(l)
c      end do
c 6100 format('l,z,p,dpr=',i2,f6.1,4(e12.5,1x))
c      write(6,*) 
c      pi=2.*asin(1.e0)
c      do j=1,nlat
c        ydeg=rad(j)*180./pi
c        write(6,6110) j,ydeg,yfac(j)
c      end do
c 6110 format('j,y,yfac=',i2,f8.2,3(1x,f9.6))
c      write(6,*) 
c ----------------------------------------------------------------
C
      NR=0
      NSETS=0
C     
C     * READ THE NEXT SETS OF FIELDS.
C
 1000 CONTINUE 

c      write(6,*) 'nsets,nr=',nsets,nr
c      jprint=20
c      write(6,*) 'jprint,lat=',jprint,rad(jprint)*180./pi
C
C     * LOOP THROUGH LEVELS AND PUT DATA INTO LAT VS HEIGHT ARRAYS,
C
      DO L=1,NLEV
        CALL GETFLD2(1,VRES(1,L),NC4TO8("ZONL"),-1,-1,-1,
     &               IBUF,MAXLAT,OK) 
        IF(.NOT.OK) THEN
          WRITE(6,6010) NSETS 
          IF(NR.EQ.0)  CALL                    XIT('MSTRFCN',-4)  
          CALL                                 XIT('MSTRFCN',0)  
        END IF
        NR=NR+1
      END DO
C
C     * INTEGRATE FROM TOP OF ATMOSPHERE (WHERE P=0 & V*=0) DOWN TO LEVEL LL
C
      DO J=1,NLAT
        DO LL=1,NLEV 
          SUM = 0.5*VRES(J,1)*DPR(1)
c          if (j.eq.jprint) then 
c            write(6,*) 'll,sum,v(1),dp(1)=',ll,sum,vres(j,1),dpr(1)
c          endif
          IF (LL.GT.1) THEN 
            DO L=2,LL
              SUM = SUM + 0.5*(VRES(J,L)+VRES(J,L-1))*DPR(L)
c              if (j.eq.jprint) then 
c                write(6,*)  '  l,sum,v(l),v(l-1),dpr(l)=',
c     $                     l,sum,vres(j,l),vres(j,l-1),dpr(l)
c              endif
            end DO
          END IF
          MPSI(J,LL) = YFAC(J)*SUM
        END DO
      END DO

c ----------------------------------------------------------------
c      do l=1,nlev
c         write(6,6211) l,zlogp(l)*1.e-3,(mpsi(j,l),j=1,nlat)
c      enddo
c 6211 format(' l,z,mpsi=',i3,f6.1,1x,32(e10.3,1x))     
c ----------------------------------------------------------------
C
C     * PUT THE RESULTS INTO OUPUT FILE.
C
      DO L=1,NLEV
        CALL SETLAB (IBUF,-1,-1,NC4TO8("MPSI"),LEV(L),-1,-1,-1,-1) 
        CALL PUTFLD2(2,MPSI(1,L),IBUF,MAXLAT)
      ENDDO
C 
C     * DO NEXT SET
C
      NSETS=NSETS+1
      GO TO 1000
C---------------------------------------------------------------------
 6010 FORMAT(' ',I6,' SETS OF MULTILEVEL RECORDS')
 6025 FORMAT(' ',A4,I10,2X,A4,I10,4I6)
      END
