!! Diagnostics
!! See moddule_diagnostic for which fields are available 

MODULE module_diagnostics

  USE output_module
  USE gridinfo_module
  USE module_interp
  USE module_arrays
  USE constants_module
  USE module_pressure


  USE module_calc_pressure
  USE module_calc_height
  USE module_calc_theta
  USE module_calc_tk
  USE module_calc_tc
  USE module_calc_td
  USE module_calc_rh
  USE module_calc_uvmet
  USE module_calc_wdir
  USE module_calc_wspd
  USE module_calc_slp
  USE module_calc_dbz
  USE module_calc_cape
  USE module_calc_clfr

  CONTAINS

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: process_diagnostics
   ! Purpose: All new calls to diagnostics routines go here
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   SUBROUTINE process_diagnostics ()

      IMPLICIT NONE

      ! Local variables
      character (len=128)                    :: cname, cunits, cdesc, c_nm
      real, allocatable, dimension(:,:)      :: SCR2
      real, pointer, dimension(:,:,:)        :: SCR3
      real, allocatable, dimension(:,:,:,:)  :: SCR4
      real, pointer, dimension(:,:,:)        :: SCR
      real, pointer, dimension(:,:,:)        :: data_out 
      integer                                :: nxout, nyout, nzout
      integer                                :: ii, jj, kk



!!! If we don't have PRES (Pa) already - calculate and keep
      IF ( .not. have_PRES ) THEN 
        IF ( iprogram == 8 .AND. have_P .AND. have_PB ) THEN 
          ALLOCATE(SCR3(west_east_dim,south_north_dim,bottom_top_dim))
          SCR3 = P + PB
          c_nm = 'PRES'
          CALL keep_arrays(c_nm, SCR3)
        END IF
        IF ( iprogram == 6 .AND. have_MU .AND. have_MUB .AND. have_ZNU .AND. have_ZNW .AND. have_PTOP) THEN 
          CALL pressure(SCR3)
          c_nm = 'PRES'
          CALL keep_arrays(c_nm, SCR3)
        END IF
        DEALLOCATE(SCR3)
      END IF


!!! If we don't have TK (K) already - calculate and keep 
      IF ( .not. have_TK ) THEN 
        IF ( have_T .AND. have_PRES ) THEN 
          ALLOCATE(SCR3(west_east_dim,south_north_dim,bottom_top_dim))
          c_nm = 'TK'
!!          SCR3 = (T+300.) * ( PRES / p0 )**RCP
          SCR3 = (T+220.) * ( PRES / p0 )**RCP
          CALL keep_arrays(c_nm, SCR3)
          DEALLOCATE(SCR3)
        END IF
      END IF



!!
!! DIAGNOSTICS GO HERE
!!



!!! Calculate pressure in hPA
        IF ( INDEX(plot_these_fields,",pressure,") /= 0 ) THEN
          IF ( have_PRES ) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            CALL calc_pressure(SCR, cname, cdesc, cunits) 
            CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            DEALLOCATE(SCR)
          END IF
        END IF


!!! Calculate height from PH and PHB 
        IF ( INDEX(plot_these_fields,",height,") /= 0 ) THEN
          IF ( have_PH .AND. have_PHB ) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            CALL calc_height(SCR, cname, cdesc, cunits) 
            CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            DEALLOCATE(SCR)
          END IF
        END IF


!!! Calculate temperature 
        IF ( INDEX(plot_these_fields,",tk,") /= 0 ) THEN
          IF ( have_TK ) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            CALL calc_tk(SCR, cname, cdesc, cunits) 
            CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            DEALLOCATE(SCR)
          END IF
        END IF
        IF ( INDEX(plot_these_fields,",tc,") /= 0 ) THEN
          IF ( have_TK ) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            CALL calc_tc(SCR, cname, cdesc, cunits) 
            CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            DEALLOCATE(SCR)
          END IF
        END IF


!!! Calculate theta
        IF ( INDEX(plot_these_fields,",theta,") /= 0 ) THEN
          IF ( have_T ) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            CALL calc_theta(SCR, cname, cdesc, cunits) 
            CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            DEALLOCATE(SCR)
          END IF
        END IF


!!! Dewpoint temperature
        IF ( INDEX(plot_these_fields,",td,") /= 0 ) THEN
          IF ( have_QV .AND. have_PRES ) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            CALL calc_td(SCR, cname, cdesc, cunits) 
            CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            DEALLOCATE(SCR)
          END IF
        END IF


