      PROGRAM EPFLUX_CMIP6
C     PROGRAM EPFLUX_CMIP6 (UPVP, UPWP, VPTP, UZ, TZ, VZ, WZ, 
C    1                EPFY, EPFZ, EPFD, VRES, WRES, 
C    B                UTENDVTEM, UTENDWTEM, OUTPUT  )
C    2                TAPE11=UPVP, TAPE12=UPWP, TAPE13=VPTP, 
C    3                TAPE14=UZ,   TAPE15=TZ,   TAPE16=VZ,   TAPE17=WZ, 
C    4                TAPE21=EPFY, TAPE22=EPFZ, TAPE23=EPFD, 
C    5                TAPE24=VRES, TAPE25=WRES, 
C    C                TAPE26=UTENDVTEM, TAPE27=UTENDWTEM, TAPE6=OUTPUT)                                                                
C     --------------------------------------------------------------
C
C DEC 20/06 - C. MCLANDRESS
C 
C AUTHOR  - C. MCLANDRESS
C
C MODIFICATIONS:
C
C           FEB 14/2018 - DIVIDE EP FLUXES BY REFERENCE DENSITY SO THAT THEY 
C                         HAVE THE UNITS REQUIRED FOR CMIP6. CHANGE CONSTANTS
C                         TO BE CONSISTENT WITH CMIP6.  
C                         RENAMED PROGRAM EPFLUX2 --> EPFLUX_CMIP6-- M. SIGMOND/C.MCLANDRESS
C
C           NOV 15/2018 - COMPUTE HORIZONTAL AND VERTICAL ADVECTIVE TENDENCY 
C                         TERMS IN TEM EQUATION: VRES*(D(UZ*COSY)/DY/ACOSY-F)  
C                         AND WRES*D(UZ)/DZ -- C.MCLANDRESS
C
C PURPOSE - COMPUTES EP FLUX DIAGNOSTICS AND RESIDUAL CIRCULATION
C           IN LOG-PRESSURE COORDINATES USING THE FORMULATION GIVEN IN 
C           "MIDDLE ATMOSPHERE DYNAMICS" (1987, P.128) BY ANDREWS, HOLTON & LEOVY.
C           THE INPUT FLUXES CAN EITHER BE TIME-VARYING (IE INSTANTANEOUS)
C           OR STATIONARY FLUXES COMPUTED FROM THE TIME-MEAN WAVE FIELDS.
C
C INPUT FILES:
C         UPVP - MERIDIONAL FLUX OF ZONAL MOMENTUM FOR ZONAL DEVIATIONS (M^2/S^2)
C         UPWP - VERTICAL (IE PRESSURE VELOCITY) FLUX OF ZONAL MOMENTUM FOR 
C                ZONAL DEVIATIONS (PA M/S^2)
C         VPTP - MERIDIONAL HEAT FLUX FOR ZONAL DEVIATIONS (K M/S)
C         UZ   - ZONAL MEAN ZONAL WIND (M/S)
C         TZ   - ZONAL MEAN TEMPERATURE (K)
C         VZ   - ZONAL MEAN MERIDIONAL WIND (M/S)
C         WZ   - ZONAL MEAN VERTICAL WIND (PRESSURE VELOCITY IN PA/S)
C
C OUTPUT FILES:
C         EPFY - MERIDIONAL COMPONENT OF EP FLUX VECTOR (M^3/S^2)
C         EPFZ - VERTICAL COMPONENT OF EP FLUX VECTOR (M^3/S^2)
C         EPFD - EP FLUX DIVERGENCE (FORCE PER UNIT MASS, 
C                IE DIVIDED BY RHO*A*COS(PHI)) (M/S^2)
C         VRES - RESIDUAL MERIDIONAL VELOCITY (M/S) 
C         WRES - RESIDUAL VERTICAL VELOCITY (M/S) 
C         UTENDVTEM - TENDENCY OF EASTWARD WIND DUE TO TEM NORTHWARD WIND 
C                     ADVECTION AND THE CORIOLIS TERM)
C         UTENDWTEM (TENDENCY OF EASTWARD WIND DUE TO TEM UPWARD WIND ADVECTION)
C----------------------------------------------------------------------------
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      PARAMETER (MAXLAT=130, MAXLEV=100)
C
      LOGICAL OK, OK1, OK2, OK3, OK4, OK5, OK6, OK7
      REAL*8 SINY(MAXLAT),COSY(MAXLAT),RAD(MAXLAT)
      REAL*8 WL(MAXLAT),WOSSL(MAXLAT)
      INTEGER LEV(MAXLEV)                                                              
      REAL PR(MAXLEV), ZLOGP(MAXLEV), DENS(MAXLEV)
      REAL TFAC(MAXLEV), WFAC(MAXLEV)
      REAL UPVP(MAXLAT,MAXLEV), UPWP(MAXLAT,MAXLEV)
      REAL VPTP(MAXLAT,MAXLEV)
      REAL UZ(MAXLAT,MAXLEV),     TZ(MAXLAT,MAXLEV)
      REAL VZ(MAXLAT,MAXLEV),     WZ(MAXLAT,MAXLEV)
      REAL UZDZ(MAXLAT,MAXLEV), TZDZ(MAXLAT,MAXLEV)
      REAL EPFY(MAXLAT,MAXLEV), EPFZ(MAXLAT,MAXLEV)
      REAL EPFD(MAXLAT,MAXLEV)
      REAL VRES(MAXLAT,MAXLEV), WRES(MAXLAT,MAXLEV)
      REAL UTENDVTEM(MAXLAT,MAXLEV)
      REAL UTENDWTEM(MAXLAT,MAXLEV)
      REAL WORKV(MAXLAT,MAXLEV),WORKW(MAXLAT,MAXLEV)
      REAL FHAT(MAXLAT,MAXLEV)
      COMMON/ICOM/IBUF(8),IDAT(MAXLAT)
      PARAMETER (OMEGA=7.29212E-5, RADIUS=6.37123E6)
      PARAMETER (PR0=1013.25E2,  HSCALE=7.E3)
      PARAMETER (RGAS=287.058,  CPGAS=1004.64E0, GRAVITY=9.80665 )
      PARAMETER (RKAPPA=RGAS/CPGAS, DENS0=PR0/(GRAVITY*HSCALE))
