!---------------------------------------------------------------------
!  Fortran program to obtain te (equivalent potential temperature)
!  and tes (saturated equivalent potential temperature)  
!  for WSM6 scheme (mp_physics=6)
!---------------------------------------------------------------------
!  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'
!
!  fields = 'XLAT,XLONG,U,V,W,T,P,PB,P_HYD,QVAPOR,QCLOUD,QRAIN,QICE,
!  QSNOW,QGRAUP,H_DIABATIC,RAINNCV,tk,rh,slp,dbz,cape'
!
!   interp_levels = 0.05, 0.3, 0.55, 0.8, 1.05, 1.3, 1.55, 1.8, 2.05, 
!   2.3, 2.55,  2.8, 3.05, 3.3, 3.55, 3.8, 4.05, 4.3, 4.55, 4.8,  
!   5.05, 5.3, 5.55, 5.8, 6.05, 6.3, 6.55, 6.8, 7.05, 7.3, 7.55, 
!   7.8, 8.05, 8.3, 8.55,  8.8,  9.05, 9.3, 9.55, 9.8, 10.05, 
!-----------------------------------------------------------------------
! Then compile and run this program 
! compile: ifort -convert big_endian -o a.out calc-te_WSM6.f 
! run:    ./a.out 
!-------------------------------------------------------------------
      PARAMETER (IMD=450,JMD=450,KMD=41,LMD=3)

!     (IMD,JMD,KMD)=> dimension of original data 
!      LMD = dimension with respect to time   

!   ORIGINAL DATA (INPUT DATA = list in namelist.ARWpost) 
      REAL  XLATD(IMD,JMD)        ! XLAT
      REAL  XLONGD(IMD,JMD)       ! XLONG
      REAL  UD(IMD,JMD,KMD)   ! U 
      REAL  VD(IMD,JMD,KMD)   ! V 
      REAL  WD(IMD,JMD,KMD)   ! W 
      REAL  TD(IMD,JMD,KMD)   ! T 
      REAL  PD(IMD,JMD,KMD)   ! P 
      REAL  PBD(IMD,JMD,KMD)  ! PB 
      REAL  PHD(IMD,JMD,KMD)  ! P_HYD 
      REAL  QVD(IMD,JMD,KMD)   ! QVAPOR 
      REAL  QCD(IMD,JMD,KMD)   ! QCLOUD 
      REAL  QRD(IMD,JMD,KMD)   ! QRAIN  
      REAL  QID(IMD,JMD,KMD)   ! QICE   
      REAL  QSD(IMD,JMD,KMD)   ! QSNOW  
      REAL  QGD(IMD,JMD,KMD)   ! QGRAUP 
!!      REAL  QNRAIND(IMD,JMD,KMD) ! QNRAIN 
      REAL  DIABAD(IMD,JMD,KMD)   ! H_DIABATIC 
      REAL  RAINNCVD(IMD,JMD)     ! RAINNCV 
      REAL  TEMPD(IMD,JMD,KMD)    ! tk (temperature in K) 
      REAL  RHD(IMD,JMD,KMD)      ! rh 
      REAL  SLPD(IMD,JMD)         ! SLP 
      REAL  DBZD(IMD,JMD,KMD)     ! dbz  
      REAL  CAPED(IMD,JMD,KMD)    ! cape 



!****Calculation of equivalent potential temperature (following G. Bryan)
      REAL  QSAT(IMD,JMD,KMD)  ! Saturation mixing ratio 
      REAL  E(IMD,JMD,KMD)       ! vapor pressure
      REAL  TK2(IMD,JMD,KMD)     ! 
      REAL  TLCL(IMD,JMD,KMD)    ! temperature at the LCL 
      REAL  TE(IMD,JMD,KMD)      ! equivalent potential temperature
      REAL  ES(IMD,JMD,KMD)      ! saturation vapor pressure 
      REAL  TLCLS(IMD,JMD,KMD)   ! 
      REAL  TES(IMD,JMD,KMD)     ! saturated equivalent potential temp



      INTEGER I,J,K,L

      irecin=0
      irec=0

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


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


      do 3000 L=1,LMD

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


      do 200 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (XLATD(i,j),i=1,imd)
 200  continue

      do 202 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (XLONGD(i,j),i=1,imd)
 202  continue

      do 205 k=1,kmd
      do 205 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (UD(i,j,k),i=1,imd)
 205  continue

      do 210 k=1,kmd
      do 210 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (VD(i,j,k),i=1,imd)
 210  continue

      do 215 k=1,kmd
      do 215 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (WD(i,j,k),i=1,imd)
 215  continue

      do 220 k=1,kmd
      do 220 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (TD(i,j,k),i=1,imd)
 220  continue


      do 225 k=1,kmd
      do 225 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (PD(i,j,k),i=1,imd)
 225  continue

      do 226 k=1,kmd
      do 226 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (PBD(i,j,k),i=1,imd)
 226  continue

      do 227 k=1,kmd
      do 227 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (PHD(i,j,k),i=1,imd)
 227  continue


      do 230 k=1,kmd
      do 230 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (QVD(i,j,k),i=1,imd)
 230  continue

      do 232 k=1,kmd
      do 232 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (QCD(i,j,k),i=1,imd)
 232  continue

      do 234 k=1,kmd
      do 234 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (QRD(i,j,k),i=1,imd)
 234  continue

      do 236 k=1,kmd
      do 236 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (QID(i,j,k),i=1,imd)
 236  continue

      do 238 k=1,kmd
      do 238 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (QSD(i,j,k),i=1,imd)
 238  continue

      do 240 k=1,kmd
      do 240 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (QGD(i,j,k),i=1,imd)
 240  continue