!!! Relative Humidity    
        IF ( INDEX(plot_these_fields,",rh,") /= 0 .OR. INDEX(plot_these_fields,",clfr,") /= 0 ) THEN
          IF ( have_TK .AND. have_QV .AND. have_PRES ) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            CALL calc_rh(SCR, cname, cdesc, cunits) 
            CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)

            !!!  Cloud fractions
            IF ( INDEX(plot_these_fields,",clfr,") /= 0 ) THEN
              IF ( ALLOCATED(SCR4) ) DEALLOCATE(SCR4)
              ALLOCATE(SCR4(west_east_dim,south_north_dim,bottom_top_dim,4))
              SCR4 = 0.0
              CALL calc_clfr(SCR4, cname, cdesc, cunits, SCR) 
              DEALLOCATE(SCR)
              ALLOCATE(SCR(west_east_dim,south_north_dim,1))
              cname = "clflo"
              cdesc = "Low Cloud Fraction"
              SCR(:,:,1) = SCR4(:,:,1,1)
              CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
              cname = "clfmi"
              cdesc = "Mid Cloud Fraction"
              SCR(:,:,1) = SCR4(:,:,1,2)
              CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
              cname = "clfhi"
              cdesc = "High Cloud Fraction"
              SCR(:,:,1) = SCR4(:,:,1,3)
              CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF

            DEALLOCATE(SCR)
          END IF
        END IF


!!! Wind speed       
        IF ( INDEX(plot_these_fields,",wspd,") /= 0 ) THEN
          IF ( have_UUU .AND. have_VVV ) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            CALL calc_wspd(SCR, cname, cdesc, cunits, 1) 
            CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            DEALLOCATE(SCR)
          END IF
        END IF


!!! Wind direction       
        IF ( INDEX(plot_these_fields,",wdir,") /= 0 ) THEN
          IF ( have_UUU .AND. have_VVV ) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            CALL calc_wdir(SCR, cname, cdesc, cunits, 1) 
            CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            DEALLOCATE(SCR)
          END IF
        END IF


!!! Wind speed  - 10m      
        IF ( INDEX(plot_these_fields,",ws10,") /= 0 ) THEN
          IF ( have_U10 .AND. have_V10 ) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,1))
            CALL calc_wspd(SCR, cname, cdesc, cunits, 0) 
            CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            DEALLOCATE(SCR)
          END IF
        END IF


!!! Wind direction  - 10m      
        IF ( INDEX(plot_these_fields,",wd10,") /= 0 ) THEN
          IF ( have_U10 .AND. have_V10 ) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,1))
            CALL calc_wdir(SCR, cname, cdesc, cunits, 0) 
            CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            DEALLOCATE(SCR)
          END IF
        END IF


!!! Rotated U and V      
        IF ( INDEX(plot_these_fields,",umet,") /= 0 .OR. &
             INDEX(plot_these_fields,",vmet,") /= 0) THEN
          IF ( have_UUU .AND. have_VVV .AND. have_XLAT .AND. have_XLONG ) THEN 

            IF ( ALLOCATED(SCR4) ) DEALLOCATE(SCR4)
            ALLOCATE(SCR4(west_east_dim,south_north_dim,bottom_top_dim,4))
            SCR4 = 0.0
            CALL calc_uvmet(SCR4, cname, cdesc, cunits) 

            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            IF ( INDEX(plot_these_fields,",umet,") /= 0 ) THEN
              cname = "umet"
              SCR(:,:,:) = SCR4(:,:,:,1)
              CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF
            IF ( INDEX(plot_these_fields,",vmet,") /= 0) THEN
              cname = "vmet"
              SCR(:,:,:) = SCR4(:,:,:,2)
              CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF
            DEALLOCATE(SCR)
            DEALLOCATE(SCR4)

          END IF
        END IF


!!! Rotated U10 and V10      
        IF ( INDEX(plot_these_fields,",u10m,") /= 0 .OR. &
             INDEX(plot_these_fields,",v10m,") /= 0) THEN
          IF ( have_U10 .AND. have_V10 .AND. have_XLAT .AND. have_XLONG ) THEN 

            IF ( ALLOCATED(SCR4) ) DEALLOCATE(SCR4)
            ALLOCATE(SCR4(west_east_dim,south_north_dim,bottom_top_dim,4))
            SCR4 = 0.0
            CALL calc_uvmet(SCR4, cname, cdesc, cunits) 

            ALLOCATE(SCR(west_east_dim,south_north_dim,1))
            IF ( INDEX(plot_these_fields,",u10m,") /= 0 ) THEN
              cname = "u10m"
              SCR(:,:,1) = SCR4(:,:,1,1)
              CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF
            IF ( INDEX(plot_these_fields,",v10m,") /= 0) THEN
              cname = "v10m"
              SCR(:,:,1) = SCR4(:,:,1,2)
              CALL interp( SCR, west_east_dim, south_north_dim, 1, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF
            DEALLOCATE(SCR)
            DEALLOCATE(SCR4)

          END IF
        END IF


!!! Sea Level Pressure  
        IF ( INDEX(plot_these_fields,",slp,") /= 0 ) THEN
          IF ( have_TK .AND. have_QV .AND. have_PRES .AND. have_PH .AND. have_PHB) THEN 
            IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)
            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            CALL calc_slp(SCR, cname, cdesc, cunits) 
            CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                         data_out, nxout, nyout, nzout, &
                         vert_array, interp_levels, number_of_zlevs) 
            CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            DEALLOCATE(SCR)
          END IF
        END IF