C---------------------------------------------------------------------
      NF=14 
      CALL JCLPNT(NF,11,12,13,14,15,16,17,21,22,23,24,25,26,27,6)
C
C     * REWIND INPUT FILES.
C
      DO 100 N=11,17
 100  REWIND N
C
C     * DETERMINE THE PRESSURE LEVELS OF FIRST INPUT FILE - IT IS ASSUMED
C     * THAT OTHER FILES ARE THE SAME.
C     
      CALL FILEV(LEV,NLEV,IBUF,11)
C
      NLAT=IBUF(5)
      NLATM=NLAT-1
      NLEVM=NLEV-1
      PI=2.*ASIN(1.E0)
      RADSP=-PI/2.
      RADNP=PI/2.
      OM2=2.*OMEGA
C
C     * ABORT IF ARRAY SIZES OR TYPE INCORRECT.
C
      IF((NLEV.LT.3).OR.(NLEV.GT.MAXLEV)) CALL   XIT('EPFLUX2',-1) 
      IF(NLAT.GT.MAXLAT) CALL                    XIT('EPFLUX2',-2) 
      IF(IBUF(1).NE.4HZONL) CALL                 XIT('EPFLUX2',-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)
C
C     * COMPUTE PRESSURE IN PASCALS,
C     * LOG-PRESSURE HEIGHT AND SCALE FACTORS FOR CONVERTING 
C     * TO POTENTIAL TEMPERATURE AND VERTICAL VELOCITY.
C
      CALL LVDCODE(PR,LEV,NLEV)
      DO L=1,NLEV
        PR(L)=PR(L)*100.E0 
        ZLOGP(L)=-HSCALE*ALOG(PR(L)/PR0)                
        TFAC(L)=(PR0/PR(L))**RKAPPA
        WFAC(L)=-HSCALE/PR(L)
        DENS(L)=DENS0*EXP(-ZLOGP(L)/HSCALE)
      END DO
