!---------------------------------------------------------------------
!  Fortran program for Coordinate transformation (by Kawashima)
!  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 transgrid_WSM6.f 
! run:    ./a.out 
!-------------------------------------------------------------------
      PARAMETER (IMD=450,JMD=450,KMD=41,LMD=10)
      PARAMETER (IM=200,JM=200,KM=41)

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

!  COORDINATE SYSTEMS 
      REAL  XD(IMD),YD(JMD)       ! X, Y location of original grid 
      REAL  X(IM,JM),Y(IM,JM)     ! X, Y location of transformed grid 
                                  ! in the original grid 

! Center point of the transformed domain in original grid (XC,YC)
      REAL  XC(LMD),YC(LMD) 

! start point of the transformed grid in the original grid  
      REAL  XS(LMD),YS(LMD)

! ANGLE = angle that the new coordinate makes with original coordinate  
! (positive clockwise, in degree)
! ANGLE CAN BE VARIED WITH RESPECT TO TIME 

      REAL  ANGLE(LMD)   ! 

      REAL  FLGX(IM,JM,IMD)   ! FLAG for interpolation 
      REAL  FLGY(IM,JM,JMD)   ! FLAG for interpolation 
      INTEGER KD(IM,JM)   ! 
      INTEGER MD(IM,JM)   ! 


!   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 (absent in WSM6)
      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 


C  TRANSFORMED DATA (OUTPUT) 
      REAL  XLAT(IM,JM)   ! original data
      REAL  XLONG(IM,JM)   ! original data
      REAL  U(IM,JM,KM)   ! U in the old x-direction 
      REAL  V(IM,JM,KM)   ! V in the old y-direction 
      REAL  UROT(IM,JM,KM)   ! U in the new x-direction 
      REAL  VROT(IM,JM,KM)   ! V in the new y-direction 
      REAL  W(IM,JM,KM)   ! W 
      REAL  T(IM,JM,KM)   ! T
      REAL  P(IM,JM,KM)   ! P 
      REAL  PB(IM,JM,KM)   ! PB
      REAL  PH(IM,JM,KM)   ! PH 
      REAL  QV(IM,JM,KM)   ! QVAPOR
      REAL  QC(IM,JM,KM)   ! QCLOUD 
      REAL  QR(IM,JM,KM)   ! QRAIN
      REAL  QI(IM,JM,KM)   ! QICE 
      REAL  QS(IM,JM,KM)   ! QSNOW
      REAL  QG(IM,JM,KM)   ! QGRAUP 
!!      REAL  QNRAIN(IM,JM,KM)   ! QNRAIN (absent in WSM6)
      REAL  RAINNCV(IM,JM)   ! RAINNCV
      REAL  DIABA(IM,JM,KM)  ! H_DIABATIC
      REAL  TEMP(IM,JM,KM)   ! tk (temperature in K) 
      REAL  RH(IM,JM,KM)     ! rh 
      REAL  SLP(IM,JM)       ! slp 
      REAL  DBZ(IM,JM,KM)    ! dbz 
      REAL  CAPE(IM,JM,KM)   ! cape

!  Additional variables (calculated in this program) 

!***** Slopeparameter lamda, N0, dbzr
      REAL  RHO(IM,JM,KM)   ! air density
      REAL  LAMR(IM,JM,KM)  ! slope parameter 
      REAL  N0R(IM,JM,KM)   ! N0 for rain
      REAL  ZER(IM,JM,KM)   ! Ze
      REAL  DBZR(IM,JM,KM)  ! dBZ calculated from model output  

!********** Moist brunt valsala frequency ***************
      REAL  QSAT(IM,JM,KM)  ! Saturation mixing ratio 
      REAL  QSATI(IM,JM,KM)  ! Saturation mixing ratio w.r.t. ice
      REAL  BVF(IM,JM,KM)   ! Squared Brunt Vaisala frequency 
      REAL  BVFI(IM,JM,KM)  ! Squared Brunt Vaisala frequency (ice)  
      REAL  BVF2(IM,JM,KM)   ! original data
      REAL  BVFI2(IM,JM,KM)   ! original data