!!! Reflectivity             
        IF ( INDEX(plot_these_fields,",dbz,") /= 0 .OR. &
             INDEX(plot_these_fields,",max_dbz,") /= 0) THEN
          IF ( have_QV .AND. have_QR .AND. have_TK .AND. have_PRES ) THEN 
            CALL calc_dbz(SCR, cname, cdesc, cunits) 

            IF ( INDEX(plot_these_fields,",dbz,") /= 0 ) THEN
              CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF

            IF ( INDEX(plot_these_fields,",max_dbz,") /= 0) THEN
              ALLOCATE(SCR2(west_east_dim,south_north_dim))
              SCR2 = 0.0
              DO jj = 1,south_north_dim
                DO ii = 1,west_east_dim
                  DO kk = 1,bottom_top_dim
                    SCR2(ii,jj) = MAX( SCR2(ii,jj) , SCR(ii,jj,kk) )
                  END DO 
                END DO 
              END DO 
              IF (ASSOCIATED(SCR)) DEALLOCATE(SCR)
              ALLOCATE(SCR(west_east_dim,south_north_dim,1))
              SCR(:,:,1) = SCR2(:,:)
              cname = "max_dbz"
              cdesc    = "Max Reflectivity"
              CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
              DEALLOCATE(SCR2)
              DEALLOCATE(SCR)
            END IF

          END IF
        END IF



!!! 3D CAPE/CIN 
        IF ( INDEX(plot_these_fields,",cape,") /= 0 .OR. &
             INDEX(plot_these_fields,",cin,") /= 0) THEN
          IF ( have_HGT .AND. have_QV .AND. have_PRES .AND. have_TK .AND. &
               have_PH .AND. have_PHB) THEN 

            IF ( ALLOCATED(SCR4) ) DEALLOCATE(SCR4)
            ALLOCATE(SCR4(west_east_dim,south_north_dim,bottom_top_dim,4))
            SCR4 = 0.0
            CALL calc_cape(SCR4, cname, cdesc, cunits,1) 

            ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim))
            IF ( INDEX(plot_these_fields,",cape,") /= 0 ) THEN
              cname = "cape"
              cdesc = "CAPE"
              cunits = "J/kg"
              SCR(:,:,:) = SCR4(:,:,:,1)
              CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF
            IF ( INDEX(plot_these_fields,",cin,") /= 0) THEN
              cname = "cin"
              cdesc = "CIN"
              cunits = "J/kg"
              SCR(:,:,:) = SCR4(:,:,:,2)
              CALL interp (SCR, west_east_dim, south_north_dim, bottom_top_dim, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF
            DEALLOCATE(SCR)
            DEALLOCATE(SCR4)

          END IF
        END IF



!!! 2D CAPE/CIN & LCL/LFC
        IF ( INDEX(plot_these_fields,",mcape,") /= 0 .OR. &
             INDEX(plot_these_fields,",mcin,") /= 0 .OR.  &
             INDEX(plot_these_fields,",lcl,") /= 0 .OR.   &
             INDEX(plot_these_fields,",lfc,") /= 0) THEN
          IF ( have_HGT .AND. have_QV .AND. have_PRES .AND. have_TK .AND. &
               have_PH .AND. have_PHB) THEN 

            IF ( ALLOCATED(SCR4) ) DEALLOCATE(SCR4)
            ALLOCATE(SCR4(west_east_dim,south_north_dim,bottom_top_dim,4))
            SCR4 = 0.0
            CALL calc_cape(SCR4, cname, cdesc, cunits,0) 

            ALLOCATE(SCR(west_east_dim,south_north_dim,1))
            IF ( INDEX(plot_these_fields,",mcape,") /= 0 ) THEN
              cname = "mcape"
              cdesc = "MCAPE"
              cunits = "J/kg"
              SCR(:,:,1) = SCR4(:,:,1,1)
              CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF
            IF ( INDEX(plot_these_fields,",mcin,") /= 0) THEN
              cname = "mcin"
              cdesc = "MCIN"
              cunits = "J/kg"
              SCR(:,:,1) = SCR4(:,:,1,2)
              CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF
            IF ( INDEX(plot_these_fields,",lcl,") /= 0) THEN
              cname = "lcl"
              cdesc = "LCL"
              cunits = "meters AGL"
              SCR(:,:,1) = SCR4(:,:,1,3)
              CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF
            IF ( INDEX(plot_these_fields,",lfc,") /= 0) THEN
              cname = "lfc"
              cdesc = "LFC"
              cunits = "meters AGL"
              SCR(:,:,1) = SCR4(:,:,1,4)
              CALL interp (SCR, west_east_dim, south_north_dim, 1, &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            END IF
            DEALLOCATE(SCR)
            DEALLOCATE(SCR4)

          END IF
        END IF


       IF ( ALLOCATED(SCR4) ) DEALLOCATE(SCR4)
       IF ( ASSOCIATED(SCR) ) DEALLOCATE(SCR)

   END SUBROUTINE process_diagnostics

END MODULE module_diagnostics