c ----------------------------------------------------------------
c      write(6,*) 'nlev=',nlev
c      do l=1,nlev
c        write(6,6100) l,zlogp(l)/1.e3,pr(l),
c     &                dens(l),wfac(l),tfac(l)
c      end do
c 6100 format('l,z,p,dens,wfac,tfac=',i2,f6.1,4(e12.5,1x))
c      write(6,*) 
c      write(6,*) 'nlat=',nlat
c      write(6,*) pi
c      do j=1,nlat
c        ydeg=rad(j)*180./pi
c        write(6,6110) j,ydeg,rad(j),cosy(j),siny(j)
c      end do
c 6110 format('j,y,rad,cos,sin=',i2,f8.2,3(1x,f6.3))
c      write(6,*) 
c ----------------------------------------------------------------
C
      NR=0
      NSETS=0
C     
C     * READ THE NEXT SETS OF FIELDS.
C
 1000 CONTINUE 
C
C     * LOOP THROUGH LEVELS AND PUT DATA INTO LAT VS HEIGHT ARRAYS,
C     * CONVERT TEMPERATURE TO POT TEMP, AND OMEGA TO LOG-PRESSURE
C     * VERTICAL VELOCITY.
C
      DO L=1,NLEV
        CALL GETFLD2(11,UPVP(1,L),NC4TO8("ZONL"),-1,-1,-1,
     &                  IBUF,MAXLAT,OK1) 
        CALL GETFLD2(12,UPWP(1,L),NC4TO8("ZONL"),-1,-1,-1,
     &                  IBUF,MAXLAT,OK2) 
        CALL GETFLD2(13,VPTP(1,L),NC4TO8("ZONL"),-1,-1,-1,
     &                  IBUF,MAXLAT,OK3) 
        CALL GETFLD2(14,UZ(1,L),NC4TO8("ZONL"),-1,-1,-1,
     &                  IBUF,MAXLAT,OK4)         
        CALL GETFLD2(15,TZ(1,L),NC4TO8("ZONL"),-1,-1,-1,
     &                  IBUF,MAXLAT,OK5)         
        CALL GETFLD2(16,VZ(1,L),NC4TO8("ZONL"),-1,-1,-1,
     &                  IBUF,MAXLAT,OK6)         
        CALL GETFLD2(17,WZ(1,L),NC4TO8("ZONL"),-1,-1,-1,
     &                  IBUF,MAXLAT,OK7)         
        OK=.FALSE.
        IF(OK1.AND.OK2.AND.OK3.AND.OK4.AND.OK5.AND.OK6.AND.OK7) 
     &  OK=.TRUE.
        IF(.NOT.OK) THEN
          WRITE(6,6010) NSETS 
          IF(NR.EQ.0)  CALL                         XIT('EPFLUX2',-4)  
          CALL                                      XIT('EPFLUX2',0)  
        END IF
        DO J=1,NLAT
          UPWP(J,L)=UPWP(J,L)*WFAC(L)
          VPTP(J,L)=VPTP(J,L)*TFAC(L)
          TZ(J,L)=TZ(J,L)*TFAC(L)
          WZ(J,L)=WZ(J,L)*WFAC(L)
        END DO
        NR=NR+1
      END DO
c ----------------------------------------------------------------
c      do l=1,nlev
c        write(6,6120) l,zlogp(l)*1.e-3,uz(10,l),vz(10,l),wz(1,l)
c      enddo
c6120  format(' l,z,u,v,w=',i3,f6.1,1x,f7.2,2(1x,e14.7))
c      write(6,*) 
c ----------------------------------------------------------------
C
C     * COMPUTE VERTICAL DERIVATIVES. ENDS USE ONE-SIDED DIFFERENCE;
C     * INTERIOR USES CENTERED DIFF.                                                                 
C
      DO J=1,NLAT
        DZ=ZLOGP(1)-ZLOGP(2)
        UZDZ(J,1)=(UZ(J,1)-UZ(J,2))/DZ
        TZDZ(J,1)=(TZ(J,1)-TZ(J,2))/DZ
        DO L=2,NLEVM 
          DZ=ZLOGP(L-1)-ZLOGP(L+1)
          UZDZ(J,L)=(UZ(J,L-1)-UZ(J,L+1))/DZ
          TZDZ(J,L)=(TZ(J,L-1)-TZ(J,L+1))/DZ
        END DO
        DZ=ZLOGP(NLEVM)-ZLOGP(NLEV)
        UZDZ(J,NLEV)=(UZ(J,NLEVM)-UZ(J,NLEV))/DZ
        TZDZ(J,NLEV)=(TZ(J,NLEVM)-TZ(J,NLEV))/DZ
      END DO
