!---------------------------------------------------------------------
!  Fortran program to calculate isarobaric Ertel's Potential Vorticity 
!  at 500, 250, and 200 hPa level 
!---------------------------------------------------------------------
!  At first, create GrADS data for required variables (by ARWpost.exe). 
!  plot, fields, and interp_levels in namelist.ARWpost  
!  for this case are 
!
!  plot='list'
!  output_root_name = './domain_1_PV_input'
!  fields = 'XLAT,XLONG,U,V,W,T,P,PB,F'
!  interp_levels = 510.,500.,490.,260.,250.,240.,210.,200.,190., 
!-----------------------------------------------------------------------
! Then compile and run this program 
! compile: ifort -convert big_endian -o a.out calc-EPV.f 
! run:    ./a.out 
!------------------------------------------------------------------
! Output data: 
! EPV500: PV at 500 hPa (PVU, 1PVU=10^{-6} m^2 s^{-1} K kg^{-1} ) 
! EPV250: PV at 250 hPa 
! EPV200: PV at 250 hPa
!-------------------------------------------------------------------
      PARAMETER (IM=210,JM=240,KM=9,LMD=19)
!     (IM,JM)=> horizontal dimension of input data
!      LMD = dimension with respect to time   
!   ORIGINAL DATA (INPUT DATA = list in namelist.ARWpost) 

C  input grads data
      REAL  XLAT(IM,JM)   ! latitude
      REAL  XLONG(IM,JM)  ! longitude
      REAL  U(IM,JM,KM)   ! U 
      REAL  V(IM,JM,KM)   ! V 
      REAL  T(IM,JM,KM)   ! T 
      REAL  F(IM,JM)      ! Coriolis parameter

! intermediate data
      REAL P(KM)   ! Pressure (Pa)
      REAL PI(KM)  ! Exner function

      REAL DTDPI500(IM,JM),DTDPI250(IM,JM),DTDPI200(IM,JM)
      REAL DUDPI500(IM,JM),DUDPI250(IM,JM),DUDPI200(IM,JM)
      REAL DVDPI500(IM,JM),DVDPI250(IM,JM),DVDPI200(IM,JM)
      REAL DTDX500(IM,JM),DTDX250(IM,JM),DTDX200(IM,JM)
      REAL DTDY500(IM,JM),DTDY250(IM,JM),DTDY200(IM,JM)
      REAL AVOR500(IM,JM),AVOR250(IM,JM),AVOR200(IM,JM)

! Output data 
      REAL EPV500(IM,JM) 
      REAL EPV250(IM,JM)   
      REAL EPV200(IM,JM)   

      INTEGER ID,JD,I,J,K,JY,M,L

!******************
      DX=18000.      ! grid spacing 
      DY=18000.      ! grid spacing 
!*******************

      irecin=0
      irec=0

! open input data file 
      open (10,file='domain_1_PV_input.dat'
     *       ,form='unformatted',access='direct',recl=im)


! open output data file 
      open (20,file='domain_1_PV_output.dat',
     *          form='unformatted',access='direct',recl=im)


      do 3000 L=1,LMD

      write(6,*)'L=',L
      write(6,*)'READ INPUT DATA'


      do 200 j=1,jm
      irecin=irecin+1   
      read(10,rec=irecin)  (XLAT(i,j),i=1,im)
 200  continue

      do 202 j=1,jm
      irecin=irecin+1   
      read(10,rec=irecin)  (XLONG(i,j),i=1,im)
 202  continue

      do 205 k=1,km
      do 205 j=1,jm
      irecin=irecin+1   
      read(10,rec=irecin)  (U(i,j,k),i=1,im)
 205  continue

      do 210 k=1,km
      do 210 j=1,jm
      irecin=irecin+1   
      read(10,rec=irecin)  (V(i,j,k),i=1,im)
 210  continue

      do 220 k=1,km
      do 220 j=1,jm
      irecin=irecin+1   
      read(10,rec=irecin)  (T(i,j,k),i=1,im)
 220  continue

      do 280 j=1,jm
      irecin=irecin+1   
      read(10,rec=irecin)  (F(i,j),i=1,im)
 280  continue

!  END  READ ORIGINAL DATA  

      RD=287.0    ! Gas constant for dry air 
      CP=1005.7   ! Heat capacity at constant pressure 
      G=9.81

      P(1)=51000.
      P(2)=50000.
      P(3)=49000.
      P(4)=26000.
      P(5)=25000.
      P(6)=24000.
      P(7)=21000.
      P(8)=20000.
      P(9)=18000.

      do 300 k=1,km
      PI(k)= (P(k)/101300.)**(RD/CP)
      write(6,*) PI(k)
 300  continue




      do 400 j=2,jm-1
      do 400 i=2,im-1