!****Calculation of equivalent potential temperature (following G. Bryan)
      REAL  E(IM,JM,KM)       ! vapor pressure
      REAL  TK2(IM,JM,KM)     ! 
      REAL  TLCL(IM,JM,KM)    ! temperature at the LCL 
      REAL  TE(IM,JM,KM)      ! equivalent potential temperature
      REAL  ES(IM,JM,KM)      ! saturation vapor pressure 
      REAL  ESI(IM,JM,KM)      ! saturation vapor pressure w.r.t. ice 
      REAL  TLCLS(IM,JM,KM)   ! 
      REAL  TES(IM,JM,KM)     ! saturated equivalent potential temp

!******** Others ********************
      REAL BUOY(IM,JM,KM)    ! buoyancy
      REAL GRADT(IM,JM,KM)   ! potential temperature gradient (in K/km) 
      REAL FGC(IM,JM,KM)     ! frontogenesis by confluence 
      REAL FGS(IM,JM,KM)     ! frontogenesis by shear 



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

!******************
      DXD=2.0/3.0       ! original X-grid spacing (km) 
      DYD=2.0/3.0       ! original Y-grid spacing (km)
      DZ0=250.          ! original Z-grid spacing (m)

      DX=2.0/3.0        ! grid spacing for new grid (km)
      DY=2.0/3.0        ! grid spacing for new grid (km) 
!*******************

      irecin=0
      irec=0

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