C
C     * WORK ARRAYS FOR RESIDUAL CIRCULATION
C
      DO J=1,NLAT
        DO L=1,NLEV
          WORKV(J,L)=DENS(L)*VPTP(J,L)/TZDZ(J,L)
          WORKW(J,L)=COSY(J)*VPTP(J,L)/TZDZ(J,L)
        END DO
      ENDDO                                                                                               
      DO J=1,NLAT
        DZ=ZLOGP(1)-ZLOGP(2)
        VRES(J,1)=(WORKV(J,1)-WORKV(J,2))/DZ
        DO L=2,NLEVM 
          DZ=ZLOGP(L-1)-ZLOGP(L+1)
          VRES(J,L)=(WORKV(J,L-1)-WORKV(J,L+1))/DZ
        END DO
        DZ=ZLOGP(NLEVM)-ZLOGP(NLEV)
        VRES(J,NLEV)=(WORKV(J,NLEVM)-WORKV(J,NLEV))/DZ
      END DO
c ----------------------------------------------------------------
c      do j=1,nlat
c        do l=1,nlev
c          uz(j,l)=omega*radius*cosy(j)
c        enddo
c      enddo
c ----------------------------------------------------------------
C
C     * COMPUTE LATITUDINAL DERIVATIVES. USE CENTERED
C     * DIFFERENCES AT ALL POINTS (FUNCTION ZERO AT POLES).
C
      DO L=1,NLEV
        DY=RAD(2)-RADSP
        FHAT(1,L)=UZ(2,L)*COSY(2)/DY
c        FHAT(1,L)=UZ(2,L)/DY
        WRES(1,L)=WORKW(2,L)/DY
        DO J=2,NLATM
          DY=RAD(J+1)-RAD(J-1)
          FHAT(J,L)=(UZ(J+1,L)*COSY(J+1)-UZ(J-1,L)*COSY(J-1))/DY
c          FHAT(J,L)=(UZ(J+1,L)-UZ(J-1,L))/DY
          WRES(J,L)=(WORKW(J+1,L)-WORKW(J-1,L))/DY
        END DO
        DY=RADNP-RAD(NLATM)
        FHAT(NLAT,L)=-UZ(NLATM,L)*COSY(NLATM)/DY
c        FHAT(NLAT,L)=-UZ(NLATM,L)/DY
        WRES(NLAT,L)=-WORKW(NLATM,L)/DY
      ENDDO                                                                                               
C
      DO J=1,NLAT
        ACOS=RADIUS*COSY(J)
        CORIOL=OM2*SINY(J)
        DO L=1,NLEV
          FHAT(J,L)=CORIOL-FHAT(J,L)/ACOS
c          FHAT(J,L)=CORIOL+UZ(J,L)*TAN(RAD(J))/RADIUS-FHAT(J,L)/RADIUS
          VRES(J,L)=VZ(J,L)-VRES(J,L)/DENS(L)
          WRES(J,L)=WZ(J,L)+WRES(J,L)/ACOS
        ENDDO                                                                                               
      ENDDO                                                                                               
c ----------------------------------------------------------------
c      do l=1,nlev
c         write(6,6211) l,zlogp(l)*1.e-3,(fhat(j,l),j=1,nlat)
c      enddo
c 6211 format(' l,z,fhat=',i3,f6.1,1x,32(e10.3,1x))     
c ----------------------------------------------------------------
C
C     * COMPUTE COMPONENTS OF EP FLUX VECTOR
C
      DO J=1,NLAT
        ACOS=RADIUS*COSY(J)
        DO L=1,NLEV
          AFCT=DENS(L)*ACOS
          EPFY(J,L)=AFCT*(UZDZ(J,L)*VPTP(J,L)/TZDZ(J,L)-UPVP(J,L))
          EPFZ(J,L)=AFCT*(FHAT(J,L)*VPTP(J,L)/TZDZ(J,L)-UPWP(J,L)) 
        END DO
      ENDDO                                                                                               