! Calculate intermediate data at 500, 250, 200hPa 
      DTDPI500(i,j)=(T(i,j,3)-T(I,j,1))/(PI(3)-PI(1))
      DTDPI250(i,j)=(T(i,j,6)-T(I,j,4))/(PI(6)-PI(4))
      DTDPI200(i,j)=(T(i,j,9)-T(I,j,7))/(PI(9)-PI(7))

      DUDPI500(i,j)=(U(i,j,3)-U(I,j,1))/(PI(3)-PI(1))
      DUDPI250(i,j)=(U(i,j,6)-U(I,j,4))/(PI(6)-PI(4))
      DUDPI200(i,j)=(U(i,j,9)-U(I,j,7))/(PI(9)-PI(7))

      DVDPI500(i,j)=(V(i,j,3)-V(I,j,1))/(PI(3)-PI(1))
      DVDPI250(i,j)=(V(i,j,6)-V(I,j,4))/(PI(6)-PI(4))
      DVDPI200(i,j)=(V(i,j,9)-V(I,j,7))/(PI(9)-PI(7))

      DTDX500(i,j)=(T(i+1,j,2)-T(I-1,j,2))/2.0/DX
      DTDX250(i,j)=(T(i+1,j,5)-T(I-1,j,5))/2.0/DX
      DTDX200(i,j)=(T(i+1,j,8)-T(I-1,j,8))/2.0/DX

      DTDY500(i,j)=(T(i,j+1,2)-T(I,j-1,2))/2.0/DY
      DTDY250(i,j)=(T(i,j+1,5)-T(I,j-1,5))/2.0/DY
      DTDY200(i,j)=(T(i,j+1,8)-T(I,j-1,8))/2.0/DY


! Absolute vorticty at 500, 250, 200 hPa

      AVOR500(i,j)=F(i,j)+(V(i+1,j,2)-V(i-1,j,2))/2.0/DX
     *                   -(U(i,j+1,2)-U(i,j-1,2))/2.0/DY

      AVOR250(i,j)=F(i,j)+(V(i+1,j,5)-V(i-1,j,5))/2.0/DX
     *                   -(U(i,j+1,5)-U(i,j-1,5))/2.0/DY

      AVOR200(i,j)=F(i,j)+(V(i+1,j,8)-V(i-1,j,8))/2.0/DX
     *                   -(U(i,j+1,8)-U(i,j-1,8))/2.0/DY

 400  continue



        do 450 j=1,jm
        do 450 i=1,im
        EPV500(i,j)=1.e30
        EPV250(i,j)=1.e30
        EPV200(i,j)=1.e30
 450    CONTINUE

! Calculate isarobaric EPV 

        do 500 j=2,jm-1
        do 500 i=2,im-1

        EPV500(i,j)=-g*RD/CP*PI(2)/p(2)*1e6
     *       *(AVOR500(i,j)*DTDPI500(i,j)
     *         -DVDPI500(i,j)*DTDX500(i,j)
     *         +DUDPI500(i,j)*DTDY500(i,j) )

        EPV250(i,j)=-g*RD/CP*PI(5)/p(5)*1e6
     *       *(AVOR250(i,j)*DTDPI250(i,j)
     *         -DVDPI250(i,j)*DTDX250(i,j)
     *         +DUDPI250(i,j)*DTDY250(i,j) )

        EPV200(i,j)=-g*RD/CP*PI(8)/p(8)*1e6
     *       *(AVOR200(i,j)*DTDPI200(i,j)
     *         -DVDPI200(i,j)*DTDX200(i,j)
     *         +DUDPI200(i,j)*DTDY200(i,j) )

 500    CONTINUE


!**** Write interpolated data and calculated data

      do 1005 j=1,jm
      irec=irec+1
 1005 write(20,rec=irec)  (XLAT(i,j),i=1,im)

      do 1006 j=1,jm
      irec=irec+1
 1006 write(20,rec=irec)  (XLONG(i,j),i=1,im)

      do 1010 j=1,jm
      irec=irec+1
 1010 write(20,rec=irec)  (EPV500(i,j),i=1,im)

      do 1020 j=1,jm
      irec=irec+1
 1020 write(20,rec=irec)  (EPV250(i,j),i=1,im)

      do 1100 j=1,jm
      irec=irec+1
 1100 write(20,rec=irec)  (EPV200(i,j),i=1,im)

 2999 continue
 3000 continue

      stop
      end 