! open output data file 
      open (20,file='transgrid.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,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


!    QNRAIN is absent in WSM6 
!      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  


!**** CALCULATE COORDINATE SYSTEM**********************
! ANGLE(L): =Angle (in degree) that the new coordinate makes 
!            with the old one (positive clockwise) 
!            ANGLE(L) can be varied with respect to time 
!
!   (XC(L),YC(L)) = location of the center point (km, km) of the 
!              new coordinate in the original coordinate system  
!              (can be varied with respect to time) 
!   (XS(L),YS(L)) = location of the start point (i,j)=(1,1) 
!   of the new coordinate in the original coordinate system 
!              (can be varied with respect to time) 
!******************************************************
      pai=3.141592
      ANGLE(L)=30.0+0.0*FLOAT(L)/FLOAT(LMD)
      XC(L)=200.0
      YC(L)=200.0

      AS=sqrt(
     *        (DX*float(IM-1)/2.0)**2.
     *       +(DY*float(JM-1)/2.0)**2.
     *           )


      ALPHA=ATAND((DX*FLOAT(IM-1)/2.0)/(DY*FLOAT(JM-1)/2.0))


      XS(L)=XC(L)-AS*SIN(pai*(ANGLE(L)+ALPHA)/180.0)
      YS(L)=YC(L)-AS*COS(pai*(ANGLE(L)+ALPHA)/180.0)
       
      DO 110 I=1,IMD
      XD(I)=DXD*FLOAT(I-1)
 110   CONTINUE

      DO 120 J=1,JMD
      YD(J)=DYD*FLOAT(J-1)
 120  CONTINUE

!  calculate the location of each grid point (X(i,j),Y(i,j)) 

      DO 130 J=1,JM
      DO 130 I=1,IM
      X(I,J)=XS(L)+DY*SIN(pai*ANGLE(L)/180.0)*FLOAT(J-1)
     *         +DX*COS(pai*ANGLE(L)/180.0)*FLOAT(I-1) 

      Y(I,J)=YS(L)+DY*COS(pai*ANGLE(L)/180.0)*FLOAT(J-1)
     *         -DX*SIN(pai*ANGLE(L)/180.0)*FLOAT(I-1)
 130  CONTINUE

      do 140,J=1,JM
      write(6,*) J,X(1,J),Y(1,J)
 140  CONTINUE


      do 150 J=1,JM
      do 150 K=1,IMD-1
      do 150 I=1,IM
      FLGX(I,J,K)=(X(I,J)-XD(K))*(XD(K+1)-X(I,J))
 150  continue

      do 160 J=1,JM 
      do 160 M=1,JMD-1
      do 160 I=1,IM
      FLGY(I,J,M)=(Y(I,J)-YD(M))*(YD(M+1)-Y(I,J))
 160  continue



!**** Interpolate the original input data into new coordinate ********

      do 300 J=1,JM
      do 300 I=1,IM
      do 300 M=1,JMD-1
      do 300 K=1,IMD-1

      IF (FLGX(I,J,K).ge.0.0.AND.FLGY(I,J,M).ge.0.0) then 

      KD(I,J)= K
      MD(I,J)= M
      end if
 300  continue


      do 340 K=1,KM
      do 340 J=1,JM
      do 340 I=1,IM

      XLAT(I,J)= XLATD(KD(I,J)+1,MD(I,J)+1)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +XLATD(KD(I,J)+1,MD(I,J))
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +XLATD(KD(I,J),MD(I,J)+1)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +XLATD(KD(I,J),MD(I,J))
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      XLONG(I,J)= XLONGD(KD(I,J)+1,MD(I,J)+1)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +XLONGD(KD(I,J)+1,MD(I,J))
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +XLONGD(KD(I,J),MD(I,J)+1)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +XLONGD(KD(I,J),MD(I,J))
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD


      U(I,J,K)= UD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +UD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +UD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +UD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      V(I,J,K)= VD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +VD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +VD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +VD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      W(I,J,K)= WD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +WD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +WD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +WD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      T(I,J,K)= TD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +TD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +TD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +TD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD


      P(I,J,K)= PD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +PD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +PD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +PD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      PB(I,J,K)= PBD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +PBD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +PBD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +PBD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      PH(I,J,K)= PHD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +PHD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +PHD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +PHD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD


      QV(I,J,K)= QVD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QVD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +QVD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QVD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      QC(I,J,K)= QCD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QCD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +QCD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QCD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      QR(I,J,K)= QRD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QRD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +QRD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QRD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      QI(I,J,K)= QID(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QID(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +QID(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QID(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      QS(I,J,K)= QSD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QSD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +QSD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QSD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      QG(I,J,K)= QGD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QGD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +QGD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +QGD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

!     QNRAIN is absent in WSM6
!      QNRAIN(I,J,K)= QNRAIND(KD(I,J)+1,MD(I,J)+1,K)
!     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
!     *       +QNRAIND(KD(I,J)+1,MD(I,J),K)
!     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
!     *       +QNRAIND(KD(I,J),MD(I,J)+1,K)
!     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
!     *       +QNRAIND(KD(I,J),MD(I,J),K)
!     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      DIABA(I,J,K)= DIABAD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +DIABAD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +DIABAD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +DIABAD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD


      TEMP(I,J,K)= TEMPD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +TEMPD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +TEMPD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +TEMPD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD


      RH(I,J,K)= RHD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +RHD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +RHD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +RHD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      SLP(I,J)= SLPD(KD(I,J)+1,MD(I,J)+1)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +SLPD(KD(I,J)+1,MD(I,J))
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +SLPD(KD(I,J),MD(I,J)+1)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +SLPD(KD(I,J),MD(I,J))
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD


      DBZ(I,J,K)= DBZD(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +DBZD(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +DBZD(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +DBZD(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD

      CAPE(I,J,K)= CAPED(KD(I,J)+1,MD(I,J)+1,K)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +CAPED(KD(I,J)+1,MD(I,J),K)
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +CAPED(KD(I,J),MD(I,J)+1,K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +CAPED(KD(I,J),MD(I,J),K)
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD


      RAINNCV(I,J)= RAINNCVD(KD(I,J)+1,MD(I,J)+1)
     *        *(X(I,J)-XD(KD(I,J)))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +RAINNCVD(KD(I,J)+1,MD(I,J))
     *        *(X(I,J)-XD(KD(I,J)))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD
     *       +RAINNCVD(KD(I,J),MD(I,J)+1)
     *       *(XD(KD(I,J)+1)-X(I,J))*(Y(I,J)-YD(MD(I,J)))/DXD/DYD
     *       +RAINNCVD(KD(I,J),MD(I,J))
     *       *(XD(KD(I,J)+1)-X(I,J))*(YD(MD(I,J)+1)-Y(I,J))/DXD/DYD




 340  continue

!      To calculate system-relative wind, specify USYS and VSYS 
!      and USE UROT=(U-USYS)....
!      and USE VROT=(V-VSYS)....
!      USYS=5.0    ! system moving speed (m/s) in x-direction
!      VSYS=5.0    ! system moving speed (m/s) in y-direction

      DO 360 K=1,KM
      DO 360 J=1,JM
      DO 360 I=1,IM

!     USE below for the calculation of system-relative wind
!      UROT(I,J,K)=(U(I,J,K)-USYS)*COS(pai*ANGLE(L)/180.0)
!     *           -(V(I,J,K)-VSYS)*SIN(pai*ANGLE(L)/180.0)
!      VROT(I,J,K)=(V(I,J,K)-VSYS)*COS(pai*ANGLE(L)/180.0)
!     *           +(U(I,J,K)-USYS)*SIN(pai*ANGLE(L)/180.0)


      UROT(I,J,K)=U(I,J,K)*COS(pai*ANGLE(L)/180.0)
     *           -V(I,J,K)*SIN(pai*ANGLE(L)/180.0)

      VROT(I,J,K)=V(I,J,K)*COS(pai*ANGLE(L)/180.0)
     *           +U(I,J,K)*SIN(pai*ANGLE(L)/180.0)

      BUOY(I,J,K)=(T(I,J,K)+300)*(1+0.608*QV(I,J,K)-QC(I,J,K) 
     *            -QI(I,J,K)-QR(I,J,K)-QG(I,J,K)
     *            -QS(I,J,K))

 360   continue



!   Values at undefined points are set as undef value(=1e30)

      THV=100.
      THT=1000.
      THP=1e6
      THQ=1.0
      THH=1.0
      THQN=1e10
      THRH=120.
      THDB=120.
      THCAPE=2000.

      do 420 k=1,km
      do 420 j=1,jm
      do 420 i=1,im

      if (abs(u(i,j,k)).gt.THV) THEN
      u(i,j,k)=1e30
      end if

      if (abs(v(i,j,k)).gt.THV) THEN
      v(i,j,k)=1e30
      end if

      if (abs(urot(i,j,k)).gt.THV) THEN
      urot(i,j,k)=1e30
      end if

      if (abs(vrot(i,j,k)).gt.THV) THEN
      vrot(i,j,k)=1e30
      end if

      if (abs(w(i,j,k)).gt.THV) THEN
      w(i,j,k)=1e30
      end if

      if (abs(t(i,j,k)).gt.THT) THEN
      t(i,j,k)=1e30
      end if

      if (abs(p(i,j,k)).gt.THP) THEN
      p(i,j,k)=1e30
      end if

      if (abs(pb(i,j,k)).gt.THP) THEN
      pb(i,j,k)=1e30
      end if

      if (abs(ph(i,j,k)).gt.THP) THEN
      ph(i,j,k)=1e30
      end if

      if (abs(qv(i,j,k)).gt.THQ) THEN
      qv(i,j,k)=1e30
      end if

      if (abs(qc(i,j,k)).gt.THQ) THEN
      qc(i,j,k)=1e30
      end if

      if (abs(qr(i,j,k)).gt.THQ) THEN
      qr(i,j,k)=1e30
      end if

      if (abs(qi(i,j,k)).gt.THQ) THEN
      qi(i,j,k)=1e30
      end if

      if (abs(qs(i,j,k)).gt.THQ) THEN
      qs(i,j,k)=1e30
      end if

      if (abs(qg(i,j,k)).gt.THQ) THEN
      qg(i,j,k)=1e30
      end if

!      QNRAIN is absent in WSM6
!      if (abs(qnrain(i,j,k)).gt.THQN) THEN
!      qnrain(i,j,k)=1e30
!      end if

      if (abs(diaba(i,j,k)).gt.THH) THEN
      diaba(i,j,k)=1e30
      end if

      if (abs(temp(i,j,k)).gt.THT) THEN
      temp(i,j,k)=1e30
      end if

      if (abs(rh(i,j,k)).gt.THRH) THEN
      RH(i,j,k)=1e30
      end if


      if (abs(dbz(i,j,k)).gt.THDB) THEN
      dbz(i,j,k)=1e30
      end if

      if (abs(cape(i,j,k)).gt.THCAPE) THEN
      cape(i,j,k)=1e30
      end if
 420  continue



!**** calculate moist brunt vaisala frequency (Durran and Klemp 1982)***
 
      CL=2.5E+6
      CLI=2.834E+6
      CP=1005.0
      CR=CP*0.286
      EPS=0.622
      THQC=1.0E-9
    
      DO 710 K=1,KM
      DO 710 J=1,JM
      DO 710 I=1,IM
      QSAT(I,J,K)=QV(I,J,K)*100./RH(I,J,K)
      ESI(I,J,K)=6.11*EXP(22.514-6150.0/TEMP(I,J,K))
 710  CONTINUE

      DO 715 K=1,KM
      DO 715 J=1,JM
      DO 715 I=1,IM
      QSATI(I,J,K)=0.622*ESI(I,J,K)/
     *      ( (P(I,J,K)+PB(I,J,K)) /100. - ESI(I,J,K))

 715  CONTINUE

      DO 720 K=2,KM-1
      DO 720 J=1,JM
      DO 720 I=1,IM

      IF ((QC(I,J,K)+QI(I,J,K)).GT.THQC) THEN
      BVF(I,J,K)=9.8*( 
     *     (1.0+CL*QSAT(I,J,K)/CR/TEMP(I,J,K) )/
     *  (1.0+EPS*CL*CL*QSAT(I,J,K)/CP/CR/TEMP(I,J,K)/TEMP(I,J,K) )
     *   *  ((T(I,J,K+1)-T(I,J,K-1))/(2.0*DZ0)/
     *      (T(I,J,K)+300.)
     *      +CL/CP/TEMP(I,J,K)
     *       *(QSAT(I,J,K+1)-QSAT(I,J,K-1))/(2.0*DZ0)  )
     *    -(QC(I,J,K+1)+QI(I,J,K+1)+QR(I,J,K+1)+QG(I,J,K+1)+QS(I,J,K+1)
     *    -QC(I,J,K-1)-QI(I,J,K-1)-QR(I,J,K-1)-QG(I,J,K-1)-QS(I,J,K-1)
     *         )/(2.0*DZ0) )
      ELSE

      BVF(I,J,K)=9.8*(BUOY(I,J,K+1)-BUOY(I,J,K-1))/(2.0*DZ0)/
     *      BUOY(I,J,K)
      END IF


      IF ((QC(I,J,K)+QI(I,J,K)).GT.THQC) THEN
      BVFI(I,J,K)=9.8*( 
     *     (1.0+CLI*QSATI(I,J,K)/CR/TEMP(I,J,K) )/
     *  (1.0+EPS*CLI*CLI*QSATI(I,J,K)/CP/CR/TEMP(I,J,K)/TEMP(I,J,K) )
     *   *  ((T(I,J,K+1)-T(I,J,K-1))/(2.0*DZ0)/
     *      (T(I,J,K)+300.)
     *      +CLI/CP/TEMP(I,J,K)
     *       *(QSATI(I,J,K+1)-QSATI(I,J,K-1))/(2.0*DZ0)  )
     *    -(QC(I,J,K+1)+QI(I,J,K+1)+QR(I,J,K+1)+QG(I,J,K+1)+QS(I,J,K+1)
     *    -QC(I,J,K-1)-QI(I,J,K-1)-QR(I,J,K-1)-QG(I,J,K-1)-QS(I,J,K-1)
     *         )/(2.0*DZ0) )
      ELSE

      BVFI(I,J,K)=9.8*(BUOY(I,J,K+1)-BUOY(I,J,K-1))/(2.0*DZ0)/
     *      BUOY(I,J,K)
      END IF


 720  CONTINUE

      DO 725 J=1,JM
      DO 725 I=1,IM
      BVF(I,J,KM)=1.e30
      BVF(I,J,1)=1.e30
      BVFI(I,J,KM)=1.e30
      BVFI(I,J,1)=1.e30
 725  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,km
      do 440 j=1,jm
      do 440 i=1,im
      E(I,J,K)=0.01*(p(i,j,k)+pb(i,j,k))*qv(i,j,k)
     *        /(0.6219718+qv(i,j,k))+1e-20

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


     
      tk2(i,j,k)=(t(i,j,k)+300.)
     *       *((p(i,j,k)+pb(i,j,k))/100000.0)**(287.04/1005.7)

 440   continue

      do 450 k=1,km
      do 450 j=1,jm
      do 450 i=1,im
      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,km
      do 460 j=1,jm
      do 460 i=1,im
      te(i,j,k)=tk2(i,j,k)
     *  *(100000.0/(p(i,j,k)+pb(i,j,k)))**(0.2854*(1.0-0.28*qv(i,j,k)))
     *     *exp(
     *          (   (3376.0/tlcl(i,j,k))-2.54 )
     *     *qv(i,j,k)*(1.0+0.81*qv(i,j,k))
     *           )


      tes(i,j,k)=tk2(i,j,k)
     *  *(100000.0/(p(i,j,k)+pb(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



! Calculation of Reflectivity factor ZE for rain based on Ferrier

      pai=3.1415
      rhow=1000.0
      rhos=100.0
      rhog=400.0

      R2=1.E-8
      DO 500 K=1,KM
      DO 500 J=1,JM
      DO 500 I=1,IM
      RHO(I,J,K)=(p(i,j,k)+pb(i,j,k))/287.0/
     *          temp(i,j,k)/(1+1.608*qv(i,j,k))


      IF (QR(I,J,K).GT.R2) THEN
!     N0R=1.e8 (constant) for WSM6
      N0R(I,J,K)=1.e8 
      LAMR(I,J,K)=(pai*rhow*N0R(i,j,k)/qr(i,j,k)/rho(i,j,k))**0.25

!      LAMR(I,J,K)=(pai*rhow*qnrain(i,j,k)/qr(i,j,k))**0.33333
!      N0R(I,J,K)=qnrain(i,j,k)
!     *        *(pai*rhow*qnrain(i,j,k)/qr(i,j,k))**0.33333
      ELSE 
      LAMR(I,J,K)=1.e30
      N0R(I,J,K)=1.e30
      END IF

 500  continue


      DO 510 K=1,KM
      DO 510 J=1,JM
      DO 510 I=1,IM
      IF (QR(I,J,K).GT.R2) THEN
      ZER(I,J,K)=720.*N0R(I,J,K)/LAMR(I,J,K)**7.0*1e18
      DBZR(I,J,K)=10*log10(720.*N0R(I,J,K)/LAMR(I,J,K)**7.0*1e18)
      ELSE 
      ZER(I,J,K)=1.e30
      DBZR(I,J,K)=1.e30
      END IF

 510  continue


!**** CALCULATE temperature gradient and frontogenesis 

! GRADT= temperature gradient (K/km)

        do 320 k=1,km
        do 320 j=2,jm-1
        do 320 i=2,im-1

        GRADT(i,j,k)=SQRT(
     *   ( (T(i+1,j,k)-T(i-1,j,k)  )/2.0/DX  )**2.0
     *  +( (T(i,j+1,k)-T(i,j-1,k)  )/2.0/DY  )**2.0
     *      )


 320    CONTINUE

        do 322 k=1,km
        do 322 j=2,jm-1
        do 322 i=2,im-1
        FGC(i,j,k)=-(T(i+1,j,k)-T(i-1,j,k))/2.0/DX *
     *   ((U(i+1,j,k)-U(i-1,j,k))/2.0/DX*(T(i+1,j,k)-T(i-1,j,k))/2.0/DX
     *   +(V(i+1,j,k)-V(i-1,j,k))/2.0/DX*(T(i,j+1,k)-T(i,j-1,k))/2.0/DY)
     *   /GRADT(i,j,k)

        FGS(i,j,k)=-(T(i,j+1,k)-T(i,j-1,k))/2.0/DY *
     *   ((U(i,j+1,k)-U(i,j-1,k))/2.0/DY*(T(i+1,j,k)-T(i-1,j,k))/2.0/DX
     *   +(V(i,j+1,k)-V(i,j-1,k))/2.0/DY*(T(i,j+1,k)-T(i,j-1,k))/2.0/DY)
     *   /GRADT(i,j,k)

 322    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 k=1,km
      do 1010 j=1,jm
      irec=irec+1
 1010 write(20,rec=irec)  (U(i,j,k),i=1,im)

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


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

      do 1200 k=1,km
      do 1200 j=1,jm
      irec=irec+1
 1200 write(20,rec=irec)  (VROT(i,j,k),i=1,im)


      do 1300 k=1,km
      do 1300 j=1,jm
      irec=irec+1
 1300 write(20,rec=irec)  (W(i,j,k),i=1,im)

      do 1350 k=1,km
      do 1350 j=1,jm
      irec=irec+1
 1350 write(20,rec=irec)  (T(i,j,k),i=1,im)


      do 1400 k=1,km
      do 1400 j=1,jm
      irec=irec+1
 1400 write(20,rec=irec)  (P(i,j,k),i=1,im)

      do 1450 k=1,km
      do 1450 j=1,jm
      irec=irec+1
 1450 write(20,rec=irec)  (PB(i,j,k),i=1,im)

      do 1475 k=1,km
      do 1475 j=1,jm
      irec=irec+1
 1475 write(20,rec=irec)  (PH(i,j,k),i=1,im)


      do 1600 k=1,km
      do 1600 j=1,jm
      irec=irec+1   
      write(20,rec=irec)  (QV(i,j,k),i=1,im)
 1600  continue

      do 1650 k=1,km
      do 1650 j=1,jm
      irec=irec+1   
      write(20,rec=irec)  (QC(i,j,k),i=1,im)
 1650 continue

      do 1700 k=1,km
      do 1700 j=1,jm
      irec=irec+1   
      write(20,rec=irec)  (QR(i,j,k),i=1,im)
 1700  continue

      do 1750 k=1,km
      do 1750 j=1,jm
      irec=irec+1   
      write(20,rec=irec)  (QI(i,j,k),i=1,im)
 1750 continue

      do 1800 k=1,km
      do 1800 j=1,jm
      irec=irec+1   
      write(20,rec=irec)  (QS(i,j,k),i=1,im)
 1800 continue

      do 1850 k=1,km
      do 1850 j=1,jm
      irec=irec+1   
      write(20,rec=irec)  (QG(i,j,k),i=1,im)
 1850 continue


! QNRAIN is absent for WSM6
!      do 1900 k=1,km
!      do 1900 j=1,jm
!      irec=irec+1
! 1900 write(20,rec=irec)  (QNRAIN(i,j,k),i=1,im)


      do 1910 k=1,km
      do 1910 j=1,jm
      irec=irec+1
 1910 write(20,rec=irec)  (DIABA(i,j,k),i=1,im)


      do 1920 k=1,km
      do 1920 j=1,jm
      irec=irec+1
 1920 write(20,rec=irec)  (TEMP(i,j,k),i=1,im)

      do 1940 k=1,km
      do 1940 j=1,jm
      irec=irec+1
 1940 write(20,rec=irec)  (RH(i,j,k),i=1,im)

      do 1960 k=1,km
      do 1960 j=1,jm
      irec=irec+1
 1960 write(20,rec=irec)  (DBZ(i,j,k),i=1,im)

      do 1980 k=1,km
      do 1980 j=1,jm
      irec=irec+1
 1980 write(20,rec=irec)  (CAPE(i,j,k),i=1,im)

      do 1990 k=1,km
      do 1990 j=1,jm
      irec=irec+1
 1990 write(20,rec=irec)  (TE(i,j,k),i=1,im)

      do 2240 k=1,km
      do 2240 j=1,jm
      irec=irec+1
 2240 write(20,rec=irec)  (TES(i,j,k),i=1,im)


      do 2010 k=1,km
      do 2010 j=1,jm
      irec=irec+1
 2010 write(20,rec=irec)  (DBZR(i,j,k),i=1,im)

      do 2230 k=1,km
      do 2230 j=1,jm
      irec=irec+1
 2230 write(20,rec=irec)  (BVF(i,j,k),i=1,im)

      do 2232 k=1,km
      do 2232 j=1,jm
      irec=irec+1
 2232 write(20,rec=irec)  (BVFI(i,j,k),i=1,im)


      do 2250 k=1,km
      do 2250 j=1,jm
      irec=irec+1
 2250 write(20,rec=irec)  (RHO(i,j,k),i=1,im)

      do 2300 k=1,km
      do 2300 j=1,jm
      irec=irec+1
 2300 write(20,rec=irec)  (GRADT(i,j,k),i=1,im)


      do 2320 k=1,km
      do 2320 j=1,jm
      irec=irec+1
 2320 write(20,rec=irec)  (FGC(i,j,k),i=1,im)

      do 2340 k=1,km
      do 2340 j=1,jm
      irec=irec+1
 2340 write(20,rec=irec)  (FGS(i,j,k),i=1,im)

      do 2350 j=1,jm
      irec=irec+1
 2350 write(20,rec=irec)  (SLP(i,j),i=1,im)

      do 2360 j=1,jm
      irec=irec+1
 2360 write(20,rec=irec)  (RAINNCV(i,j),i=1,im)



 2999 continue
 3000 continue

      stop
      end 