C
C     * EP FLUX DIVERGENCE    
C
      DO J=1,NLAT
        DZ=ZLOGP(1)-ZLOGP(2)
        EPFD(J,1)=(EPFZ(J,1)-EPFZ(J,2))/DZ
        DO L=2,NLEVM 
          DZ=ZLOGP(L-1)-ZLOGP(L+1)
          EPFD(J,L)=(EPFZ(J,L-1)-EPFZ(J,L+1))/DZ
        END DO
        DZ=ZLOGP(NLEVM)-ZLOGP(NLEV)
        EPFD(J,NLEV)=(EPFZ(J,NLEVM)-EPFZ(J,NLEV))/DZ
      END DO
      DO L=1,NLEV
        ACOSDY=RADIUS*COSY(1)*(RAD(2)-RADSP)
        EPFD(1,L)=EPFD(1,L) +
     &        EPFY(2,L)*COSY(2)/ACOSDY
        DO J=2,NLATM
          ACOSDY=RADIUS*COSY(J)*(RAD(J+1)-RAD(J-1))
          EPFD(J,L)=EPFD(J,L) +
     &       (EPFY(J+1,L)*COSY(J+1)-EPFY(J-1,L)*COSY(J-1))/ACOSDY
        END DO
        ACOSDY=RADIUS*COSY(NLAT)*(RADNP-RAD(NLATM))
        EPFD(NLAT,L)=EPFD(NLAT,L)
     &       -EPFY(NLATM,L)*COSY(NLATM)/ACOSDY
      ENDDO
C
C     * CONVERT EP FLUX DIVERGENCE TO ZONAL WIND TENDENCY AND DIVIDE EP FLUXES
C     * BY REFERENCE DENSITY SO THAT THEY HAVE THE UNITS REQUIRED FOR CMIP6
      DO J=1,NLAT
        ACOS=RADIUS*COSY(J)
        DO L=1,NLEV
          EPFD(J,L)=EPFD(J,L)/(DENS(L)*ACOS)
          EPFY(J,L)=EPFY(J,L)/DENS(L)
          EPFZ(J,L)=EPFZ(J,L)/DENS(L)
        ENDDO
      ENDDO
C
C     ADVECTIVE TENDENCY TERMS IN TEM EQN
C
      DO J=1,NLAT
        DO L=1,NLEV
           UTENDVTEM(J,L)=VRES(J,L)*FHAT(J,L)
           UTENDWTEM(J,L)=WRES(J,L)*UZDZ(J,L)
        ENDDO
      ENDDO
C
c ----------------------------------------------------------------
c      do l=1,nlev
c        write(6,6150) l,zlogp(l)*1.e-3,vres(10,l),wres(1,l),
c     &                epfd(10,l)*86400.
c      enddo
c6150  format(' l,z,vr,wr,epfd=',i3,f6.1,1x,3(e14.7,1x))
c      write(6,*) 
c ----------------------------------------------------------------
C
C     * PUT THE RESULTS INTO OUPUT FILES.
C
      DO L=1,NLEV
        CALL SETLAB (IBUF,-1,-1,NC4TO8("EPFY"),LEV(L),-1,-1,-1,-1) 
        CALL PUTFLD2(21,EPFY(1,L),IBUF,MAXLAT)
        CALL SETLAB (IBUF,-1,-1,NC4TO8("EPFZ"),LEV(L),-1,-1,-1,-1) 
        CALL PUTFLD2(22,EPFZ(1,L),IBUF,MAXLAT) 
        CALL SETLAB (IBUF,-1,-1,NC4TO8("EPFD"),LEV(L),-1,-1,-1,-1) 
        CALL PUTFLD2(23,EPFD(1,L),IBUF,MAXLAT)
        CALL SETLAB (IBUF,-1,-1,NC4TO8("VRES"),LEV(L),-1,-1,-1,-1)  
        CALL PUTFLD2(24,VRES(1,L),IBUF,MAXLAT)
        CALL SETLAB (IBUF,-1,-1,NC4TO8("WRES"),LEV(L),-1,-1,-1,-1)  
        CALL PUTFLD2(25,WRES(1,L),IBUF,MAXLAT)
        CALL SETLAB (IBUF,-1,-1,NC4TO8("UTHA"),LEV(L),-1,-1,-1,-1)  
        CALL PUTFLD2(26,UTENDVTEM(1,L),IBUF,MAXLAT)
        CALL SETLAB (IBUF,-1,-1,NC4TO8("UTVA"),LEV(L),-1,-1,-1,-1)  
        CALL PUTFLD2(27,UTENDWTEM(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