!      do 242 k=1,kmd
!      do 242 j=1,jmd
!      irecin=irecin+1   
!      read(10,rec=irecin)  (QNRAIND(i,j,k),i=1,imd)
! 242  continue

      do 244 k=1,kmd
      do 244 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (DIABAD(i,j,k),i=1,imd)
 244  continue


      do 250 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (RAINNCVD(i,j),i=1,imd)
 250  continue



      do 255 k=1,kmd
      do 255 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (TEMPD(i,j,k),i=1,imd)
 255  continue

      do 260 k=1,kmd
      do 260 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (RHD(i,j,k),i=1,imd)
 260  continue


      do 262 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (SLPD(i,j),i=1,imd)
 262  continue

      do 270 k=1,kmd
      do 270 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (DBZD(i,j,k),i=1,imd)
 270  continue

      do 280 k=1,kmd
      do 280 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (CAPED(i,j,k),i=1,imd)
 280  continue


!  END  READ ORIGINAL DATA  



! calculation of equvalent potential temperature


      DO 710 k=1,kmd
      DO 710 j=1,jmd
      DO 710 i=1,imd
      QSAT(I,J,K)=QVD(I,J,K)*100./RHD(I,J,K)
 710  CONTINUE

! calculation of equvalent potential temperature
!'define e=0.01*(p+pb)*qvapor/(0.6219718+qvapor)+1e-20'
!'define tk2=(t+300.)*pow(((p+pb)/100000.0),(287.04/1005.7))'
!'define tlcl=55.0+2840.0/(3.5*log(tk2)-log(e)-4.805)'
!'define the=tk2*pow(100000.0/(p+pb),0.2854*(1.0-0.28*qvapor))
!      *exp(  ((3376.0/tlcl)-2.54)*qvapor*(1.0+0.81*qvapor)   )'
!  pow(x,y)(C) =x**y (fortran)


      do 440 k=1,kmd
      do 440 j=1,jmd
      do 440 i=1,imd

      E(I,J,K)=0.01*(pd(i,j,k)+pbd(i,j,k))*qvd(i,j,k)
     *        /(0.6219718+qvd(i,j,k))+1e-20

      ES(I,J,K)=0.01*(pd(i,j,k)+pbd(i,j,k))*qsat(i,j,k)
     *        /(0.6219718+qsat(i,j,k))+1e-20

     
      tk2(i,j,k)=(td(i,j,k)+300.)
     *       *((pd(i,j,k)+pbd(i,j,k))/100000.0)**(287.04/1005.7)

 440   continue

      do 450 k=1,kmd
      do 450 j=1,jmd
      do 450 i=1,imd
      tlcl(i,j,k)=55.0+2840.0/(3.5*log(tk2(i,j,k))-log(e(i,j,k))-4.805)

      tlcls(i,j,k)=55.0+2840.0/
     *          (3.5*log(tk2(i,j,k))-log(es(i,j,k))-4.805)

 450   continue

      do 460 k=1,kmd
      do 460 j=1,jmd
      do 460 i=1,imd

      te(i,j,k)=tk2(i,j,k)
     *  *(100000.0/(pd(i,j,k)+pbd(i,j,k)))**(0.2854*(1.0-
     *      0.28*qvd(i,j,k)))
     *     *exp(
     *          (   (3376.0/tlcl(i,j,k))-2.54 )
     *     *qvd(i,j,k)*(1.0+0.81*qvd(i,j,k))
     *           )


      tes(i,j,k)=tk2(i,j,k)
     *  *(100000.0/(pd(i,j,k)+pbd(i,j,k)))**(0.2854*
     *                    (1.0-0.28*qsat(i,j,k)))
     *     *exp(
     *          (   (3376.0/tlcls(i,j,k))-2.54 )
     *     *qsat(i,j,k)*(1.0+0.81*qsat(i,j,k))
     *           )

 460   continue


!**** Write temp, potential temp, and calculated te and tes

! temperature 

      do 1010 k=1,kmd
      do 1010 j=1,jmd
      irec=irec+1
 1010 write(20,rec=irec)  (TEMPD(i,j,k),i=1,imd)

! potential temperature (total)
      do 1020 k=1,kmd
      do 1020 j=1,jmd
      irec=irec+1
 1020 write(20,rec=irec)  (TD(i,j,k)+300.,i=1,imd)


! equivalent potential temperature 
      do 1030 k=1,kmd
      do 1030 j=1,jmd
      irec=irec+1
 1030 write(20,rec=irec)  (TE(i,j,k),i=1,imd)


! saturated equivalent potential temperature 
      do 1040 k=1,kmd
      do 1040 j=1,jmd
      irec=irec+1
 1040 write(20,rec=irec)  (TES(i,j,k),i=1,imd)


 2999 continue
 3000 continue

      stop
      end 
