fv_diagnostics.F90

Go to the documentation of this file.

!***********************************************************************
!*                   GNU Lesser General Public License                 
!*
!* This file is part of the FV3 dynamical core.
!*
!* The FV3 dynamical core is free software: you can redistribute it 
!* and/or modify it under the terms of the
!* GNU Lesser General Public License as published by the
!* Free Software Foundation, either version 3 of the License, or 
!* (at your option) any later version.
!*
!* The FV3 dynamical core is distributed in the hope that it will be 
!* useful, but WITHOUT ANYWARRANTY; without even the implied warranty 
!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
!* See the GNU General Public License for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with the FV3 dynamical core.  
!* If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************

module fv_diagnostics_mod

! <table>
! <tr>
!     <th>Module Name</th>
!     <th>Functions Included</th>
!   </tr>
!   <tr>
!     <td>a2b_edge_mod</td>
!     <td>a2b_ord2, a2b_ord4</td>
!   </tr>
!   <tr>
!     <td>constants_mod</td>
!     <td>grav, rdgas, rvgas, pi=>pi_8, radius, kappa, WTMAIR, WTMCO2,
!         omega, hlv, cp_air, cp_vapor</td>
!   </tr>
!   <tr>
!     <td>diag_manager_mod</td>
!     <td>diag_axis_init, register_diag_field,
!         register_static_field, send_data, diag_grid_init</td>
!   </tr>
!   <tr>
!     <td>field_manager_mod</td>
!     <td>MODEL_ATMOS</td>
!   </tr>
!   <tr>
!     <td>fms_mod</td>
!     <td>write_version_number</td>
!   </tr>
!   <tr>
!     <td>fms_io_mod</td>
!     <td>set_domain, nullify_domain, write_version_number</td>
!   </tr>
!   <tr>
!     <td>fv_arrays_mod</td>
!     <td>fv_atmos_type, fv_grid_type, fv_diag_type, fv_grid_bounds_type, 
!         R_GRIDmax_step</td>
!   </tr>
!   <tr>
!     <td>fv_eta_mod</td>
!     <td>get_eta_level, gw_1d</td>
!   </tr>
!   <tr>
!     <td>fv_grid_utils_mod</td>
!     <td> g_sum</td>
!   </tr>
!   <tr>
!     <td>fv_io_mod</td>
!     <td>fv_io_read_tracers</td>
!   </tr>
!   <tr>
!     <td>fv_mp_mod</td>
!     <td>mp_reduce_sum, mp_reduce_min, mp_reduce_max, is_master</td>
!   </tr>
!   <tr>
!     <td>fv_mapz_mod</td>
!     <td>E_Flux, moist_cv</td>
!   </tr>
!   <tr>
!     <td>fv_sg_mod</td>
!     <td>qsmith</td>
!   </tr>
!   <tr>
!     <td>fv_surf_map_mod</td>
!     <td>zs_g</td>
!   </tr>
!   <tr>
!     <td>fv_timing_mod</td>
!     <td>timing_on, timing_off</td>
!   </tr>
!   <tr>
!     <td>gfdl_cloud_microphys_mod</td>
!     <td>wqs1, qsmith_init</td>
!   </tr>
!   <tr>
!     <td>mpp_mod</td>
!     <td>mpp_error, FATAL, stdlog, mpp_pe, mpp_root_pe, mpp_sum, mpp_max, NOTE</td>
!   </tr>
!   <tr>
!     <td>mpp_domains_mod</td>
!     <td>domain2d, mpp_update_domains, DGRID_NE</td>
!   </tr>>
!   <tr>
!     <td>sat_vapor_pres_mod</td>
!     <td>compute_qs, lookup_es</td>
!   </tr>
!   <tr>
!     <td>time_manager_mod</td>
!     <td>time_type, get_date, get_time</td>
!   </tr>
!   <tr>
!     <td>tracer_manager_mod</td>
!     <td>get_tracer_names, get_number_tracers, get_tracer_index, set_tracer_profile</td>
!   </tr>
! </table>

 use constants_mod,      only: grav, rdgas, rvgas, pi=>pi_8, radius, kappa, wtmair, wtmco2, &
                               omega, hlv, cp_air, cp_vapor
 use fms_mod,            only: write_version_number
 use fms_io_mod,         only: set_domain, nullify_domain, write_version_number
 use time_manager_mod,   only: time_type, get_date, get_time
 use mpp_domains_mod,    only: domain2d, mpp_update_domains, dgrid_ne
 use diag_manager_mod,   only: diag_axis_init, register_diag_field, &
                               register_static_field, send_data, diag_grid_init
 use fv_arrays_mod,      only: fv_atmos_type, fv_grid_type, fv_diag_type, fv_grid_bounds_type, & 
                               r_grid
 !!! CLEANUP needs rem oval?
 use fv_mapz_mod,        only: e_flux, moist_cv
 use fv_mp_mod,          only: mp_reduce_sum, mp_reduce_min, mp_reduce_max, is_master
 use fv_eta_mod,         only: get_eta_level, gw_1d
 use fv_grid_utils_mod,  only: g_sum
 use a2b_edge_mod,       only: a2b_ord2, a2b_ord4
 use fv_surf_map_mod,    only: zs_g
 use fv_sg_mod,          only: qsmith

 use tracer_manager_mod, only: get_tracer_names, get_number_tracers, get_tracer_index
 use field_manager_mod,  only: model_atmos
 use mpp_mod,            only: mpp_error, fatal, stdlog, mpp_pe, mpp_root_pe, mpp_sum, mpp_max, note
 use sat_vapor_pres_mod, only: compute_qs, lookup_es

 use fv_arrays_mod, only: max_step 
#ifndef GFS_PHYS
 use gfdl_cloud_microphys_mod, only: wqs1, qsmith_init
#endif
#ifdef MULTI_GASES
 use multi_gases_mod,  only:  virq, virqd, vicpqd, vicvqd, num_gas
#endif

 implicit none
 private


 real, parameter:: missing_value = -1.e10
 real, parameter:: missing_value2 = -1.e3 
 real, parameter:: missing_value3 = 1.e10 
 real :: ginv
 real :: pk0
 logical master
 character(len=3) :: gn = ''

! private (to this module) diag:

 type(time_type) :: fv_time
 type(fv_diag_type), pointer :: idiag

 logical :: module_is_initialized=.false.
 logical :: prt_minmax =.false.
 logical :: m_calendar
 integer  sphum, liq_wat, ice_wat, cld_amt    ! GFDL physics
 integer  rainwat, snowwat, graupel
 integer :: istep, mp_top
 real    :: ptop
 real, parameter    ::     rad2deg = 180./pi

! tracers
 character(len=128)   :: tname
 character(len=256)   :: tlongname, tunits
 real :: sphum_ll_fix = 0.
 real :: qcly0 ! initial value for terminator test

 public :: fv_diag_init, fv_time, fv_diag, prt_mxm, prt_maxmin, range_check!, id_divg, id_te
 public :: prt_mass, prt_minmax, ppme, fv_diag_init_gn, z_sum, sphum_ll_fix, eqv_pot, qcly0, gn
 public :: prt_height, prt_gb_nh_sh, interpolate_vertical, rh_calc, get_height_field, dbzcalc
 public :: max_vv,get_vorticity,max_uh 
 public :: max_vorticity,max_vorticity_hy1,bunkers_vector
 public :: helicity_relative_caps 

 integer, parameter :: nplev = 31
 integer :: levs(nplev)

! version number of this module
! Include variable "version" to be written to log file.
#include<file_version.h>

contains

 subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref)
    type(fv_atmos_type), intent(inout), target :: atm(:)
    integer, intent(out) :: axes(4)
    type(time_type), intent(in) :: time
    integer,         intent(in) :: npx, npy, npz
    real, intent(in):: p_ref

    real, allocatable :: grid_xt(:), grid_yt(:), grid_xe(:), grid_ye(:), grid_xn(:), grid_yn(:)
    real, allocatable :: grid_x(:),  grid_y(:)
    real              :: vrange(2), vsrange(2), wrange(2), trange(2), slprange(2), rhrange(2), skrange(2)
    real, allocatable :: a3(:,:,:)
    real              :: pfull(npz)
    real              :: hyam(npz), hybm(npz)

    !These id_* are not needed later since they are for static data which is not used elsewhere
    integer :: id_bk, id_pk, id_area, id_lon, id_lat, id_lont, id_latt, id_phalf, id_pfull
    integer :: id_hyam, id_hybm
    integer :: id_plev
    integer :: i, j, k, n, ntileme, id_xt, id_yt, id_x, id_y, id_xe, id_ye, id_xn, id_yn
    integer :: isc, iec, jsc, jec
    logical :: used

    character(len=64) :: plev
    character(len=64) :: field
    integer              :: ntprog
    integer              :: unit

    integer :: ncnst
    integer :: axe2(3)


    call write_version_number ( 'FV_DIAGNOSTICS_MOD', version )
    idiag => atm(1)%idiag

! For total energy diagnostics:
    idiag%steps = 0
    idiag%efx = 0.;       idiag%efx_sum = 0.
    idiag%mtq = 0.;       idiag%mtq_sum = 0.

    ncnst = atm(1)%ncnst
    m_calendar = atm(1)%flagstruct%moist_phys
    call set_domain(atm(1)%domain)  ! Set domain so that diag_manager can access tile information

    sphum   = get_tracer_index(model_atmos, 'sphum')
    liq_wat = get_tracer_index(model_atmos, 'liq_wat')
    ice_wat = get_tracer_index(model_atmos, 'ice_wat')

    rainwat = get_tracer_index(model_atmos, 'rainwat')
    snowwat = get_tracer_index(model_atmos, 'snowwat')
    graupel = get_tracer_index(model_atmos, 'graupel')
    cld_amt = get_tracer_index(model_atmos, 'cld_amt')

! valid range for some fields

!!!  This will need mods for more than 1 tile per pe  !!!

    vsrange = (/ -200.,  200. /)  ! surface (lowest layer) winds

    vrange = (/ -330.,  330. /)  ! winds
    wrange = (/ -100.,  100. /)  ! vertical wind
   rhrange = (/  -10.,  150. /)  ! RH
#ifdef HIWPP
    trange = (/    5.,  350. /)  ! temperature
#else
    trange = (/  100.,  350. /)  ! temperature
#endif
    slprange = (/800.,  1200./)  ! sea-level-pressure
    skrange  = (/ -10000000.0,  10000000.0 /)  ! dissipation estimate for SKEB

    ginv = 1./grav
     if (atm(1)%grid_number == 1) fv_time = time

    allocate ( idiag%phalf(npz+1) )
    call get_eta_level(atm(1)%npz, p_ref, pfull, idiag%phalf, atm(1)%ak, atm(1)%bk, 0.01)

    mp_top = 1
    do k=1,npz
       if ( pfull(k) > 30.e2 ) then
            mp_top = k
            exit
       endif 
    enddo
    if ( is_master() ) write(*,*) 'mp_top=', mp_top, 'pfull=', pfull(mp_top)

!   allocate(grid_xt(npx-1), grid_yt(npy-1), grid_xe(npx), grid_ye(npy-1), grid_xn(npx-1), grid_yn(npy))
    allocate(grid_xt(npx-1), grid_yt(npy-1))
    grid_xt = (/ (i, i=1,npx-1) /)
    grid_yt = (/ (j, j=1,npy-1) /)
!   grid_xe = (/ (i, i=1,npx) /)
!   grid_ye = (/ (j, j=1,npy-1) /)
!   grid_xn = (/ (i, i=1,npx-1) /)
!   grid_yn = (/ (j, j=1,npy) /)

    allocate(grid_x(npx), grid_y(npy))
    grid_x = (/ (i, i=1,npx) /)
    grid_y = (/ (j, j=1,npy) /)

    n=1
    isc = atm(n)%bd%isc; iec = atm(n)%bd%iec
    jsc = atm(n)%bd%jsc; jec = atm(n)%bd%jec

    ! Send diag_manager the grid informtaion
    call diag_grid_init(domain=atm(n)%domain, &
         &              glo_lon=rad2deg*atm(n)%gridstruct%grid(isc:iec+1,jsc:jec+1,1), &
         &              glo_lat=rad2deg*atm(n)%gridstruct%grid(isc:iec+1,jsc:jec+1,2), &
         &              aglo_lon=rad2deg*atm(n)%gridstruct%agrid(isc-1:iec+1,jsc-1:jec+1,1), &
         &              aglo_lat=rad2deg*atm(n)%gridstruct%agrid(isc-1:iec+1,jsc-1:jec+1,2))

    ntileme = size(atm(:))
    if (ntileme > 1) call mpp_error(fatal, "fv_diag_init can only be called with one grid at a time.")

!    do n = 1, ntileMe
    n = 1
       field = 'grid'

       id_xt = diag_axis_init('grid_xt',grid_xt,'degrees_E','x','T-cell longitude', &
                           set_name=trim(field),domain2=atm(n)%Domain, tile_count=n)
       id_yt = diag_axis_init('grid_yt',grid_yt,'degrees_N','y','T-cell latitude',  &
                           set_name=trim(field), domain2=atm(n)%Domain, tile_count=n)
!  Don't need these right now
!      id_xe = diag_axis_init ('grid_xe',grid_xe,'degrees_E','x','E-cell longitude', &
!                              set_name=trim(field),Domain2=Domain, tile_count=n)
!      id_ye = diag_axis_init ('grid_ye',grid_ye,'degrees_N','y','E-cell latitude',  &
!                              set_name=trim(field), Domain2=Domain, tile_count=n)
!      id_xn = diag_axis_init ('grid_xn',grid_xn,'degrees_E','x','N-cell longitude', &
!                              set_name=trim(field),Domain2=Domain, aux='geolon_n, geolat_n', tile_count=n)
!      id_yn = diag_axis_init ('grid_yn',grid_yn,'degrees_N','y','N-cell latitude',  &
!                              set_name=trim(field), Domain2=Domain, tile_count=n)

       id_x = diag_axis_init('grid_x',grid_x,'degrees_E','x','cell corner longitude', &
                           set_name=trim(field),domain2=atm(n)%Domain, tile_count=n)
       id_y = diag_axis_init('grid_y',grid_y,'degrees_N','y','cell corner latitude',  &
                           set_name=trim(field), domain2=atm(n)%Domain, tile_count=n)

!    end do
!   deallocate(grid_xt, grid_yt, grid_xe, grid_ye, grid_xn, grid_yn)
    deallocate(grid_xt, grid_yt)
    deallocate(grid_x,  grid_y )

    id_phalf = diag_axis_init('phalf', idiag%phalf, 'mb', 'z', &
            'ref half pressure level', direction=-1, set_name="dynamics")
    id_pfull = diag_axis_init('pfull', pfull, 'mb', 'z', &
            'ref full pressure level', direction=-1, set_name="dynamics", edges=id_phalf)

!---- register static fields -------

    id_bk    = register_static_field( "dynamics", 'bk', (/id_phalf/), &
         'vertical coordinate sigma value', 'none' )

    id_pk    = register_static_field( "dynamics", 'pk', (/id_phalf/), &
         'pressure part of the hybrid coordinate', 'pascal' )

    id_hyam    = register_static_field( "dynamics", 'hyam', (/id_pfull/), &
         'vertical coordinate A value', '1E-5 Pa' )

    id_hybm    = register_static_field( "dynamics", 'hybm', (/id_pfull/), &
         'vertical coordinate B value', 'none' )

!--- Send static data

    if ( id_bk > 0 )    used = send_data( id_bk,atm(1)%bk, time )
    if ( id_pk > 0 )    used = send_data( id_pk,atm(1)%ak, time )
    if ( id_hyam > 0 ) then
         do k=1,npz
            hyam(k) = 0.5 * ( atm(1)%ak(k) + atm(1)%ak(k+1) ) * 1.e-5
         enddo
         used = send_data( id_hyam, hyam, time )
    endif
    if ( id_hybm > 0 ) then
         do k=1,npz
            hybm(k) = 0.5 * ( atm(1)%bk(k) + atm(1)%bk(k+1) )
         enddo
         used = send_data( id_hybm, hybm, time )
    endif

!   Approach will need modification if we wish to write values on other than A grid.
    axes(1) = id_xt
    axes(2) = id_yt
    axes(3) = id_pfull
    axes(4) = id_phalf

! Selected pressure levels
! SJL note: 31 is enough here; if you need more levels you should do it OFF line
! do not add more to prevent the model from slowing down too much.
    levs = (/1,2,3,5,7,10,20,30,50,70,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,925,950,975,1000/)

    id_plev = diag_axis_init('plev', levs(:)*1.0, 'mb', 'z', &
            'actual pressure level', direction=-1, set_name="dynamics")

    axe2(1) = id_xt
    axe2(2) = id_yt
    axe2(3) = id_plev

!---- register time independent fields -------

!    do n = 1, ntileMe
    n = 1
       field= 'dynamics'
       id_lon  = register_static_field( trim(field), 'grid_lon', (/id_x,id_y/),  &
                                         'longitude', 'degrees_E' )
       id_lat  = register_static_field( trim(field), 'grid_lat', (/id_x,id_y/),  &
                                         'latitude', 'degrees_N' )
       id_lont = register_static_field( trim(field), 'grid_lont', (/id_xt,id_yt/),  &
                                         'longitude', 'degrees_E' )
       id_latt = register_static_field( trim(field), 'grid_latt', (/id_xt,id_yt/),  &
                                         'latitude', 'degrees_N' )
       id_area = register_static_field( trim(field), 'area', axes(1:2),  &
                                         'cell area', 'm**2' )
#ifndef DYNAMICS_ZS
       idiag%id_zsurf = register_static_field( trim(field), 'zsurf', axes(1:2),  &
                                         'surface height', 'm' )
#endif
       idiag%id_zs = register_static_field( trim(field), 'zs', axes(1:2),  &
                                        'Original Mean Terrain', 'm' )
! 3D hybrid_z fields:
       idiag%id_ze = register_static_field( trim(field), 'ze', axes(1:3),  &
                                        'Hybrid_Z_surface', 'm' )
       idiag%id_oro = register_static_field( trim(field), 'oro', axes(1:2),  &
                                        'Land/Water Mask', 'none' )
       idiag%id_sgh = register_static_field( trim(field), 'sgh', axes(1:2),  &
                                        'Terrain Standard deviation', 'm' )
!       idiag%id_ts = register_static_field ( trim(field), 'ts', axes(1:2),  &
!                                        'Skin temperature', 'K' )

!--------------------
! Initial conditions:
!--------------------
       idiag%ic_ps  = register_static_field( trim(field), 'ps_ic', axes(1:2),  &
                                         'initial surface pressure', 'Pa' )
       idiag%ic_ua = register_static_field( trim(field), 'ua_ic', axes(1:3),        &
            'zonal wind', 'm/sec' )
       idiag%ic_va = register_static_field( trim(field), 'va_ic', axes(1:3),        &
            'meridional wind', 'm/sec' )
       idiag%ic_ppt= register_static_field( trim(field), 'ppt_ic', axes(1:3),        &
            'potential temperature perturbation', 'K' )
       idiag%ic_sphum  = register_static_field( trim(field), 'sphum_ic', axes(1:2),  &
                                         'initial surface pressure', 'Pa' )

!    end do

    master = (mpp_pe()==mpp_root_pe())

    n=1
    isc = atm(n)%bd%isc; iec = atm(n)%bd%iec
    jsc = atm(n)%bd%jsc; jec = atm(n)%bd%jec

    allocate ( idiag%zsurf(isc:iec,jsc:jec) )

    do j=jsc,jec
       do i=isc,iec
          idiag%zsurf(i,j) = ginv * atm(n)%phis(i,j)
       enddo
    enddo

!--- Send time independent data

!    do n = 1, ntileMe
    n = 1
       isc = atm(n)%bd%isc; iec = atm(n)%bd%iec
       jsc = atm(n)%bd%jsc; jec = atm(n)%bd%jec
       if (id_lon  > 0) used = send_data(id_lon,  rad2deg*atm(n)%gridstruct%grid(isc:iec+1,jsc:jec+1,1), time)
       if (id_lat  > 0) used = send_data(id_lat,  rad2deg*atm(n)%gridstruct%grid(isc:iec+1,jsc:jec+1,2), time)
       if (id_lont > 0) used = send_data(id_lont, rad2deg*atm(n)%gridstruct%agrid(isc:iec,jsc:jec,1), time)
       if (id_latt > 0) used = send_data(id_latt, rad2deg*atm(n)%gridstruct%agrid(isc:iec,jsc:jec,2), time)
       if (id_area > 0) used = send_data(id_area, atm(n)%gridstruct%area(isc:iec,jsc:jec), time)
#ifndef DYNAMICS_ZS
       if (idiag%id_zsurf > 0) used = send_data(idiag%id_zsurf, idiag%zsurf, time)
#endif
       if ( atm(n)%flagstruct%fv_land ) then
         if (idiag%id_zs  > 0) used = send_data(idiag%id_zs , zs_g, time)
         if (idiag%id_oro > 0) used = send_data(idiag%id_oro, atm(n)%oro(isc:iec,jsc:jec), time)
         if (idiag%id_sgh > 0) used = send_data(idiag%id_sgh, atm(n)%sgh(isc:iec,jsc:jec), time)
       endif

       if ( atm(n)%flagstruct%ncep_ic ) then
         if (idiag%id_ts > 0) used = send_data(idiag%id_ts, atm(n)%ts(isc:iec,jsc:jec), time)
       endif

       if ( atm(n)%flagstruct%hybrid_z .and. idiag%id_ze > 0 ) &
                      used = send_data(idiag%id_ze, atm(n)%ze0(isc:iec,jsc:jec,1:npz), time)

       if (idiag%ic_ps > 0) used = send_data(idiag%ic_ps, atm(n)%ps(isc:iec,jsc:jec)*ginv, time)

       if(idiag%ic_ua > 0) used=send_data(idiag%ic_ua, atm(n)%ua(isc:iec,jsc:jec,:), time)
       if(idiag%ic_va > 0) used=send_data(idiag%ic_va, atm(n)%va(isc:iec,jsc:jec,:), time)

       pk0 = 1000.e2 ** kappa
       if(idiag%ic_ppt> 0) then
! Potential temperature
          allocate ( idiag%pt1(npz) )
          allocate ( a3(isc:iec,jsc:jec,npz) )
#ifdef TEST_GWAVES
          call gw_1d(npz, 1000.e2, atm(n)%ak, atm(n)%ak, atm(n)%ak(1), 10.e3, idiag%pt1)
#else
          idiag%pt1 = 0.
#endif
          do k=1,npz
          do j=jsc,jec
             do i=isc,iec
                a3(i,j,k) =  (atm(n)%pt(i,j,k)/atm(n)%pkz(i,j,k) - idiag%pt1(k)) * pk0
             enddo
          enddo
          enddo
          used=send_data(idiag%ic_ppt, a3, time)
          deallocate ( a3 )
          deallocate ( idiag%pt1 )
       endif
!    end do

!--------------------------------------------------------------
! Register main prognostic fields: ps, (u,v), t, omega (dp/dt)
!--------------------------------------------------------------

    allocate(idiag%id_tracer(ncnst))
    allocate(idiag%id_tracer_dmmr(ncnst))
    allocate(idiag%id_tracer_dvmr(ncnst))
    allocate(idiag%w_mr(ncnst))
    idiag%id_tracer(:)      = 0
    idiag%id_tracer_dmmr(:) = 0 
    idiag%id_tracer_dvmr(:) = 0 
    idiag%w_mr(:) = 0.e0

    allocate(idiag%id_u(nplev))
    allocate(idiag%id_v(nplev))
    allocate(idiag%id_t(nplev))
    allocate(idiag%id_h(nplev))
    allocate(idiag%id_q(nplev))
    allocate(idiag%id_omg(nplev))
    idiag%id_u(:)   = 0
    idiag%id_v(:)   = 0
    idiag%id_t(:)   = 0
    idiag%id_h(:)   = 0
    idiag%id_q(:)   = 0
    idiag%id_omg(:) = 0

!    do n = 1, ntileMe
    n = 1
       field= 'dynamics'

#ifdef DYNAMICS_ZS
       idiag%id_zsurf = register_diag_field( trim(field), 'zsurf', axes(1:2), time,           &
                                       'surface height', 'm')
#endif
!-------------------
! Surface pressure
!-------------------
       idiag%id_ps = register_diag_field( trim(field), 'ps', axes(1:2), time,           &
            'surface pressure', 'Pa', missing_value=missing_value )

!-------------------
! Mountain torque
!-------------------
       idiag%id_mq = register_diag_field( trim(field), 'mq', axes(1:2), time,           &
            'mountain torque', 'Hadleys per unit area', missing_value=missing_value )
!-------------------
! Angular momentum
!-------------------
       idiag%id_aam = register_diag_field( trim(field), 'aam', axes(1:2), time,           &
            'angular momentum', 'kg*m^2/s', missing_value=missing_value )
       idiag%id_amdt = register_diag_field( trim(field), 'amdt', axes(1:2), time,           &
            'angular momentum error', 'kg*m^2/s^2', missing_value=missing_value )

!
      do i=1,nplev
        write(plev,'(I5)') levs(i)
! Height:
        idiag%id_h(i)   = register_diag_field(trim(field), 'z'//trim(adjustl(plev)), axes(1:2), time, &
                                    trim(adjustl(plev))//'-mb height', 'm', missing_value=missing_value)
! u-wind:
        idiag%id_u(i)   = register_diag_field(trim(field), 'u'//trim(adjustl(plev)), axes(1:2), time, &
                                    trim(adjustl(plev))//'-mb u', 'm/s', missing_value=missing_value)
! v-wind:
        idiag%id_v(i)   = register_diag_field(trim(field), 'v'//trim(adjustl(plev)), axes(1:2), time, &
                                    trim(adjustl(plev))//'-mb v', 'm/s', missing_value=missing_value)
! Temperature (K):
        idiag%id_t(i)   = register_diag_field(trim(field), 't'//trim(adjustl(plev)), axes(1:2), time, &
                                    trim(adjustl(plev))//'-mb temperature', 'K', missing_value=missing_value)
! specific humidity:
        idiag%id_q(i)   = register_diag_field(trim(field), 'q'//trim(adjustl(plev)), axes(1:2), time, &
                                    trim(adjustl(plev))//'-mb specific humidity', 'kg/kg', missing_value=missing_value)
! Omega (Pa/sec)
        idiag%id_omg(i) = register_diag_field(trim(field), 'omg'//trim(adjustl(plev)), axes(1:2), time, &
                                    trim(adjustl(plev))//'-mb omega', 'Pa/s', missing_value=missing_value)
      enddo

      idiag%id_u_plev = register_diag_field( trim(field), 'u_plev', axe2(1:3), time,        &
           'zonal wind', 'm/sec', missing_value=missing_value, range=vrange )
      idiag%id_v_plev = register_diag_field( trim(field), 'v_plev', axe2(1:3), time,        &
           'meridional wind', 'm/sec', missing_value=missing_value, range=vrange )
      idiag%id_t_plev = register_diag_field( trim(field), 't_plev', axe2(1:3), time,        &
           'temperature', 'K', missing_value=missing_value, range=trange )
      idiag%id_h_plev = register_diag_field( trim(field), 'h_plev', axe2(1:3), time,        &
           'height', 'm', missing_value=missing_value )
      idiag%id_q_plev = register_diag_field( trim(field), 'q_plev', axe2(1:3), time,        &
           'specific humidity', 'kg/kg', missing_value=missing_value )
      idiag%id_omg_plev = register_diag_field( trim(field), 'omg_plev', axe2(1:3), time,        &
           'omega', 'Pa/s', missing_value=missing_value )

      ! flag for calculation of geopotential
      if ( all(idiag%id_h(minloc(abs(levs-10)))>0)  .or. all(idiag%id_h(minloc(abs(levs-50)))>0)  .or. &
           all(idiag%id_h(minloc(abs(levs-100)))>0) .or. all(idiag%id_h(minloc(abs(levs-200)))>0) .or. &
           all(idiag%id_h(minloc(abs(levs-250)))>0) .or. all(idiag%id_h(minloc(abs(levs-300)))>0) .or. &
           all(idiag%id_h(minloc(abs(levs-500)))>0) .or. all(idiag%id_h(minloc(abs(levs-700)))>0) .or. &
           all(idiag%id_h(minloc(abs(levs-850)))>0) .or. all(idiag%id_h(minloc(abs(levs-1000)))>0) ) then
           idiag%id_hght = 1
      else
           idiag%id_hght = 0
      endif
!-----------------------------
! mean temp between 300-500 mb
!-----------------------------
      idiag%id_tm = register_diag_field(trim(field), 'tm', axes(1:2),  time,   &
                                   'mean 300-500 mb temp', 'K', missing_value=missing_value )

!-------------------
! Sea-level-pressure
!-------------------
       idiag%id_slp = register_diag_field(trim(field), 'slp', axes(1:2),  time,   &
                                     'sea-level pressure', 'mb', missing_value=missing_value,  &
                                      range=slprange )
!----------------------------------
! Bottom level pressure for masking
!----------------------------------
       idiag%id_pmask = register_diag_field(trim(field), 'pmask', axes(1:2),  time,   &
                                     'masking pressure at lowest level', 'mb',   &
                                      missing_value=missing_value )
!------------------------------------------
! Fix for Bottom level pressure for masking
!------------------------------------------
       idiag%id_pmaskv2 = register_diag_field(trim(field), 'pmaskv2', axes(1:2), time,&
            & 'masking pressure at lowest level', 'mb', missing_value=missing_value)
                                     
!-------------------
! Hurricane scales:
!-------------------
! Net effects: ~ intensity * freq
       idiag%id_c15 = register_diag_field(trim(field), 'cat15', axes(1:2),  time,   &
                                     'de-pression < 1000', 'mb', missing_value=missing_value)
       idiag%id_c25 = register_diag_field(trim(field), 'cat25', axes(1:2),  time,   &
                                     'de-pression < 980', 'mb', missing_value=missing_value)
       idiag%id_c35 = register_diag_field(trim(field), 'cat35', axes(1:2),  time,   &
                                     'de-pression < 964', 'mb', missing_value=missing_value)
       idiag%id_c45 = register_diag_field(trim(field), 'cat45', axes(1:2),  time,   &
                                     'de-pression < 944', 'mb', missing_value=missing_value)
! Frequency:
       idiag%id_f15 = register_diag_field(trim(field), 'f15', axes(1:2),  time,   &
                                     'Cat15 frequency', 'none', missing_value=missing_value)
       idiag%id_f25 = register_diag_field(trim(field), 'f25', axes(1:2),  time,   &
                                     'Cat25 frequency', 'none', missing_value=missing_value)
       idiag%id_f35 = register_diag_field(trim(field), 'f35', axes(1:2),  time,   &
                                     'Cat35 frequency', 'none', missing_value=missing_value)
       idiag%id_f45 = register_diag_field(trim(field), 'f45', axes(1:2),  time,   &
                                     'Cat45 frequency', 'none', missing_value=missing_value)
!-------------------
! A grid winds (lat-lon)
!-------------------
       if (atm(n)%flagstruct%write_3d_diags) then
          idiag%id_ua = register_diag_field( trim(field), 'ucomp', axes(1:3), time,        &
               'zonal wind', 'm/sec', missing_value=missing_value, range=vrange )
          idiag%id_va = register_diag_field( trim(field), 'vcomp', axes(1:3), time,        &
               'meridional wind', 'm/sec', missing_value=missing_value, range=vrange)
          if ( .not. atm(n)%flagstruct%hydrostatic )                                        &
               idiag%id_w = register_diag_field( trim(field), 'w', axes(1:3), time,        &
               'vertical wind', 'm/sec', missing_value=missing_value, range=wrange )

          idiag%id_pt   = register_diag_field( trim(field), 'temp', axes(1:3), time,       &
               'temperature', 'K', missing_value=missing_value, range=trange )
          idiag%id_ppt  = register_diag_field( trim(field), 'ppt', axes(1:3), time,       &
               'potential temperature perturbation', 'K', missing_value=missing_value )
          idiag%id_theta_e = register_diag_field( trim(field), 'theta_e', axes(1:3), time,       &
               'theta_e', 'K', missing_value=missing_value )
          idiag%id_omga = register_diag_field( trim(field), 'omega', axes(1:3), time,      &
               'omega', 'Pa/s', missing_value=missing_value )
          idiag%id_divg  = register_diag_field( trim(field), 'divg', axes(1:3), time,      &
               'mean divergence', '1/s', missing_value=missing_value )
          ! diagnotic output for skeb testing
          idiag%id_diss = register_diag_field( trim(field), 'diss_est', axes(1:3), time,    &
               'random', 'none', missing_value=missing_value, range=skrange )

          idiag%id_rh = register_diag_field( trim(field), 'rh', axes(1:3), time,        &
               'Relative Humidity', '%', missing_value=missing_value )
          !            'Relative Humidity', '%', missing_value=missing_value, range=rhrange )
          idiag%id_delp = register_diag_field( trim(field), 'delp', axes(1:3), time,        &
               'pressure thickness', 'pa', missing_value=missing_value )
          if ( .not. atm(n)%flagstruct%hydrostatic )                                        &
               idiag%id_delz = register_diag_field( trim(field), 'delz', axes(1:3), time,        &
               'height thickness', 'm', missing_value=missing_value )
          if( atm(n)%flagstruct%hydrostatic ) then 
             idiag%id_pfhy = register_diag_field( trim(field), 'pfhy', axes(1:3), time,        &
                  'hydrostatic pressure', 'pa', missing_value=missing_value )
          else
             idiag%id_pfnh = register_diag_field( trim(field), 'pfnh', axes(1:3), time,        &
                  'non-hydrostatic pressure', 'pa', missing_value=missing_value )
          endif
          idiag%id_zratio = register_diag_field( trim(field), 'zratio', axes(1:3), time,        &
               'nonhydro_ratio', 'n/a', missing_value=missing_value )
          !--------------------
          ! 3D Condensate
          !--------------------
          idiag%id_qn = register_diag_field( trim(field), 'qn', axes(1:3), time,       &
               'cloud condensate', 'kg/m/s^2', missing_value=missing_value )
          idiag%id_qp = register_diag_field( trim(field), 'qp', axes(1:3), time,       &
               'precip condensate', 'kg/m/s^2', missing_value=missing_value )
          ! fast moist phys tendencies:
          idiag%id_mdt = register_diag_field( trim(field), 'mdt', axes(1:3), time,       &
               'DT/Dt: fast moist phys', 'deg/sec', missing_value=missing_value )
          idiag%id_qdt = register_diag_field( trim(field), 'qdt', axes(1:3), time,       &
               'Dqv/Dt: fast moist phys', 'kg/kg/sec', missing_value=missing_value )
          idiag%id_dbz = register_diag_field( trim(field), 'reflectivity', axes(1:3), time, &
               'Stoelinga simulated reflectivity', 'dBz', missing_value=missing_value)

          !--------------------
          ! Relative vorticity
          !--------------------
          idiag%id_vort = register_diag_field( trim(field), 'vort', axes(1:3), time,       &
               'vorticity', '1/s', missing_value=missing_value )
          !--------------------
          ! Potential vorticity
          !--------------------
          idiag%id_pv = register_diag_field( trim(field), 'pv', axes(1:3), time,       &
               'potential vorticity', '1/s', missing_value=missing_value )

       endif

! Total energy (only when moist_phys = .T.)
       idiag%id_te    = register_diag_field( trim(field), 'te', axes(1:2), time,      &
            'Total Energy', 'J/kg', missing_value=missing_value )
! Total Kinetic energy
       idiag%id_ke    = register_diag_field( trim(field), 'ke', axes(1:2), time,      &
            'Total KE', 'm^2/s^2', missing_value=missing_value )
       idiag%id_ws     = register_diag_field( trim(field), 'ws', axes(1:2), time,        &
            'Terrain W', 'm/s', missing_value=missing_value )
       idiag%id_maxdbz = register_diag_field( trim(field), 'max_reflectivity', axes(1:2), time, &
                'Stoelinga simulated maximum (composite) reflectivity', 'dBz', missing_value=missing_value)
       idiag%id_basedbz = register_diag_field( trim(field), 'base_reflectivity', axes(1:2), time, &
                'Stoelinga simulated base (1 km AGL) reflectivity', 'dBz', missing_value=missing_value)
       idiag%id_dbz4km = register_diag_field( trim(field), '4km_reflectivity', axes(1:2), time, &
                'Stoelinga simulated base reflectivity', 'dBz', missing_value=missing_value)
       idiag%id_dbztop = register_diag_field( trim(field), 'echo_top', axes(1:2), time, &
                'Echo top ( <= 18.5 dBz )', 'm', missing_value=missing_value2)
       idiag%id_dbz_m10C = register_diag_field( trim(field), 'm10C_reflectivity', axes(1:2), time, &
                'Reflectivity at -10C level', 'm', missing_value=missing_value)

!--------------------------
! Extra surface diagnistics:
!--------------------------
! Surface (lowest layer) vorticity: for tropical cyclones diag.
       idiag%id_vorts = register_diag_field( trim(field), 'vorts', axes(1:2), time,       &
            'surface vorticity', '1/s', missing_value=missing_value )
       idiag%id_us = register_diag_field( trim(field), 'us', axes(1:2), time,        &
            'surface u-wind', 'm/sec', missing_value=missing_value, range=vsrange )
       idiag%id_vs = register_diag_field( trim(field), 'vs', axes(1:2), time,        &
            'surface v-wind', 'm/sec', missing_value=missing_value, range=vsrange )
       idiag%id_tq = register_diag_field( trim(field), 'tq', axes(1:2), time,        &
            'Total water path', 'kg/m**2', missing_value=missing_value )
       idiag%id_iw = register_diag_field( trim(field), 'iw', axes(1:2), time,        &
            'Ice water path', 'kg/m**2', missing_value=missing_value )
       idiag%id_lw = register_diag_field( trim(field), 'lw', axes(1:2), time,        &
            'Liquid water path', 'kg/m**2', missing_value=missing_value )
       idiag%id_ts = register_diag_field( trim(field), 'ts', axes(1:2), time,  &
                                        'Skin temperature', 'K' )
       idiag%id_tb = register_diag_field( trim(field), 'tb', axes(1:2), time,  &
                                        'lowest layer temperature', 'K' )
       idiag%id_ctt = register_diag_field( trim(field), 'ctt', axes(1:2), time,  &
                                        'cloud_top temperature', 'K', missing_value=missing_value3 )
       idiag%id_ctp = register_diag_field( trim(field), 'ctp', axes(1:2), time,  &
                                        'cloud_top pressure', 'hPa' , missing_value=missing_value3 )
       idiag%id_ctz = register_diag_field( trim(field), 'ctz', axes(1:2), time,  &
                                        'cloud_top height', 'hPa' , missing_value=missing_value2 )
       idiag%id_cape = register_diag_field( trim(field), 'cape', axes(1:2), time,  &
                                        'Convective available potential energy (surface-based)', 'J/kg' , missing_value=missing_value )
       idiag%id_cin = register_diag_field( trim(field), 'cin', axes(1:2), time,  &
                                        'Convective inhibition (surface-based)', 'J/kg' , missing_value=missing_value )
#ifdef HIWPP
       idiag%id_acl = register_diag_field( trim(field), 'acl', axes(1:2), time,        &
            'Column-averaged Cl mixing ratio', 'kg/kg', missing_value=missing_value )
       idiag%id_acl2 = register_diag_field( trim(field), 'acl2', axes(1:2), time,        &
            'Column-averaged Cl2 mixing ratio', 'kg/kg', missing_value=missing_value )
       idiag%id_acly = register_diag_field( trim(field), 'acly', axes(1:2), time,        &
            'Column-averaged total chlorine mixing ratio', 'kg/kg', missing_value=missing_value )
#endif

!--------------------------
! 850-mb vorticity
!--------------------------
       idiag%id_vort850 = register_diag_field( trim(field), 'vort850', axes(1:2), time,       &
                           '850-mb vorticity', '1/s', missing_value=missing_value )

       idiag%id_vort200 = register_diag_field( trim(field), 'vort200', axes(1:2), time,       &
                           '200-mb vorticity', '1/s', missing_value=missing_value )

! Cubed_2_latlon interpolation is more accurate, particularly near the poles, using
! winds speed (a scalar), rather than wind vectors or kinetic energy directly.
       idiag%id_s200 = register_diag_field( trim(field), 's200', axes(1:2), time,       &
                           '200-mb wind_speed', 'm/s', missing_value=missing_value )
       idiag%id_sl12 = register_diag_field( trim(field), 'sl12', axes(1:2), time,       &
                           '12th L wind_speed', 'm/s', missing_value=missing_value )
       idiag%id_sl13 = register_diag_field( trim(field), 'sl13', axes(1:2), time,       &
                           '13th L wind_speed', 'm/s', missing_value=missing_value )
! Selceted (HIWPP) levels of non-precip condensates:
       idiag%id_qn200 = register_diag_field( trim(field), 'qn200', axes(1:2), time,       &
               '200mb condensate', 'kg/m/s^2', missing_value=missing_value )
       idiag%id_qn500 = register_diag_field( trim(field), 'qn500', axes(1:2), time,       &
               '500mb condensate', 'kg/m/s^2', missing_value=missing_value )
       idiag%id_qn850 = register_diag_field( trim(field), 'qn850', axes(1:2), time,       &
               '850mb condensate', 'kg/m/s^2', missing_value=missing_value )

       idiag%id_vort500 = register_diag_field( trim(field), 'vort500', axes(1:2), time,       &
                           '500-mb vorticity', '1/s', missing_value=missing_value )

       idiag%id_rain5km = register_diag_field( trim(field), 'rain5km', axes(1:2), time,       &
                           '5-km AGL liquid water', 'kg/kg', missing_value=missing_value )
!--------------------------
! w on height or pressure levels
!--------------------------
       if( .not. atm(n)%flagstruct%hydrostatic ) then
          idiag%id_w200 = register_diag_field( trim(field), 'w200', axes(1:2), time,       &
                              '200-mb w-wind', 'm/s', missing_value=missing_value )
          idiag%id_w500 = register_diag_field( trim(field), 'w500', axes(1:2), time,       &
                              '500-mb w-wind', 'm/s', missing_value=missing_value )
          idiag%id_w700 = register_diag_field( trim(field), 'w700', axes(1:2), time,       &
                              '700-mb w-wind', 'm/s', missing_value=missing_value )

          idiag%id_w850 = register_diag_field( trim(field), 'w850', axes(1:2), time,       &
                           '850-mb w-wind', 'm/s', missing_value=missing_value )
          idiag%id_w5km = register_diag_field( trim(field), 'w5km', axes(1:2), time,       &
                           '5-km AGL w-wind', 'm/s', missing_value=missing_value )
          idiag%id_w2500m = register_diag_field( trim(field), 'w2500m', axes(1:2), time,       &
                           '2.5-km AGL w-wind', 'm/s', missing_value=missing_value )
          idiag%id_w1km = register_diag_field( trim(field), 'w1km', axes(1:2), time,       &
                           '1-km AGL w-wind', 'm/s', missing_value=missing_value )
          idiag%id_wmaxup = register_diag_field( trim(field), 'wmaxup', axes(1:2), time,       &
                           'column-maximum updraft', 'm/s', missing_value=missing_value )
          idiag%id_wmaxdn = register_diag_field( trim(field), 'wmaxdn', axes(1:2), time,       &
                           'column-maximum downdraft', 'm/s', missing_value=missing_value )

       endif

! helicity
       idiag%id_x850 = register_diag_field( trim(field), 'x850', axes(1:2), time,       &
                           '850-mb vertical comp. of helicity', 'm/s**2', missing_value=missing_value )
!       idiag%id_x03 = register_diag_field ( trim(field), 'x03', axes(1:2), Time,       &
!                           '0-3 km vertical comp. of helicity', 'm**2/s**2', missing_value=missing_value )
!       idiag%id_x25 = register_diag_field ( trim(field), 'x25', axes(1:2), Time,       &
!                           '2-5 km vertical comp. of helicity', 'm**2/s**2', missing_value=missing_value )

! Storm Relative Helicity
       idiag%id_srh1 = register_diag_field( trim(field), 'srh01', axes(1:2), time,       &
                           '0-1 km Storm Relative Helicity', 'm/s**2', missing_value=missing_value )
       idiag%id_srh3 = register_diag_field( trim(field), 'srh03', axes(1:2), time,       &
                           '0-3 km Storm Relative Helicity', 'm/s**2', missing_value=missing_value )
       idiag%id_ustm = register_diag_field( trim(field), 'ustm', axes(1:2), time,       &
                           'u Component of Storm Motion', 'm/s', missing_value=missing_value )
       idiag%id_vstm = register_diag_field( trim(field), 'vstm', axes(1:2), time,       &
                           'v Component of Storm Motion', 'm/s', missing_value=missing_value )

       idiag%id_srh25 = register_diag_field( trim(field), 'srh25', axes(1:2), time,       &
                           '2-5 km Storm Relative Helicity', 'm/s**2', missing_value=missing_value )
       
       if( .not. atm(n)%flagstruct%hydrostatic ) then
          idiag%id_uh03 = register_diag_field( trim(field), 'uh03', axes(1:2), time,       &
                           '0-3 km Updraft Helicity', 'm/s**2', missing_value=missing_value )
          idiag%id_uh25 = register_diag_field( trim(field), 'uh25', axes(1:2), time,       &
                           '2-5 km Updraft Helicity', 'm/s**2', missing_value=missing_value )
       endif
! TC test winds at 100 m
       if( .not. atm(n)%flagstruct%hydrostatic )                                          &
       idiag%id_w100m = register_diag_field( trim(field), 'w100m', axes(1:2), time,       &
                        '100-m AGL w-wind', 'm/s', missing_value=missing_value )
       idiag%id_u100m = register_diag_field( trim(field), 'u100m', axes(1:2), time,       &
                        '100-m AGL u-wind', 'm/s', missing_value=missing_value )
       idiag%id_v100m = register_diag_field( trim(field), 'v100m', axes(1:2), time,       &
                        '100-m AGL v-wind', 'm/s', missing_value=missing_value )
!--------------------------
! relative humidity (physics definition):
!--------------------------
       idiag%id_rh10 = register_diag_field( trim(field), 'rh10', axes(1:2), time,       &
                           '10-mb relative humidity', '%', missing_value=missing_value )
       idiag%id_rh50 = register_diag_field( trim(field), 'rh50', axes(1:2), time,       &
                           '50-mb relative humidity', '%', missing_value=missing_value )
       idiag%id_rh100 = register_diag_field( trim(field), 'rh100', axes(1:2), time,       &
                           '100-mb relative humidity', '%', missing_value=missing_value )
       idiag%id_rh200 = register_diag_field( trim(field), 'rh200', axes(1:2), time,       &
                           '200-mb relative humidity', '%', missing_value=missing_value )
       idiag%id_rh250 = register_diag_field( trim(field), 'rh250', axes(1:2), time,       &
                           '250-mb relative humidity', '%', missing_value=missing_value )
       idiag%id_rh300 = register_diag_field( trim(field), 'rh300', axes(1:2), time,       &
                           '300-mb relative humidity', '%', missing_value=missing_value )
       idiag%id_rh500 = register_diag_field( trim(field), 'rh500', axes(1:2), time,       &
                           '500-mb relative humidity', '%', missing_value=missing_value )
       idiag%id_rh700 = register_diag_field( trim(field), 'rh700', axes(1:2), time,       &
                           '700-mb relative humidity', '%', missing_value=missing_value )
       idiag%id_rh850 = register_diag_field( trim(field), 'rh850', axes(1:2), time,       &
                           '850-mb relative humidity', '%', missing_value=missing_value )
       idiag%id_rh925 = register_diag_field( trim(field), 'rh925', axes(1:2), time,       &
                           '925-mb relative humidity', '%', missing_value=missing_value )
       idiag%id_rh1000 = register_diag_field( trim(field), 'rh1000', axes(1:2), time,       &
                           '1000-mb relative humidity', '%', missing_value=missing_value )
!--------------------------
! Dew Point
!--------------------------
       idiag%id_dp10 = register_diag_field( trim(field), 'dp10', axes(1:2), time,       &
                           '10-mb dew point', 'K', missing_value=missing_value )
       idiag%id_dp50 = register_diag_field( trim(field), 'dp50', axes(1:2), time,       &
                           '50-mb dew point', 'K', missing_value=missing_value )
       idiag%id_dp100 = register_diag_field( trim(field), 'dp100', axes(1:2), time,       &
                           '100-mb dew point', 'K', missing_value=missing_value )
       idiag%id_dp200 = register_diag_field( trim(field), 'dp200', axes(1:2), time,       &
                           '200-mb dew point', 'K', missing_value=missing_value )
       idiag%id_dp250 = register_diag_field( trim(field), 'dp250', axes(1:2), time,       &
                           '250-mb dew point', 'K', missing_value=missing_value )
       idiag%id_dp300 = register_diag_field( trim(field), 'dp300', axes(1:2), time,       &
                           '300-mb dew point', 'K', missing_value=missing_value )
       idiag%id_dp500 = register_diag_field( trim(field), 'dp500', axes(1:2), time,       &
                           '500-mb dew point', 'K', missing_value=missing_value )
       idiag%id_dp700 = register_diag_field( trim(field), 'dp700', axes(1:2), time,       &
                           '700-mb dew point', 'K', missing_value=missing_value )
       idiag%id_dp850 = register_diag_field( trim(field), 'dp850', axes(1:2), time,       &
                           '850-mb dew point', 'K', missing_value=missing_value )
       idiag%id_dp925 = register_diag_field( trim(field), 'dp925', axes(1:2), time,       &
                           '925-mb dew point', 'K', missing_value=missing_value )
       idiag%id_dp1000 = register_diag_field( trim(field), 'dp1000', axes(1:2), time,       &
                           '1000-mb dew point', 'K', missing_value=missing_value )
!--------------------------
! relative humidity (CMIP definition):
!--------------------------
       idiag%id_rh10_cmip = register_diag_field( trim(field), 'rh10_cmip', axes(1:2), time,       &
                           '10-mb relative humidity (CMIP)', '%', missing_value=missing_value )
       idiag%id_rh50_cmip = register_diag_field( trim(field), 'rh50_cmip', axes(1:2), time,       &
                           '50-mb relative humidity (CMIP)', '%', missing_value=missing_value )
       idiag%id_rh100_cmip = register_diag_field( trim(field), 'rh100_cmip', axes(1:2), time,       &
                           '100-mb relative humidity (CMIP)', '%', missing_value=missing_value )
       idiag%id_rh250_cmip = register_diag_field( trim(field), 'rh250_cmip', axes(1:2), time,       &
                           '250-mb relative humidity (CMIP)', '%', missing_value=missing_value )
       idiag%id_rh300_cmip = register_diag_field( trim(field), 'rh300_cmip', axes(1:2), time,       &
                           '300-mb relative humidity (CMIP)', '%', missing_value=missing_value )
       idiag%id_rh500_cmip = register_diag_field( trim(field), 'rh500_cmip', axes(1:2), time,       &
                           '500-mb relative humidity (CMIP)', '%', missing_value=missing_value )
       idiag%id_rh700_cmip = register_diag_field( trim(field), 'rh700_cmip', axes(1:2), time,       &
                           '700-mb relative humidity (CMIP)', '%', missing_value=missing_value )
       idiag%id_rh850_cmip = register_diag_field( trim(field), 'rh850_cmip', axes(1:2), time,       &
                           '850-mb relative humidity (CMIP)', '%', missing_value=missing_value )
       idiag%id_rh925_cmip = register_diag_field( trim(field), 'rh925_cmip', axes(1:2), time,       &
                           '925-mb relative humidity (CMIP)', '%', missing_value=missing_value )
       idiag%id_rh1000_cmip = register_diag_field( trim(field), 'rh1000_cmip', axes(1:2), time,       &
                           '1000-mb relative humidity (CMIP)', '%', missing_value=missing_value )

       if (atm(n)%flagstruct%write_3d_diags) then
          do i=1, ncnst
             !--------------------
             ! Tracer diagnostics:
             !--------------------
             call get_tracer_names ( model_atmos, i, tname, tlongname, tunits )
             idiag%id_tracer(i) = register_diag_field( field, trim(tname),  &
                  axes(1:3), time, trim(tlongname), &
                  trim(tunits), missing_value=missing_value)
             if (master) then
                if (idiag%id_tracer(i) > 0) then
                   unit = stdlog()
                   write(unit,'(a,a,a,a)') &
                        & 'Diagnostics available for tracer ',trim(tname), &
                        ' in module ', trim(field)
                end if
             endif
             !----------------------------------
             ! ESM Tracer dmmr/dvmr diagnostics:
             !   for specific elements only
             !----------------------------------
             !---co2
             if (trim(tname).eq.'co2') then
                idiag%w_mr(:) = wtmco2
                idiag%id_tracer_dmmr(i) = register_diag_field( field, trim(tname)//'_dmmr',  &
                     axes(1:3), time, trim(tlongname)//" (dry mmr)",           &
                     trim(tunits), missing_value=missing_value)
                idiag%id_tracer_dvmr(i) = register_diag_field( field, trim(tname)//'_dvmr',  &
                     axes(1:3), time, trim(tlongname)//" (dry vmr)",           &
                     'mol/mol', missing_value=missing_value)
                if (master) then
                   unit = stdlog()
                   if (idiag%id_tracer_dmmr(i) > 0) then
                      write(unit,'(a,a,a,a)') 'Diagnostics available for '//trim(tname)//' dry mmr ', &
                           trim(tname)//'_dmmr', ' in module ', trim(field)
                   end if
                   if (idiag%id_tracer_dvmr(i) > 0) then
                      write(unit,'(a,a,a,a)') 'Diagnostics available for '//trim(tname)//' dry vmr ', &
                           trim(tname)//'_dvmr', ' in module ', trim(field)
                   end if
                endif
             endif
             !---end co2

          enddo
       endif

       if ( atm(1)%flagstruct%consv_am .or. idiag%id_mq > 0 .or. idiag%id_amdt > 0 )  then
          allocate ( idiag%zxg(isc:iec,jsc:jec) )
          ! Initialize gradient of terrain for mountain torque computation:
          call init_mq(atm(n)%phis, atm(n)%gridstruct, &
               npx, npy, isc, iec, jsc, jec, atm(n)%ng)
       endif

!    end do


#ifdef TEST_TRACER
        call prt_mass(npz, atm(n)%ncnst, isc, iec, jsc, jec, atm(n)%ng, max(1,atm(n)%flagstruct%nwat),    &
                      atm(n)%ps, atm(n)%delp, atm(n)%q, atm(n)%gridstruct%area_64, atm(n)%domain)
#else
        call prt_mass(npz, atm(n)%ncnst, isc, iec, jsc, jec, atm(n)%ng, atm(n)%flagstruct%nwat,    &
                      atm(n)%ps, atm(n)%delp, atm(n)%q, atm(n)%gridstruct%area_64, atm(n)%domain)
#endif

    call nullify_domain()  ! Nullify  set_domain info

    module_is_initialized=.true.
    istep = 0
#ifndef GFS_PHYS
    if(idiag%id_theta_e >0 ) call qsmith_init
#endif
 end subroutine fv_diag_init


 subroutine init_mq(phis, gridstruct, npx, npy, is, ie, js, je, ng)
    integer, intent(in):: npx, npy, is, ie, js, je, ng
    real, intent(in):: phis(is-ng:ie+ng, js-ng:je+ng)
    type(fv_grid_type), intent(IN), target :: gridstruct

! local:
    real zs(is-ng:ie+ng, js-ng:je+ng)
    real zb(is-ng:ie+ng, js-ng:je+ng)
    real pdx(3,is:ie,js:je+1)
    real pdy(3,is:ie+1,js:je)
    integer i, j, n

    real, pointer :: rarea(:,:)
    real, pointer, dimension(:,:) :: dx, dy
    real(kind=R_GRID), pointer, dimension(:,:,:) :: en1, en2, vlon, vlat
    real, pointer, dimension(:,:,:) :: agrid

    rarea => gridstruct%rarea
    dx    => gridstruct%dx
    dy    => gridstruct%dy
    en1   => gridstruct%en1
    en2   => gridstruct%en2
    agrid => gridstruct%agrid
    vlon  => gridstruct%vlon
    vlat  => gridstruct%vlat

!   do j=js,je
!      do i=is,ie
    do j=js-ng,je+ng
       do i=is-ng,ie+ng
          zs(i,j) = phis(i,j) / grav
       enddo
    enddo
!   call mpp_update_domains( zs, domain )

!   call a2b_ord2(zs, zb, gridstruct, npx, npy, is, ie, js, je, ng)
    call a2b_ord4(zs, zb, gridstruct, npx, npy, is, ie, js, je, ng)

    do j=js,je+1
       do i=is,ie
          do n=1,3
             pdx(n,i,j) = 0.5*(zb(i,j)+zb(i+1,j))*dx(i,j)*en1(n,i,j)
          enddo
       enddo
    enddo
    do j=js,je
       do i=is,ie+1
          do n=1,3
             pdy(n,i,j) = 0.5*(zb(i,j)+zb(i,j+1))*dy(i,j)*en2(n,i,j)
          enddo
       enddo
    enddo

! Compute "volume-mean" gradient by Green's theorem
    do j=js,je
       do i=is,ie
          idiag%zxg(i,j) = vlon(i,j,1)*(pdx(1,i,j+1)-pdx(1,i,j)-pdy(1,i,j)+pdy(1,i+1,j))  &
                         + vlon(i,j,2)*(pdx(2,i,j+1)-pdx(2,i,j)-pdy(2,i,j)+pdy(2,i+1,j))  &
                         + vlon(i,j,3)*(pdx(3,i,j+1)-pdx(3,i,j)-pdy(3,i,j)+pdy(3,i+1,j))
! dF/d(lamda) = radius*cos(agrid(i,j,2)) * dF/dx, F is a scalar
!                                                       ________________________
          idiag%zxg(i,j) =  idiag%zxg(i,j)*rarea(i,j) * radius*cos(agrid(i,j,2))
!                                                       ^^^^^^^^^^^^^^^^^^^^^^^^
       enddo
    enddo

 end subroutine init_mq

 subroutine fv_diag(Atm, zvir, Time, print_freq)

    type(fv_atmos_type), intent(inout) :: atm(:)
    type(time_type),     intent(in) :: time
    real,                intent(in):: zvir
    integer,             intent(in):: print_freq

    integer :: isc, iec, jsc, jec, n, ntileme
    integer :: isd, ied, jsd, jed, npz, itrac
    integer :: ngc, nwater

    real, allocatable :: a2(:,:),a3(:,:,:), wk(:,:,:), wz(:,:,:), ucoor(:,:,:), vcoor(:,:,:)
    real, allocatable :: ustm(:,:), vstm(:,:)
    real, allocatable :: slp(:,:), depress(:,:), ws_max(:,:), tc_count(:,:)
    real, allocatable :: u2(:,:), v2(:,:), x850(:,:), var1(:,:), var2(:,:), var3(:,:)
    real, allocatable :: dmmr(:,:,:), dvmr(:,:,:)
    real height(2)
    real:: plevs(nplev), pout(nplev)
    integer:: idg(nplev), id1(nplev)
    real    :: tot_mq, tmp, sar, slon, slat
    real    :: a1d(atm(1)%npz)
!   real    :: t_gb, t_nh, t_sh, t_eq, area_gb, area_nh, area_sh, area_eq
    logical :: do_cs_intp
    logical :: used
    logical :: bad_range
    integer i,j,k, yr, mon, dd, hr, mn, days, seconds, nq, theta_d
    character(len=128)   :: tname
    real, parameter:: ws_0 = 16.   ! minimum max_wind_speed within the 7x7 search box
    real, parameter:: ws_1 = 20.
    real, parameter:: vort_c0= 2.2e-5 
    logical, allocatable :: storm(:,:), cat_crt(:,:)
    real :: tmp2, pvsum, e2, einf, qm, mm, maxdbz, allmax, rgrav
    integer :: cl, cl2
    !!! CLEANUP: does it really make sense to have this routine loop over Atm% anymore? We assume n=1 below anyway

! cat15: SLP<1000; srf_wnd>ws_0; vort>vort_c0
! cat25: SLP< 980; srf_wnd>ws_1; vort>vort_c0
! cat35: SLP< 964; srf_wnd>ws_1; vort>vort_c0
! cat45: SLP< 944; srf_wnd>ws_1; vort>vort_c0

    height(1) = 5.e3      ! for computing 5-km "pressure"
    height(2) = 0.        ! for sea-level pressure

    do i=1,nplev
       pout(i) = levs(i) * 1.e2
       plevs(i) = log( pout(i) )
    enddo

    ntileme = size(atm(:))
    n = 1
    isc = atm(n)%bd%isc; iec = atm(n)%bd%iec
    jsc = atm(n)%bd%jsc; jec = atm(n)%bd%jec
    ngc = atm(n)%ng
    npz = atm(n)%npz
    ptop = atm(n)%ak(1)
    nq = size (atm(n)%q,4)

    isd = atm(n)%bd%isd; ied = atm(n)%bd%ied
    jsd = atm(n)%bd%jsd; jed = atm(n)%bd%jed

    if( idiag%id_c15>0 ) then
        allocate (   storm(isc:iec,jsc:jec) )
        allocate ( depress(isc:iec,jsc:jec) )
        allocate (  ws_max(isc:iec,jsc:jec) )
        allocate ( cat_crt(isc:iec,jsc:jec) )
        allocate (tc_count(isc:iec,jsc:jec) )
    endif

    if( idiag%id_x850>0 ) then
        allocate ( x850(isc:iec,jsc:jec) )
    endif

    fv_time = time
    call set_domain(atm(1)%domain)

    if ( m_calendar ) then
         call get_date(fv_time, yr, mon, dd, hr, mn, seconds)
         if( print_freq == 0 ) then
                 prt_minmax = .false.
         elseif( print_freq < 0 ) then
                 istep = istep + 1
                 prt_minmax = mod(istep, -print_freq) == 0
         else
                 prt_minmax = mod(hr, print_freq) == 0 .and. mn==0 .and. seconds==0
         endif
     else
         call get_time (fv_time, seconds,  days)
         if( print_freq == 0 ) then
                 prt_minmax = .false.
         elseif( print_freq < 0 ) then
                 istep = istep + 1
                 prt_minmax = mod(istep, -print_freq) == 0
         else
                 prt_minmax = mod(seconds, 3600*print_freq) == 0
         endif
     endif

     if(prt_minmax) then
         if ( m_calendar ) then
              if(master) write(*,*) yr, mon, dd, hr, mn, seconds
         else
              if(master) write(*,*) days, seconds
         endif
     endif

    allocate ( a2(isc:iec,jsc:jec) )

    if( prt_minmax ) then

        call prt_mxm('ZS', idiag%zsurf,     isc, iec, jsc, jec, 0,   1, 1.0, atm(n)%gridstruct%area_64, atm(n)%domain)
        call prt_maxmin('PS', atm(n)%ps, isc, iec, jsc, jec, ngc, 1, 0.01)

#ifdef HIWPP
        allocate(var2(isc:iec,jsc:jec))
        !hemispheric max/min pressure
        do j=jsc,jec
        do i=isc,iec
           slat = rad2deg*atm(n)%gridstruct%agrid(i,j,2)
           if (slat >= 0.) then
              a2(i,j) = atm(n)%ps(i,j)
              var2(i,j) = 101300.
           else
              a2(i,j) = 101300.
              var2(i,j) = atm(n)%ps(i,j)
           endif
        enddo
        enddo
        call prt_maxmin('NH PS', a2, isc, iec, jsc, jec, 0, 1, 0.01)
        call prt_maxmin('SH PS', var2, isc, iec, jsc, jec, 0, 1, 0.01)

        deallocate(var2)
#endif

#ifdef TEST_TRACER
        call prt_mass(npz, nq, isc, iec, jsc, jec, ngc, max(1,atm(n)%flagstruct%nwat),    &
                      atm(n)%ps, atm(n)%delp, atm(n)%q, atm(n)%gridstruct%area_64, atm(n)%domain)
#else
        call prt_mass(npz, nq, isc, iec, jsc, jec, ngc, atm(n)%flagstruct%nwat,    &
                      atm(n)%ps, atm(n)%delp, atm(n)%q, atm(n)%gridstruct%area_64, atm(n)%domain)
#endif

#ifndef SW_DYNAMICS
        if (atm(n)%flagstruct%consv_te > 1.e-5) then
           idiag%steps = idiag%steps + 1
           idiag%efx_sum = idiag%efx_sum + e_flux
           if ( idiag%steps <= max_step ) idiag%efx(idiag%steps) = e_flux
           if (master)  then
              write(*,*) 'ENG Deficit (W/m**2)', trim(gn), '=', e_flux
           endif


        endif
        if ( .not. atm(n)%flagstruct%hydrostatic )   &
          call nh_total_energy(isc, iec, jsc, jec, isd, ied, jsd, jed, npz,  &
                               atm(n)%w, atm(n)%delz, atm(n)%pt, atm(n)%delp,  &
                               atm(n)%q, atm(n)%phis, atm(n)%gridstruct%area, atm(n)%domain, &
                               sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, atm(n)%flagstruct%nwat,     &
                               atm(n)%ua, atm(n)%va, atm(n)%flagstruct%moist_phys, a2)
#endif
        call prt_maxmin('UA_top', atm(n)%ua(isc:iec,jsc:jec,1),    &
                        isc, iec, jsc, jec, 0, 1, 1.)
        call prt_maxmin('UA', atm(n)%ua, isc, iec, jsc, jec, ngc, npz, 1.)
        call prt_maxmin('VA', atm(n)%va, isc, iec, jsc, jec, ngc, npz, 1.)

        if ( .not. atm(n)%flagstruct%hydrostatic ) then
          call prt_maxmin('W ', atm(n)%w , isc, iec, jsc, jec, ngc, npz, 1.)
          call prt_maxmin('Bottom w', atm(n)%w(isc:iec,jsc:jec,npz), isc, iec, jsc, jec, 0, 1, 1.)
          do j=jsc,jec
             do i=isc,iec
                a2(i,j) = -atm(n)%w(i,j,npz)/atm(n)%delz(i,j,npz)
             enddo
          enddo
          call prt_maxmin('Bottom: w/dz', a2, isc, iec, jsc, jec, 0, 1, 1.)

          if ( atm(n)%flagstruct%hybrid_z ) call prt_maxmin('Hybrid_ZTOP (km)', atm(n)%ze0(isc:iec,jsc:jec,1), &
                                                 isc, iec, jsc, jec, 0, 1, 1.e-3)
          call prt_maxmin('DZ (m)', atm(n)%delz(isc:iec,jsc:jec,1:npz),    &
                          isc, iec, jsc, jec, 0, npz, 1.)
          call prt_maxmin('Bottom DZ (m)', atm(n)%delz(isc:iec,jsc:jec,npz),    &
                          isc, iec, jsc, jec, 0, 1, 1.)
!         call prt_maxmin('Top DZ (m)', Atm(n)%delz(isc:iec,jsc:jec,1),    &
!                         isc, iec, jsc, jec, 0, 1, 1.)
        endif

#ifndef SW_DYNAMICS
        call prt_maxmin('TA', atm(n)%pt,   isc, iec, jsc, jec, ngc, npz, 1.)
!       call prt_maxmin('Top: TA', Atm(n)%pt(isc:iec,jsc:jec,  1), isc, iec, jsc, jec, 0, 1, 1.)
!       call prt_maxmin('Bot: TA', Atm(n)%pt(isc:iec,jsc:jec,npz), isc, iec, jsc, jec, 0, 1, 1.)
        call prt_maxmin('OM', atm(n)%omga, isc, iec, jsc, jec, ngc, npz, 1.)
#endif

    elseif ( atm(n)%flagstruct%range_warn ) then
         call range_check('DELP', atm(n)%delp, isc, iec, jsc, jec, ngc, npz, atm(n)%gridstruct%agrid,    &
                           0.01*ptop, 200.e2, bad_range)
         call range_check('UA', atm(n)%ua, isc, iec, jsc, jec, ngc, npz, atm(n)%gridstruct%agrid,   &
                           -250., 250., bad_range)
         call range_check('VA', atm(n)%va, isc, iec, jsc, jec, ngc, npz, atm(n)%gridstruct%agrid,   &
                           -250., 250., bad_range)
#ifndef SW_DYNAMICS
         call range_check('TA', atm(n)%pt, isc, iec, jsc, jec, ngc, npz, atm(n)%gridstruct%agrid,   &
#ifdef HIWPP
                           130., 350., bad_range) !DCMIP ICs have very low temperatures
#else
                           150., 350., bad_range)
#endif
#endif

      endif

    allocate ( u2(isc:iec,jsc:jec) )
    allocate ( v2(isc:iec,jsc:jec) )
    allocate ( wk(isc:iec,jsc:jec,npz) )
    if ( any(idiag%id_tracer_dmmr > 0) .or. any(idiag%id_tracer_dvmr > 0) ) then
        allocate ( dmmr(isc:iec,jsc:jec,1:npz) )
        allocate ( dvmr(isc:iec,jsc:jec,1:npz) )
    endif

!    do n = 1, ntileMe
    n = 1

#ifdef DYNAMICS_ZS
       if(idiag%id_zsurf > 0)  used=send_data(idiag%id_zsurf, idiag%zsurf, time)
#endif
       if(idiag%id_ps > 0) used=send_data(idiag%id_ps, atm(n)%ps(isc:iec,jsc:jec), time)

       if(idiag%id_c15>0 .or. idiag%id_c25>0 .or. idiag%id_c35>0 .or. idiag%id_c45>0) then
          call wind_max(isc, iec, jsc, jec ,isd, ied, jsd, jed, atm(n)%ua(isc:iec,jsc:jec,npz),   &
                        atm(n)%va(isc:iec,jsc:jec,npz), ws_max, atm(n)%domain)
          do j=jsc,jec
             do i=isc,iec
                if( abs(atm(n)%gridstruct%agrid(i,j,2)*rad2deg)<45.0 .and.     &
                    atm(n)%phis(i,j)*ginv<500.0 .and. ws_max(i,j)>ws_0 ) then
                    storm(i,j) = .true.
                else
                    storm(i,j) = .false.
                endif
             enddo
          enddo
       endif

       if ( idiag%id_vort200>0 .or. idiag%id_vort500>0 .or. idiag%id_vort850>0 .or. idiag%id_vorts>0   &
            .or. idiag%id_vort>0 .or. idiag%id_pv>0 .or. idiag%id_rh>0 .or. idiag%id_x850>0 .or.  &
            idiag%id_uh03>0 .or. idiag%id_uh25>0) then
          call get_vorticity(isc, iec, jsc, jec, isd, ied, jsd, jed, npz, atm(n)%u, atm(n)%v, wk, &
               atm(n)%gridstruct%dx, atm(n)%gridstruct%dy, atm(n)%gridstruct%rarea)

          if(idiag%id_vort >0) used=send_data(idiag%id_vort,  wk, time)
          if(idiag%id_vorts>0) used=send_data(idiag%id_vorts, wk(isc:iec,jsc:jec,npz), time)

          if(idiag%id_c15>0) then
             do j=jsc,jec
                do i=isc,iec
                   if ( storm(i,j) )    &
                   storm(i,j) = (atm(n)%gridstruct%agrid(i,j,2)>0. .and. wk(i,j,npz)> vort_c0) .or. &
                                (atm(n)%gridstruct%agrid(i,j,2)<0. .and. wk(i,j,npz)<-vort_c0) 
                enddo
             enddo
          endif

          if( idiag%id_vort200>0 ) then
             call interpolate_vertical(isc, iec, jsc, jec, npz,   &
                                       200.e2, atm(n)%peln, wk, a2)
             used=send_data(idiag%id_vort200, a2, time)
          endif
          if( idiag%id_vort500>0 ) then
             call interpolate_vertical(isc, iec, jsc, jec, npz,   &
                                       500.e2, atm(n)%peln, wk, a2)
             used=send_data(idiag%id_vort500, a2, time)
          endif

          if(idiag%id_vort850>0 .or. idiag%id_c15>0 .or. idiag%id_x850>0) then
             call interpolate_vertical(isc, iec, jsc, jec, npz,   &
                                       850.e2, atm(n)%peln, wk, a2)
             used=send_data(idiag%id_vort850, a2, time)
             if ( idiag%id_x850>0 ) x850(:,:) = a2(:,:) 

             if(idiag%id_c15>0) then
             do j=jsc,jec
                do i=isc,iec
                   if ( storm(i,j) )    &
                     storm(i,j) = (atm(n)%gridstruct%agrid(i,j,2)>0. .and. a2(i,j)> vort_c0) .or.     &
                                  (atm(n)%gridstruct%agrid(i,j,2)<0. .and. a2(i,j)<-vort_c0) 
                enddo
             enddo
             endif

          endif

          if( .not. atm(n)%flagstruct%hydrostatic ) then

             if ( idiag%id_uh03 > 0 ) then
                call updraft_helicity(isc, iec, jsc, jec, ngc, npz, zvir, sphum, a2, &
                     atm(n)%w, wk, atm(n)%delz, atm(n)%q,   &
                     atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 0., 3.e3)
                used = send_data( idiag%id_uh03, a2, time )
                if(prt_minmax) then
                   do j=jsc,jec
                      do i=isc,iec
                         tmp = rad2deg * atm(n)%gridstruct%agrid(i,j,1)
                         tmp2 = rad2deg * atm(n)%gridstruct%agrid(i,j,2)
                         if (  tmp2<25. .or. tmp2>50.    &
                              .or. tmp<235. .or. tmp>300. ) then
                            a2(i,j) = 0.
                         endif
                      enddo
                   enddo
                   call prt_maxmin('UH over CONUS', a2, isc, iec, jsc, jec, 0,   1, 1.)
                endif
             endif
             if ( idiag%id_uh25 > 0 ) then
                call updraft_helicity(isc, iec, jsc, jec, ngc, npz, zvir, sphum, a2, &
                     atm(n)%w, wk, atm(n)%delz, atm(n)%q,   &
                     atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 2.e3, 5.e3)
                used = send_data( idiag%id_uh25, a2, time )
             endif
          endif


          if ( idiag%id_srh1 > 0 .or. idiag%id_srh3 > 0 .or. idiag%id_srh25 > 0 .or. idiag%id_ustm > 0 .or. idiag%id_vstm > 0) then
              allocate(ustm(isc:iec,jsc:jec), vstm(isc:iec,jsc:jec))

              call bunkers_vector(isc, iec, jsc, jec, ngc, npz, zvir, sphum, ustm, vstm, &
                   atm(n)%ua, atm(n)%va, atm(n)%delz, atm(n)%q,   &
                   atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav)

              if ( idiag%id_ustm > 0 ) then
                 used = send_data( idiag%id_ustm, ustm, time )
              endif
              if ( idiag%id_vstm > 0 ) then
                 used = send_data( idiag%id_vstm, vstm, time )
              endif

              if ( idiag%id_srh1 > 0 ) then
                 call helicity_relative_caps(isc, iec, jsc, jec, ngc, npz, zvir, sphum, a2, ustm, vstm, &
                      atm(n)%ua, atm(n)%va, atm(n)%delz, atm(n)%q,   &
                      atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 0., 1.e3)
                 used = send_data( idiag%id_srh1, a2, time )
                 if(prt_minmax) then
                    do j=jsc,jec
                       do i=isc,iec
                          tmp = rad2deg * atm(n)%gridstruct%agrid(i,j,1)
                          tmp2 = rad2deg * atm(n)%gridstruct%agrid(i,j,2)
                          if (  tmp2<25. .or. tmp2>50.    &
                               .or. tmp<235. .or. tmp>300. ) then
                             a2(i,j) = 0.
                          endif
                       enddo
                    enddo
                    call prt_maxmin('SRH (0-1 km) over CONUS', a2, isc, iec, jsc, jec, 0,   1, 1.)
                 endif
              endif

              if ( idiag%id_srh3 > 0 ) then
                 call helicity_relative_caps(isc, iec, jsc, jec, ngc, npz, zvir, sphum, a2, ustm, vstm, &
                      atm(n)%ua, atm(n)%va, atm(n)%delz, atm(n)%q,   &
                      atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 0., 3e3)
                 used = send_data( idiag%id_srh3, a2, time )
                 if(prt_minmax) then
                    do j=jsc,jec
                       do i=isc,iec
                          tmp = rad2deg * atm(n)%gridstruct%agrid(i,j,1)
                          tmp2 = rad2deg * atm(n)%gridstruct%agrid(i,j,2)
                          if (  tmp2<25. .or. tmp2>50.    &
                               .or. tmp<235. .or. tmp>300. ) then
                             a2(i,j) = 0.
                          endif
                       enddo
                    enddo
                    call prt_maxmin('SRH (0-3 km) over CONUS', a2, isc, iec, jsc, jec, 0,   1, 1.)
                 endif
              endif

              if ( idiag%id_srh25 > 0 ) then
                 call helicity_relative_caps(isc, iec, jsc, jec, ngc, npz, zvir, sphum, a2, ustm, vstm, &
                      atm(n)%ua, atm(n)%va, atm(n)%delz, atm(n)%q,   &
                      atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 2.e3, 5e3)
                 used = send_data( idiag%id_srh25, a2, time )
                 if(prt_minmax) then
                    do j=jsc,jec
                       do i=isc,iec
                          tmp = rad2deg * atm(n)%gridstruct%agrid(i,j,1)
                          tmp2 = rad2deg * atm(n)%gridstruct%agrid(i,j,2)
                          if (  tmp2<25. .or. tmp2>50.    &
                               .or. tmp<235. .or. tmp>300. ) then
                             a2(i,j) = 0.
                          endif
                       enddo
                    enddo
                    call prt_maxmin('SRH (2-5 km) over CONUS', a2, isc, iec, jsc, jec, 0,   1, 1.)
                 endif
              endif

              deallocate(ustm, vstm)
          endif


          if ( idiag%id_pv > 0 ) then
! Note: this is expensive computation.
              call pv_entropy(isc, iec, jsc, jec, ngc, npz, wk,    &
                              atm(n)%gridstruct%f0, atm(n)%pt, atm(n)%pkz, atm(n)%delp, grav)
              used = send_data( idiag%id_pv, wk, time )
              if (prt_minmax) call prt_maxmin('PV', wk, isc, iec, jsc, jec, 0, 1, 1.)
          endif

       endif


       
!!$       if ( idiag%id_srh > 0 ) then
!!$          call helicity_relative(isc, iec, jsc, jec, ngc, npz, zvir, sphum, a2, &
!!$               Atm(n)%ua, Atm(n)%va, Atm(n)%delz, Atm(n)%q,   &
!!$               Atm(n)%flagstruct%hydrostatic, Atm(n)%pt, Atm(n)%peln, Atm(n)%phis, grav, 0., 3.e3)
!!$          used = send_data ( idiag%id_srh, a2, Time )
!!$          if(prt_minmax) then
!!$             do j=jsc,jec
!!$                do i=isc,iec
!!$                   tmp = rad2deg * Atm(n)%gridstruct%agrid(i,j,1)
!!$                   tmp2 = rad2deg * Atm(n)%gridstruct%agrid(i,j,2)
!!$                   if (  tmp2<25. .or. tmp2>50.    &
!!$                        .or. tmp<235. .or. tmp>300. ) then
!!$                      a2(i,j) = 0.
!!$                   endif
!!$                enddo
!!$             enddo
!!$             call prt_maxmin('SRH over CONUS', a2, isc, iec, jsc, jec, 0,   1, 1.)
!!$          endif
!!$       endif

!!$       if ( idiag%id_srh25 > 0 ) then
!!$          call helicity_relative(isc, iec, jsc, jec, ngc, npz, zvir, sphum, a2, &
!!$               Atm(n)%ua, Atm(n)%va, Atm(n)%delz, Atm(n)%q,   &
!!$               Atm(n)%flagstruct%hydrostatic, Atm(n)%pt, Atm(n)%peln, Atm(n)%phis, grav, 2.e3, 5.e3)
!!$          used = send_data ( idiag%id_srh25, a2, Time )
!!$       endif


       ! Relative Humidity
       if ( idiag%id_rh > 0 ) then
          ! Compute FV mean pressure
          do k=1,npz
             do j=jsc,jec
                do i=isc,iec
                   a2(i,j) = atm(n)%delp(i,j,k)/(atm(n)%peln(i,k+1,j)-atm(n)%peln(i,k,j))
                enddo
             enddo
#ifdef MULTI_GASES
             call qsmith((iec-isc+1)*(jec-jsc+1), npz,                       &
                   (iec-isc+1)*(jec-jsc+1), 1, atm(n)%pt(isc:iec,jsc:jec,k), &
                  a2, atm(n)%q(isc:iec,jsc:jec,k,sphum), wk(isc,jsc,k))
#else
             call qsmith(iec-isc+1, jec-jsc+1, 1, atm(n)%pt(isc:iec,jsc:jec,k),   &
                  a2, atm(n)%q(isc:iec,jsc:jec,k,sphum), wk(isc,jsc,k))
#endif
             do j=jsc,jec
                do i=isc,iec
                   wk(i,j,k) = 100.*atm(n)%q(i,j,k,sphum)/wk(i,j,k)
                enddo
             enddo
          enddo
          used = send_data( idiag%id_rh, wk, time )
          if(prt_minmax) then
             call prt_maxmin('RH_sf (%)', wk(isc:iec,jsc:jec,npz), isc, iec, jsc, jec, 0,   1, 1.)
             call prt_maxmin('RH_3D (%)', wk, isc, iec, jsc, jec, 0, npz, 1.)
          endif
       endif

       ! rel hum from physics at selected press levels (for IPCC)
       if (idiag%id_rh50>0  .or. idiag%id_rh100>0 .or. idiag%id_rh200>0 .or. idiag%id_rh250>0 .or. &
           idiag%id_rh300>0 .or. idiag%id_rh500>0 .or. idiag%id_rh700>0 .or. idiag%id_rh850>0 .or. &
           idiag%id_rh925>0 .or. idiag%id_rh1000>0 .or.                                            &
           idiag%id_dp50>0  .or. idiag%id_dp100>0 .or. idiag%id_dp200>0 .or. idiag%id_dp250>0 .or. &
           idiag%id_dp300>0 .or. idiag%id_dp500>0 .or. idiag%id_dp700>0 .or. idiag%id_dp850>0 .or. &
           idiag%id_dp925>0 .or. idiag%id_dp1000>0) then
           ! compute mean pressure
           do k=1,npz
               do j=jsc,jec
               do i=isc,iec
                   a2(i,j) = atm(n)%delp(i,j,k)/(atm(n)%peln(i,k+1,j)-atm(n)%peln(i,k,j))
               enddo
               enddo
               call rh_calc (a2, atm(n)%pt(isc:iec,jsc:jec,k), &
                             atm(n)%q(isc:iec,jsc:jec,k,sphum), wk(isc:iec,jsc:jec,k))
           enddo
           if (idiag%id_rh50>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 50.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh50, a2, time)
           endif
           if (idiag%id_rh100>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 100.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh100, a2, time)
           endif
           if (idiag%id_rh200>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 200.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh200, a2, time)
           endif
           if (idiag%id_rh250>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 250.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh250, a2, time)
           endif
           if (idiag%id_rh300>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 300.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh300, a2, time)
           endif
           if (idiag%id_rh500>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 500.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh500, a2, time)
           endif
           if (idiag%id_rh700>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 700.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh700, a2, time)
           endif
           if (idiag%id_rh850>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 850.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh850, a2, time)
           endif
           if (idiag%id_rh925>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 925.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh925, a2, time)
           endif
           if (idiag%id_rh1000>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 1000.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh1000, a2, time)
           endif

           if (idiag%id_dp50>0  .or. idiag%id_dp100>0 .or. idiag%id_dp200>0 .or. idiag%id_dp250>0 .or. &
                idiag%id_dp300>0 .or. idiag%id_dp500>0 .or. idiag%id_dp700>0 .or. idiag%id_dp850>0 .or. &
                idiag%id_dp925>0 .or. idiag%id_dp1000>0 ) then

              if (allocated(a3)) deallocate(a3)
              allocate(a3(isc:iec,jsc:jec,1:npz))
              !compute dew point (K)
              !using formula at https://cals.arizona.edu/azmet/dewpoint.html
              do k=1,npz
              do j=jsc,jec
              do i=isc,iec
                 tmp = ( log(max(wk(i,j,k)*1.e-2,1.e-2)) + 17.27 * ( atm(n)%pt(i,j,k) - 273.14 )/ ( -35.84 + atm(n)%pt(i,j,k)) ) / 17.27 
                 a3(i,j,k) = 273.14 + 237.3*tmp/ ( 1. - tmp )
              enddo
              enddo
              enddo

              if (idiag%id_dp50>0) then
                 call interpolate_vertical(isc, iec, jsc, jec, npz, 50.e2, atm(n)%peln, a3, a2)
                 used=send_data(idiag%id_dp50, a2, time)
              endif
              if (idiag%id_dp100>0) then
                 call interpolate_vertical(isc, iec, jsc, jec, npz, 100.e2, atm(n)%peln, a3, a2)
                 used=send_data(idiag%id_dp100, a2, time)
              endif
              if (idiag%id_dp200>0) then
                 call interpolate_vertical(isc, iec, jsc, jec, npz, 200.e2, atm(n)%peln, a3, a2)
                 used=send_data(idiag%id_dp200, a2, time)
              endif
              if (idiag%id_dp250>0) then
                 call interpolate_vertical(isc, iec, jsc, jec, npz, 250.e2, atm(n)%peln, a3, a2)
                 used=send_data(idiag%id_dp250, a2, time)
              endif
              if (idiag%id_dp300>0) then
                 call interpolate_vertical(isc, iec, jsc, jec, npz, 300.e2, atm(n)%peln, a3, a2)
                 used=send_data(idiag%id_dp300, a2, time)
              endif
              if (idiag%id_dp500>0) then
                 call interpolate_vertical(isc, iec, jsc, jec, npz, 500.e2, atm(n)%peln, a3, a2)
                 used=send_data(idiag%id_dp500, a2, time)
              endif
              if (idiag%id_dp700>0) then
                 call interpolate_vertical(isc, iec, jsc, jec, npz, 700.e2, atm(n)%peln, a3, a2)
                 used=send_data(idiag%id_dp700, a2, time)
              endif
              if (idiag%id_dp850>0) then
                 call interpolate_vertical(isc, iec, jsc, jec, npz, 850.e2, atm(n)%peln, a3, a2)
                 used=send_data(idiag%id_dp850, a2, time)
              endif
              if (idiag%id_dp925>0) then
                 call interpolate_vertical(isc, iec, jsc, jec, npz, 925.e2, atm(n)%peln, a3, a2)
                 used=send_data(idiag%id_dp925, a2, time)
              endif
              if (idiag%id_dp1000>0) then
                 call interpolate_vertical(isc, iec, jsc, jec, npz, 1000.e2, atm(n)%peln, a3, a2)
                 used=send_data(idiag%id_dp1000, a2, time)
              endif
              deallocate(a3)

           endif

       endif

       ! rel hum (CMIP definition) at selected press levels  (for IPCC)
       if (idiag%id_rh10_cmip>0  .or. idiag%id_rh50_cmip>0  .or. idiag%id_rh100_cmip>0 .or. &
           idiag%id_rh250_cmip>0 .or. idiag%id_rh300_cmip>0 .or. idiag%id_rh500_cmip>0 .or. &
           idiag%id_rh700_cmip>0 .or. idiag%id_rh850_cmip>0 .or. idiag%id_rh925_cmip>0 .or. &
           idiag%id_rh1000_cmip>0) then
           ! compute mean pressure
           do k=1,npz
               do j=jsc,jec
               do i=isc,iec
                   a2(i,j) = atm(n)%delp(i,j,k)/(atm(n)%peln(i,k+1,j)-atm(n)%peln(i,k,j))
               enddo
               enddo
               call rh_calc (a2, atm(n)%pt(isc:iec,jsc:jec,k), &
                             atm(n)%q(isc:iec,jsc:jec,k,sphum), wk(isc:iec,jsc:jec,k), do_cmip=.true.)
           enddo
           if (idiag%id_rh10_cmip>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 10.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh10_cmip, a2, time)
           endif
           if (idiag%id_rh50_cmip>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 50.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh50_cmip, a2, time)
           endif
           if (idiag%id_rh100_cmip>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 100.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh100_cmip, a2, time)
           endif
           if (idiag%id_rh250_cmip>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 250.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh250_cmip, a2, time)
           endif
           if (idiag%id_rh300_cmip>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 300.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh300_cmip, a2, time)
           endif
           if (idiag%id_rh500_cmip>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 500.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh500_cmip, a2, time)
           endif
           if (idiag%id_rh700_cmip>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 700.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh700_cmip, a2, time)
           endif
           if (idiag%id_rh850_cmip>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 850.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh850_cmip, a2, time)
           endif
           if (idiag%id_rh925_cmip>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 925.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh925_cmip, a2, time)
           endif
           if (idiag%id_rh1000_cmip>0) then
               call interpolate_vertical(isc, iec, jsc, jec, npz, 1000.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
               used=send_data(idiag%id_rh1000_cmip, a2, time)
           endif
       endif

       if(idiag%id_c25>0 .or. idiag%id_c35>0 .or. idiag%id_c45>0) then
          do j=jsc,jec
             do i=isc,iec
                if ( storm(i,j) .and. ws_max(i,j)>ws_1 ) then
                     cat_crt(i,j) = .true.
                else
                     cat_crt(i,j) = .false.
                endif
             enddo
          enddo
       endif



       if( idiag%id_slp>0 .or. idiag%id_tm>0 .or. idiag%id_hght>0 .or. idiag%id_c15>0 .or. idiag%id_ctz>0 ) then

          allocate ( wz(isc:iec,jsc:jec,npz+1) )
          call get_height_field(isc, iec, jsc, jec, ngc, npz, atm(n)%flagstruct%hydrostatic, atm(n)%delz,  &
                                wz, atm(n)%pt, atm(n)%q, atm(n)%peln, zvir)
          if( prt_minmax )   &
          call prt_mxm('ZTOP',wz(isc:iec,jsc:jec,1), isc, iec, jsc, jec, 0, 1, 1.e-3, atm(n)%gridstruct%area_64, atm(n)%domain)
!         call prt_maxmin('ZTOP', wz(isc:iec,jsc:jec,1), isc, iec, jsc, jec, 0, 1, 1.E-3)

          if(idiag%id_slp > 0) then
! Cumpute SLP (pressure at height=0)
          allocate ( slp(isc:iec,jsc:jec) )
          call get_pressure_given_height(isc, iec, jsc, jec, ngc, npz, wz, 1, height(2),   &
                                        atm(n)%pt(:,:,npz), atm(n)%peln, slp, 0.01)
          used = send_data(idiag%id_slp, slp, time)
             if( prt_minmax ) then
             call prt_maxmin('SLP', slp, isc, iec, jsc, jec, 0, 1, 1.)
! US Potential Landfall TCs (PLT):
                 do j=jsc,jec
                    do i=isc,iec
                       a2(i,j) = 1015.
                       slon = rad2deg*atm(n)%gridstruct%agrid(i,j,1)
                       slat = rad2deg*atm(n)%gridstruct%agrid(i,j,2)
                       if ( slat>15. .and. slat<40. .and. slon>270. .and. slon<290. ) then
                            a2(i,j) = slp(i,j)
                       endif
                    enddo
                 enddo
                 call prt_maxmin('ATL SLP', a2, isc, iec, jsc, jec, 0,   1, 1.)
             endif
          endif

! Compute H3000 and/or H500
          if( idiag%id_tm>0 .or. idiag%id_hght>0 .or. idiag%id_ppt>0) then

             allocate( a3(isc:iec,jsc:jec,nplev) )

             idg(:) = idiag%id_h(:)

             if ( idiag%id_tm>0 ) then
                  idg(minloc(abs(levs-300))) = 1  ! 300-mb
                  idg(minloc(abs(levs-500))) = 1  ! 500-mb
             else
                  idg(minloc(abs(levs-300))) = idiag%id_h(minloc(abs(levs-300)))
                  idg(minloc(abs(levs-500))) = idiag%id_h(minloc(abs(levs-500)))
             endif

             call get_height_given_pressure(isc, iec, jsc, jec, npz, wz, nplev, idg, plevs, atm(n)%peln, a3)
             ! reset 
             idg(minloc(abs(levs-300))) = idiag%id_h(minloc(abs(levs-300)))
             idg(minloc(abs(levs-500))) = idiag%id_h(minloc(abs(levs-500)))

             do i=1,nplev
                if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
             enddo

             if (idiag%id_h_plev>0)  then
               id1(:) = 1
               call get_height_given_pressure(isc, iec, jsc, jec, npz, wz, nplev, id1, plevs, atm(n)%peln, a3)
               used=send_data(idiag%id_h_plev, a3(isc:iec,jsc:jec,:), time)
             endif

             if( prt_minmax ) then
  
                if(all(idiag%id_h(minloc(abs(levs-100)))>0))  &
                call prt_mxm('Z100',a3(isc:iec,jsc:jec,11),isc,iec,jsc,jec,0,1,1.e-3,atm(n)%gridstruct%area_64,atm(n)%domain)

                if(all(idiag%id_h(minloc(abs(levs-500)))>0))  then
                   if (.not. atm(n)%neststruct%nested) then
#ifdef TO_BE_DELETED
                   t_eq = 0.   ;    t_nh = 0.;    t_sh = 0.;    t_gb = 0.
                   area_eq = 0.; area_nh = 0.; area_sh = 0.; area_gb = 0.
                   do j=jsc,jec
                      do i=isc,iec
                         slat = atm(n)%gridstruct%agrid(i,j,2)*rad2deg
                         area_gb = area_gb + atm(n)%gridstruct%area(i,j)
                         t_gb = t_gb + a3(i,j,19)*atm(n)%gridstruct%area(i,j)
                         if( (slat>-20. .and. slat<20.) ) then
! Tropics:
                              area_eq = area_eq + atm(n)%gridstruct%area(i,j)
                                 t_eq =    t_eq + a3(i,j,19)*atm(n)%gridstruct%area(i,j)
                         elseif( slat>=20. .and. slat<80. ) then
! NH
                              area_nh = area_nh + atm(n)%gridstruct%area(i,j)
                                 t_nh =    t_nh + a3(i,j,19)*atm(n)%gridstruct%area(i,j)
                         elseif( slat<=-20. .and. slat>-80. ) then
! SH
                              area_sh = area_sh + atm(n)%gridstruct%area(i,j)
                                 t_sh =    t_sh + a3(i,j,19)*atm(n)%gridstruct%area(i,j)
                         endif
                      enddo
                   enddo
                   call mp_reduce_sum(area_gb)
                   call mp_reduce_sum(   t_gb)
                   call mp_reduce_sum(area_nh)
                   call mp_reduce_sum(   t_nh)
                   call mp_reduce_sum(area_sh)
                   call mp_reduce_sum(   t_sh)
                   call mp_reduce_sum(area_eq)
                   call mp_reduce_sum(   t_eq)
                   if (master) write(*,*) 'Z500 GB_NH_SH_EQ=', t_gb/area_gb, t_nh/area_nh, t_sh/area_sh, t_eq/area_eq
#endif
!                  call prt_mxm('Z500',a3(isc:iec,jsc:jec,19),isc,iec,jsc,jec,0,1,1.,Atm(n)%gridstruct%area_64,Atm(n)%domain)
                   call prt_gb_nh_sh('fv_GFS Z500', isc,iec, jsc,jec, a3(isc,jsc,19), atm(n)%gridstruct%area_64(isc:iec,jsc:jec),   &
                                     atm(n)%gridstruct%agrid_64(isc:iec,jsc:jec,2))
                   endif
                endif

             endif

             ! mean virtual temp 300mb to 500mb
             if( idiag%id_tm>0 ) then
                 do j=jsc,jec
                    do i=isc,iec
                       a2(i,j) = grav*(a3(i,j,15)-a3(i,j,19))/(rdgas*(plevs(19)-plevs(15)))
                    enddo
                 enddo
                 used = send_data( idiag%id_tm, a2, time )
             endif

            if(idiag%id_c15>0 .or. idiag%id_c25>0 .or. idiag%id_c35>0 .or. idiag%id_c45>0) then
             do j=jsc,jec
                do i=isc,iec
! Minimum warm core:
                   if ( storm(i,j) ) then
                        if( a2(i,j)<254.0 .or. atm(n)%pt(i,j,npz)<281.0 ) Then
                              storm(i,j) = .false.
                            cat_crt(i,j) = .false.
                        endif
                   endif
                enddo
             enddo
! Cat 1-5:
             do j=jsc,jec
                do i=isc,iec
                   if ( storm(i,j) .and. slp(i,j)<1000.0 ) then
                         depress(i,j) = 1000. - slp(i,j)
                        tc_count(i,j) = 1.
                   else
                         depress(i,j) = 0.
                        tc_count(i,j) = 0.
                   endif
                enddo
             enddo
             used = send_data(idiag%id_c15, depress, time)
             if(idiag%id_f15>0) used = send_data(idiag%id_f15, tc_count, time)
             if(prt_minmax) then
                do j=jsc,jec
                   do i=isc,iec
                      slon = rad2deg*atm(n)%gridstruct%agrid(i,j,1)
                      slat = rad2deg*atm(n)%gridstruct%agrid(i,j,2)
! Western Pac: negative; positive elsewhere
                      if ( slat>0. .and. slat<40. .and. slon>110. .and. slon<180. ) then
                           depress(i,j) = -depress(i,j)
                      endif
                   enddo
                enddo
                call prt_maxmin('Depress', depress, isc, iec, jsc, jec, 0,   1, 1.)
                do j=jsc,jec
                   do i=isc,iec
                      if ( atm(n)%gridstruct%agrid(i,j,2)<0.) then
! Excluding the SH cyclones
                           depress(i,j) = 0.
                      endif
                   enddo
                enddo
                call prt_maxmin('NH Deps', depress, isc, iec, jsc, jec, 0,   1, 1.)

! ATL basin cyclones
                do j=jsc,jec
                   do i=isc,iec
                      tmp = rad2deg * atm(n)%gridstruct%agrid(i,j,1)
                      if ( tmp<280. ) then
                           depress(i,j) = 0.
                      endif
                   enddo
                enddo
                call prt_maxmin('ATL Deps', depress, isc, iec, jsc, jec, 0,   1, 1.)
             endif
            endif

! Cat 2-5:
            if(idiag%id_c25>0) then
             do j=jsc,jec
                do i=isc,iec
                   if ( cat_crt(i,j) .and. slp(i,j)<980.0 ) then
                        depress(i,j) = 980. - slp(i,j)
                       tc_count(i,j) = 1.
                   else
                        depress(i,j) = 0.
                       tc_count(i,j) = 0.
                   endif
                enddo
             enddo
             used = send_data(idiag%id_c25, depress, time)
             if(idiag%id_f25>0) used = send_data(idiag%id_f25, tc_count, time)
            endif

! Cat 3-5:
            if(idiag%id_c35>0) then
             do j=jsc,jec
                do i=isc,iec
                   if ( cat_crt(i,j) .and. slp(i,j)<964.0 ) then
                        depress(i,j) = 964. - slp(i,j)
                       tc_count(i,j) = 1.
                   else
                        depress(i,j) = 0.
                       tc_count(i,j) = 0.
                   endif
                enddo
             enddo
             used = send_data(idiag%id_c35, depress, time)
             if(idiag%id_f35>0) used = send_data(idiag%id_f35, tc_count, time)
            endif

! Cat 4-5:
            if(idiag%id_c45>0) then
             do j=jsc,jec
                do i=isc,iec
                   if ( cat_crt(i,j) .and. slp(i,j)<944.0 ) then
                        depress(i,j) = 944. - slp(i,j)
                       tc_count(i,j) = 1.
                   else
                        depress(i,j) = 0.
                       tc_count(i,j) = 0.
                   endif
                enddo
             enddo
             used = send_data(idiag%id_c45, depress, time)
             if(idiag%id_f45>0) used = send_data(idiag%id_f45, tc_count, time)
            endif

            if (idiag%id_c15>0) then
                deallocate(depress)
                deallocate(cat_crt)
                deallocate(storm)
                deallocate(ws_max)
                deallocate(tc_count)
            endif

            if(idiag%id_slp>0 )  deallocate( slp )

!           deallocate( a3 )
          endif

!        deallocate ( wz )
      endif

! Temperature:
       idg(:) = idiag%id_t(:)

       do_cs_intp = .false.
       do i=1,nplev
          if ( idg(i)>0 ) then
               do_cs_intp = .true.
               exit
          endif
       enddo

       if ( do_cs_intp ) then  ! log(pe) as the coordinaite for temp re-construction
          if(.not. allocated (a3) ) allocate( a3(isc:iec,jsc:jec,nplev) )
          call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%pt(isc:iec,jsc:jec,:), nplev,    &
                                plevs, wz, atm(n)%peln, idg, a3, 1)
          do i=1,nplev
             if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
          enddo
          if ( all(idiag%id_t(minloc(abs(levs-100)))>0) .and. prt_minmax ) then
             call prt_mxm('T100:', a3(isc:iec,jsc:jec,11), isc, iec, jsc, jec, 0, 1, 1.,   &
                          atm(n)%gridstruct%area_64, atm(n)%domain)
             if (.not. atm(n)%neststruct%nested)  then
                tmp = 0.
                sar = 0.
                !            Compute mean temp at 100 mb near EQ
                do j=jsc,jec
                   do i=isc,iec
                      slat = atm(n)%gridstruct%agrid(i,j,2)*rad2deg
                      if( (slat>-10.0 .and. slat<10.) ) then
                         sar = sar + atm(n)%gridstruct%area(i,j)
                         tmp = tmp + a3(i,j,11)*atm(n)%gridstruct%area(i,j)
                      endif
                   enddo
                enddo
                call mp_reduce_sum(sar)
                call mp_reduce_sum(tmp)
                if ( sar > 0. ) then
                   if (master) write(*,*) 'Tropical [10s,10n] mean T100 =', tmp/sar
                else
                   if (master) write(*,*) 'Warning: problem computing tropical mean T100'
                endif
             endif
          endif
          if ( all(idiag%id_t(minloc(abs(levs-200)))>0) .and. prt_minmax ) then
             call prt_mxm('T200:', a3(isc:iec,jsc:jec,13), isc, iec, jsc, jec, 0, 1, 1.,   &
                          atm(n)%gridstruct%area_64, atm(n)%domain)
             if (.not. atm(n)%neststruct%nested) then
                tmp = 0.
                sar = 0.
                do j=jsc,jec
                   do i=isc,iec
                      slat = atm(n)%gridstruct%agrid(i,j,2)*rad2deg
                      if( (slat>-20 .and. slat<20) ) then
                         sar = sar + atm(n)%gridstruct%area(i,j)
                         tmp = tmp + a3(i,j,13)*atm(n)%gridstruct%area(i,j)
                      endif
                   enddo
                enddo
                call mp_reduce_sum(sar)
                call mp_reduce_sum(tmp)
                if ( sar > 0. ) then
                   if (master) write(*,*) 'Tropical [-20.,20.] mean T200 =', tmp/sar
                endif
             endif
          endif
          deallocate( a3 )
       endif

       if (idiag%id_t_plev>0) then
         if(.not. allocated (a3) ) allocate( a3(isc:iec,jsc:jec,nplev) )
         id1(:) = 1
         call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%pt(isc:iec,jsc:jec,:), nplev,    &
                               plevs, wz, atm(n)%peln, id1, a3, 1)
         used=send_data(idiag%id_t_plev, a3(isc:iec,jsc:jec,:), time)
         deallocate( a3 )
       endif

       if(idiag%id_mq > 0)  then
          do j=jsc,jec
             do i=isc,iec
! zxg * surface pressure * 1.e-18--> Hadleys per unit area
! Unit Hadley = 1.E18 kg m**2 / s**2
                a2(i,j) = -1.e-18 * atm(n)%ps(i,j)*idiag%zxg(i,j)
             enddo
          enddo
          used = send_data(idiag%id_mq, a2, time)
          if( prt_minmax ) then
              tot_mq  = g_sum( atm(n)%domain, a2, isc, iec, jsc, jec, ngc, atm(n)%gridstruct%area_64, 0) 
              idiag%mtq_sum = idiag%mtq_sum + tot_mq
              if ( idiag%steps <= max_step ) idiag%mtq(idiag%steps) = tot_mq
              if(master) write(*,*) 'Total (global) mountain torque (Hadleys)=', tot_mq
          endif
       endif

       if (idiag%id_ts > 0) used = send_data(idiag%id_ts, atm(n)%ts(isc:iec,jsc:jec), time)

       if ( idiag%id_tq>0 ) then
          nwater = atm(1)%flagstruct%nwat
          a2 = 0.
          do k=1,npz
          do j=jsc,jec
             do i=isc,iec
!               a2(i,j) = a2(i,j) + Atm(n)%q(i,j,k,1)*Atm(n)%delp(i,j,k)
                a2(i,j) = a2(i,j) + sum(atm(n)%q(i,j,k,1:nwater))*atm(n)%delp(i,j,k)
             enddo
          enddo
          enddo
          used = send_data(idiag%id_tq, a2*ginv, time)
       endif
#ifdef HIWPP
       cl  = get_tracer_index(model_atmos, 'Cl')
       cl2 = get_tracer_index(model_atmos, 'Cl2')
       if (cl > 0 .and. cl2 > 0) then
        allocate(var2(isc:iec,jsc:jec))
          var2 = 0.
          do k=1,npz
          do j=jsc,jec
          do i=isc,iec
             var2(i,j) = var2(i,j) + atm(n)%delp(i,j,k)
          enddo
          enddo
          enddo

          if ( idiag%id_acl > 0 ) then
             a2 = 0.
             einf = 0.
             qm = 0.
             do k=1,npz
                do j=jsc,jec
                   do i=isc,iec
                      a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,cl)*atm(n)%delp(i,j,k) ! moist mass
                   enddo
                enddo
             enddo
             !Convert to mean mixing ratio
             do j=jsc,jec
             do i=isc,iec
                a2(i,j) = a2(i,j) / var2(i,j)
             enddo
             enddo
             used = send_data(idiag%id_acl, a2, time)
          endif
          if ( idiag%id_acl2 > 0 ) then
             a2 = 0.
             einf = 0.
             qm = 0.
             do k=1,npz
                do j=jsc,jec
                   do i=isc,iec
                      a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,cl2)*atm(n)%delp(i,j,k) ! moist mass
                   enddo
                enddo
             enddo
             !Convert to mean mixing ratio
             do j=jsc,jec
             do i=isc,iec
                a2(i,j) = a2(i,j) / var2(i,j)
             enddo
             enddo
             used = send_data(idiag%id_acl2, a2, time)
          endif
          if ( idiag%id_acly > 0 ) then
             a2 = 0.
             einf = 0.
             qm = 0.
             e2 = 0.
             do k=1,npz
                do j=jsc,jec
                   do i=isc,iec
                      mm = (atm(n)%q(i,j,k,cl)+2.*atm(n)%q(i,j,k,cl2))*atm(n)%delp(i,j,k) ! moist mass
                      a2(i,j) = a2(i,j) + mm
                      qm = qm + mm*atm(n)%gridstruct%area_64(i,j)
                   enddo
                enddo
             enddo
             !Convert to mean mixing ratio
             do j=jsc,jec
             do i=isc,iec
                a2(i,j) = a2(i,j) / var2(i,j)
             enddo
             enddo
             used = send_data(idiag%id_acly, a2, time)
             do j=jsc,jec
                do i=isc,iec
                   e2 = e2 + ((a2(i,j) - qcly0)**2)*atm(n)%gridstruct%area_64(i,j)
                   einf = max(einf, abs(a2(i,j) - qcly0))
                enddo
             enddo
             if (prt_minmax .and. .not. atm(n)%neststruct%nested) then
                call mp_reduce_sum(qm)
                call mp_reduce_max(einf)
                call mp_reduce_sum(e2)
                if (master) then
                   write(*,*) ' TERMINATOR TEST: '
                   write(*,*) '      chlorine mass: ', real(qm)/(4.*pi*radius*radius)
                   write(*,*) '             L2 err: ', sqrt(e2)/sqrt(4.*pi*radius*radius)/qcly0
                   write(*,*) '            max err: ', einf/qcly0
                endif
             endif
          endif

          deallocate(var2)

       endif
#endif
       if ( idiag%id_iw>0 ) then
          a2 = 0.
          if (ice_wat > 0) then
             do k=1,npz
             do j=jsc,jec
             do i=isc,iec
                a2(i,j) = a2(i,j) + atm(n)%delp(i,j,k) *       &
                                    atm(n)%q(i,j,k,ice_wat)
             enddo
             enddo
             enddo
          endif
          if (snowwat > 0) then
             do k=1,npz
             do j=jsc,jec
             do i=isc,iec
                a2(i,j) = a2(i,j) + atm(n)%delp(i,j,k) *      &
                                    atm(n)%q(i,j,k,snowwat)
             enddo
             enddo
             enddo
          endif
          if (graupel > 0) then
             do k=1,npz
             do j=jsc,jec
             do i=isc,iec
                a2(i,j) = a2(i,j) + atm(n)%delp(i,j,k) *      &
                                    atm(n)%q(i,j,k,graupel)
             enddo
             enddo
             enddo
          endif
          used = send_data(idiag%id_iw, a2*ginv, time)
       endif
       if ( idiag%id_lw>0 ) then
          a2 = 0.
          if (liq_wat > 0) then
             do k=1,npz
             do j=jsc,jec
             do i=isc,iec
                a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,liq_wat)*atm(n)%delp(i,j,k)
             enddo
             enddo
             enddo
          endif
          if (rainwat > 0) then
             do k=1,npz
             do j=jsc,jec
             do i=isc,iec
                a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,rainwat)*atm(n)%delp(i,j,k)
             enddo
             enddo
             enddo
          endif
          used = send_data(idiag%id_lw, a2*ginv, time)
       endif

! Cloud top temperature & cloud top press:
       if ( (idiag%id_ctt>0 .or. idiag%id_ctp>0 .or. idiag%id_ctz>0).and. atm(n)%flagstruct%nwat==6) then
            allocate ( var1(isc:iec,jsc:jec) )
            allocate ( var2(isc:iec,jsc:jec) )
!$OMP parallel do default(shared) private(tmp)
            do j=jsc,jec
               do i=isc,iec
                  do k=2,npz
                     tmp = atm(n)%q(i,j,k,liq_wat)+atm(n)%q(i,j,k,rainwat)+atm(n)%q(i,j,k,ice_wat)+  &
                           atm(n)%q(i,j,k,snowwat)+atm(n)%q(i,j,k,graupel) 
                     if( tmp>5.e-6 ) then
                         a2(i,j) = atm(n)%pt(i,j,k)
                         var1(i,j) = 0.01*atm(n)%pe(i,k,j)
                         var2(i,j) = wz(i,j,k) - wz(i,j,npz+1) ! height AGL
                         exit
                     elseif( k==npz ) then
                        a2(i,j) = missing_value2
                        var1(i,j) = missing_value2
                        var2(i,j) = missing_value2
!!$                           a2(i,j) = Atm(n)%pt(i,j,k)
!!$                         var1(i,j) = 0.01*Atm(n)%pe(i,k+1,j)   ! surface pressure
                     endif
                  enddo
               enddo
            enddo
          if ( idiag%id_ctt>0 ) then
               used = send_data(idiag%id_ctt, a2, time)
               if(prt_minmax) call prt_maxmin('Cloud_top_T (K)', a2, isc, iec, jsc, jec, 0, 1, 1.)
          endif
          if ( idiag%id_ctp>0 ) then
               used = send_data(idiag%id_ctp, var1, time)
               if(prt_minmax) call prt_maxmin('Cloud_top_P (mb)', var1, isc, iec, jsc, jec, 0, 1, 1.)
          endif
          deallocate ( var1 )
          if ( idiag%id_ctz>0 ) then
               used = send_data(idiag%id_ctz, var2, time)
               if(prt_minmax) call prt_maxmin('Cloud_top_z (m)', var2, isc, iec, jsc, jec, 0, 1, 1.)
          endif
          deallocate ( var2 )
       endif

! Condensates:
       if ( idiag%id_qn>0 .or. idiag%id_qn200>0 .or. idiag%id_qn500>0 .or. idiag%id_qn850>0 ) then
!$OMP parallel do default(shared)
          do k=1,npz
          do j=jsc,jec
          do i=isc,iec
             wk(i,j,k) = 0.
          enddo
          enddo
          enddo
          if (liq_wat > 0) then
!$OMP parallel do default(shared)
             do k=1,npz
             do j=jsc,jec
             do i=isc,iec
                wk(i,j,k) = wk(i,j,k) + atm(n)%q(i,j,k,liq_wat)*atm(n)%delp(i,j,k)
             enddo
             enddo
             enddo
          endif
          if (ice_wat > 0) then
!$OMP parallel do default(shared)
             do k=1,npz
             do j=jsc,jec
             do i=isc,iec
                wk(i,j,k) = wk(i,j,k) + atm(n)%q(i,j,k,ice_wat)*atm(n)%delp(i,j,k)
             enddo
             enddo
             enddo
          endif
          if ( idiag%id_qn>0 ) used = send_data(idiag%id_qn, wk, time)
          if ( idiag%id_qn200>0 ) then
            call interpolate_vertical(isc, iec, jsc, jec, npz, 200.e2, atm(n)%peln, wk, a2)
            used=send_data(idiag%id_qn200, a2, time)
          endif
          if ( idiag%id_qn500>0 ) then
            call interpolate_vertical(isc, iec, jsc, jec, npz, 500.e2, atm(n)%peln, wk, a2)
            used=send_data(idiag%id_qn500, a2, time)
          endif
          if ( idiag%id_qn850>0 ) then
            call interpolate_vertical(isc, iec, jsc, jec, npz, 850.e2, atm(n)%peln, wk, a2)
            used=send_data(idiag%id_qn850, a2, time)
          endif
       endif
! Total 3D condensates
       if ( idiag%id_qp>0 ) then
!$OMP parallel do default(shared)
          do k=1,npz
          do j=jsc,jec
          do i=isc,iec
             wk(i,j,k) = 0.
          enddo
          enddo
          enddo
          if (rainwat > 0) then
!$OMP parallel do default(shared)
             do k=1,npz
             do j=jsc,jec
             do i=isc,iec
                wk(i,j,k) = wk(i,j,k) + atm(n)%q(i,j,k,rainwat)*atm(n)%delp(i,j,k)
             enddo
             enddo
             enddo
          endif
          if (snowwat > 0) then
!$OMP parallel do default(shared)
             do k=1,npz
             do j=jsc,jec
             do i=isc,iec
                wk(i,j,k) = wk(i,j,k) + atm(n)%q(i,j,k,snowwat)*atm(n)%delp(i,j,k)
             enddo
             enddo
             enddo
          endif
          if (graupel > 0) then
!$OMP parallel do default(shared)
             do k=1,npz
             do j=jsc,jec
             do i=isc,iec
                wk(i,j,k) = wk(i,j,k) + atm(n)%q(i,j,k,graupel)*atm(n)%delp(i,j,k)
             enddo
             enddo
             enddo
          endif
          used = send_data(idiag%id_qp, wk, time)
       endif

       if(idiag%id_us > 0 .and. idiag%id_vs > 0) then
          u2(:,:) = atm(n)%ua(isc:iec,jsc:jec,npz)
          v2(:,:) = atm(n)%va(isc:iec,jsc:jec,npz)
          do j=jsc,jec
             do i=isc,iec
                a2(i,j) = sqrt(u2(i,j)**2 + v2(i,j)**2)
             enddo
          enddo
          used=send_data(idiag%id_us, u2, time)
          used=send_data(idiag%id_vs, v2, time)
          if(prt_minmax) call prt_maxmin('Surf_wind_speed', a2, isc, iec, jsc, jec, 0, 1, 1.)
       endif

       if(idiag%id_tb > 0) then
          a2(:,:) = atm(n)%pt(isc:iec,jsc:jec,npz)
          used=send_data(idiag%id_tb, a2, time)
          if( prt_minmax )   &
          call prt_mxm('T_bot:', a2, isc, iec, jsc, jec, 0, 1, 1., atm(n)%gridstruct%area_64, atm(n)%domain)
       endif

       if(idiag%id_ua > 0) used=send_data(idiag%id_ua, atm(n)%ua(isc:iec,jsc:jec,:), time)
       if(idiag%id_va > 0) used=send_data(idiag%id_va, atm(n)%va(isc:iec,jsc:jec,:), time)

       if(idiag%id_ke > 0) then
          a2(:,:) = 0.
          do k=1,npz
          do j=jsc,jec
             do i=isc,iec
                a2(i,j) = a2(i,j) + atm(n)%delp(i,j,k)*(atm(n)%ua(i,j,k)**2+atm(n)%va(i,j,k)**2)
             enddo
          enddo
          enddo
! Mass weighted KE 
          do j=jsc,jec
             do i=isc,iec
                a2(i,j) = 0.5*a2(i,j)/(atm(n)%ps(i,j)-ptop)
             enddo
          enddo
          used=send_data(idiag%id_ke, a2, time)
          if(prt_minmax) then
             tot_mq  = g_sum( atm(n)%domain, a2, isc, iec, jsc, jec, ngc, atm(n)%gridstruct%area_64, 1) 
             if (master) write(*,*) 'SQRT(2.*KE; m/s)=', sqrt(2.*tot_mq)  
          endif
       endif


#ifdef GFS_PHYS
       if(idiag%id_delp > 0 .or. idiag%id_cape > 0 .or. idiag%id_cin > 0 .or. ((.not. atm(n)%flagstruct%hydrostatic) .and. idiag%id_pfnh > 0)) then
          do k=1,npz
            do j=jsc,jec
            do i=isc,iec         
              wk(i,j,k) = atm(n)%delp(i,j,k)*(1.-sum(atm(n)%q(i,j,k,2:atm(n)%flagstruct%nwat)))
            enddo
            enddo
          enddo
          if (idiag%id_delp > 0) used=send_data(idiag%id_delp, wk, time)
       endif

       if( ( (.not. atm(n)%flagstruct%hydrostatic) .and. idiag%id_pfnh > 0) .or. idiag%id_cape > 0 .or. idiag%id_cin > 0) then
           do k=1,npz
             do j=jsc,jec
             do i=isc,iec         
#ifdef MULTI_GASES
                 wk(i,j,k) = -wk(i,j,k)/(atm(n)%delz(i,j,k)*grav)*rdgas*  &
                             atm(n)%pt(i,j,k)*virq(atm(n)%q(i,j,k,1:num_gas))
#else
                 wk(i,j,k) = -wk(i,j,k)/(atm(n)%delz(i,j,k)*grav)*rdgas*  &
                             atm(n)%pt(i,j,k)*(1.+zvir*atm(n)%q(i,j,k,sphum))
#endif
             enddo
             enddo
           enddo
!           if (prt_minmax) then
!              call prt_maxmin(' PFNH (mb)', wk(isc:iec,jsc:jec,1), isc, iec, jsc, jec, 0, npz, 1.E-2)
!           endif
           used=send_data(idiag%id_pfnh, wk, time)
       endif
#else
       if(idiag%id_delp > 0) used=send_data(idiag%id_delp, atm(n)%delp(isc:iec,jsc:jec,:), time)

       if( (.not. atm(n)%flagstruct%hydrostatic) .and. (idiag%id_pfnh > 0 .or.  idiag%id_cape > 0 .or. idiag%id_cin > 0)) then
           do k=1,npz
             do j=jsc,jec
             do i=isc,iec
#ifdef MULTI_GASES
                 wk(i,j,k) = -atm(n)%delp(i,j,k)/(atm(n)%delz(i,j,k)*grav)*rdgas*          &
                              atm(n)%pt(i,j,k)*virq(atm(n)%q(i,j,k,1:num_gas))
#else
                 wk(i,j,k) = -atm(n)%delp(i,j,k)/(atm(n)%delz(i,j,k)*grav)*rdgas*          &
                              atm(n)%pt(i,j,k)*(1.+zvir*atm(n)%q(i,j,k,sphum))
#endif
             enddo
             enddo
           enddo
           used=send_data(idiag%id_pfnh, wk, time)
       endif
#endif

      if( atm(n)%flagstruct%hydrostatic .and. (idiag%id_pfhy > 0 .or. idiag%id_cape > 0 .or. idiag%id_cin > 0) ) then
          do k=1,npz
            do j=jsc,jec
            do i=isc,iec         
                wk(i,j,k) = 0.5 *(atm(n)%pe(i,k,j)+atm(n)%pe(i,k+1,j))
            enddo
            enddo
          enddo
          used=send_data(idiag%id_pfhy, wk, time)
      endif

       if (idiag%id_cape > 0 .or. idiag%id_cin > 0) then
          !wk here contains layer-mean pressure

          allocate(var2(isc:iec,jsc:jec))
          allocate(a3(isc:iec,jsc:jec,npz))
          call eqv_pot(a3, atm(n)%pt, atm(n)%delp, atm(n)%delz, atm(n)%peln, atm(n)%pkz, atm(n)%q(isd:ied,jsd:jed,1:npz,sphum), &
                       isc, iec, jsc, jec, ngc, npz, atm(n)%flagstruct%hydrostatic, atm(n)%flagstruct%moist_phys)

!$OMP parallel do default(shared)
          do j=jsc,jec
          do i=isc,iec
             a2(i,j) = 0.
             var2(i,j) = 0.

             call getcape(npz, wk(i,j,:), atm(n)%pt(i,j,:), -atm(n)%delz(i,j,:), atm(n)%q(i,j,:,sphum), a3(i,j,:), a2(i,j), var2(i,j), source_in=1)
          enddo
          enddo

          if (idiag%id_cape > 0) then
             if (prt_minmax) then
                call prt_maxmin(' CAPE (J/kg)', a2, isc,iec,jsc,jec, 0, 1, 1.)
             endif
             used=send_data(idiag%id_cape, a2, time)
          endif
          if (idiag%id_cin > 0) then
             if (prt_minmax) then
                call prt_maxmin(' CIN (J/kg)', var2, isc,iec,jsc,jec, 0, 1, 1.)
             endif
             used=send_data(idiag%id_cin, var2, time)
          endif

          deallocate(var2)
          deallocate(a3)

       endif


       if((.not. atm(n)%flagstruct%hydrostatic) .and. idiag%id_delz > 0) then
          do k=1,npz
            do j=jsc,jec
            do i=isc,iec
               wk(i,j,k) = -atm(n)%delz(i,j,k)
            enddo
            enddo
          enddo
          used=send_data(idiag%id_delz, wk, time)
       endif
 

! pressure for masking p-level fields
! incorrectly defines a2 to be ps (in mb).
       if (idiag%id_pmask>0) then
            do j=jsc,jec
            do i=isc,iec
                a2(i,j) = exp((atm(n)%peln(i,npz+1,j)+atm(n)%peln(i,npz+1,j))*0.5)*0.01
               !a2(i,j) = Atm(n)%delp(i,j,k)/(Atm(n)%peln(i,k+1,j)-Atm(n)%peln(i,k,j))*0.01
            enddo
            enddo
            used=send_data(idiag%id_pmask, a2, time)
       endif
! fix for pressure for masking p-level fields
! based on lowest-level pfull
! define pressure at lowest level the same as interpolate_vertical (in mb)
       if (idiag%id_pmaskv2>0) then
            do j=jsc,jec
            do i=isc,iec
                a2(i,j) = exp((atm(n)%peln(i,npz,j)+atm(n)%peln(i,npz+1,j))*0.5)*0.01
            enddo
            enddo
            used=send_data(idiag%id_pmaskv2, a2, time)
       endif

       if ( idiag%id_u100m>0 .or. idiag%id_v100m>0 .or. idiag%id_w100m>0 .or. idiag%id_w5km>0 .or. idiag%id_w2500m>0 .or. idiag%id_w1km>0 .or. idiag%id_basedbz>0 .or. idiag%id_dbz4km>0) then
          if (.not.allocated(wz)) allocate ( wz(isc:iec,jsc:jec,npz+1) )
          if ( atm(n)%flagstruct%hydrostatic) then
             rgrav = 1. / grav
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,wz,npz,Atm,n,rgrav)
            do j=jsc,jec
               do i=isc,iec
                  wz(i,j,npz+1) = 0.
!                  wz(i,j,npz+1) = Atm(n)%phis(i,j)/grav
               enddo
               do k=npz,1,-1
                  do i=isc,iec
                     wz(i,j,k) = wz(i,j,k+1) - (rdgas*rgrav)*atm(n)%pt(i,j,k)*(atm(n)%peln(i,k,j) - atm(n)%peln(i,k+1,j))
                  enddo
               enddo
            enddo
          else
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,wz,npz,Atm,n)
            do j=jsc,jec
               do i=isc,iec
                  wz(i,j,npz+1) = 0.
!                  wz(i,j,npz+1) = Atm(n)%phis(i,j)/grav
               enddo
               do k=npz,1,-1
                  do i=isc,iec
                     wz(i,j,k) = wz(i,j,k+1) - atm(n)%delz(i,j,k)
                  enddo
               enddo
            enddo
         endif
            if( prt_minmax )   &
                 call prt_maxmin('ZTOP', wz(isc:iec,jsc:jec,1)+atm(n)%phis(isc:iec,jsc:jec)/grav, isc, iec, jsc, jec, 0, 1, 1.e-3)
       endif

       if ( idiag%id_rain5km>0 ) then
            rainwat = get_tracer_index(model_atmos, 'rainwat')
            call interpolate_z(isc, iec, jsc, jec, npz, 5.e3, wz, atm(n)%q(isc:iec,jsc:jec,:,rainwat), a2)
            used=send_data(idiag%id_rain5km, a2, time)
            if(prt_minmax) call prt_maxmin('rain5km', a2, isc, iec, jsc, jec, 0, 1, 1.)
       endif
       if ( idiag%id_w5km>0 ) then
            call interpolate_z(isc, iec, jsc, jec, npz, 5.e3, wz, atm(n)%w(isc:iec,jsc:jec,:), a2)
            used=send_data(idiag%id_w5km, a2, time)
            if(prt_minmax) call prt_maxmin('W5km', a2, isc, iec, jsc, jec, 0, 1, 1.)
       endif
       if ( idiag%id_w2500m>0 ) then
            call interpolate_z(isc, iec, jsc, jec, npz, 2.5e3, wz, atm(n)%w(isc:iec,jsc:jec,:), a2)
            used=send_data(idiag%id_w2500m, a2, time)
            if(prt_minmax) call prt_maxmin('W2500m', a2, isc, iec, jsc, jec, 0, 1, 1.)
       endif
       if ( idiag%id_w1km>0 ) then
            call interpolate_z(isc, iec, jsc, jec, npz, 1.e3, wz, atm(n)%w(isc:iec,jsc:jec,:), a2)
            used=send_data(idiag%id_w1km, a2, time)
            if(prt_minmax) call prt_maxmin('W1km', a2, isc, iec, jsc, jec, 0, 1, 1.)
       endif
       if ( idiag%id_w100m>0 ) then
            call interpolate_z(isc, iec, jsc, jec, npz, 100., wz, atm(n)%w(isc:iec,jsc:jec,:), a2)
            used=send_data(idiag%id_w100m, a2, time)
            if(prt_minmax) call prt_maxmin('w100m', a2, isc, iec, jsc, jec, 0, 1, 1.)
       endif
       if ( idiag%id_u100m>0 ) then
            call interpolate_z(isc, iec, jsc, jec, npz, 100., wz, atm(n)%ua(isc:iec,jsc:jec,:), a2)
            used=send_data(idiag%id_u100m, a2, time)
            if(prt_minmax) call prt_maxmin('u100m', a2, isc, iec, jsc, jec, 0, 1, 1.)
       endif
       if ( idiag%id_v100m>0 ) then
            call interpolate_z(isc, iec, jsc, jec, npz, 100., wz, atm(n)%va(isc:iec,jsc:jec,:), a2)
            used=send_data(idiag%id_v100m, a2, time)
            if(prt_minmax) call prt_maxmin('v100m', a2, isc, iec, jsc, jec, 0, 1, 1.)
       endif

       if ( rainwat > 0 .and. (idiag%id_dbz>0 .or. idiag%id_maxdbz>0 .or. idiag%id_basedbz>0 .or. idiag%id_dbz4km>0 .or. idiag%id_dbztop>0 .or. idiag%id_dbz_m10C>0)) then

          if (.not. allocated(a3)) allocate(a3(isc:iec,jsc:jec,npz))

!          call dbzcalc_smithxue(Atm(n)%q, Atm(n)%pt, Atm(n)%delp, Atm(n)%peln, Atm(n)%delz, &
          call dbzcalc(atm(n)%q, atm(n)%pt, atm(n)%delp, atm(n)%peln, atm(n)%delz, &
               a3, a2, allmax, atm(n)%bd, npz, atm(n)%ncnst, atm(n)%flagstruct%hydrostatic, &
               zvir, .false., .false., .false., .true. ) ! GFDL MP has constant N_0 intercept

          if (idiag%id_dbz > 0) used=send_data(idiag%id_dbz, a3, time)
          if (idiag%id_maxdbz > 0) used=send_data(idiag%id_maxdbz, a2, time)

          if (idiag%id_basedbz > 0) then
             !interpolate to 1km dbz
             call cs_interpolator(isc, iec, jsc, jec, npz, a3, 1000., wz, a2, -20.)
             used=send_data(idiag%id_basedbz, a2, time)
             if(prt_minmax) call prt_maxmin('Base_dBz', a2, isc, iec, jsc, jec, 0, 1, 1.)
          endif

          if (idiag%id_dbz4km > 0) then
             !interpolate to 1km dbz
             call cs_interpolator(isc, iec, jsc, jec, npz, a3, 4000., wz, a2, -20.)
             used=send_data(idiag%id_dbz4km, a2, time)
          endif
          if (idiag%id_dbztop > 0) then
             do j=jsc,jec
             do i=isc,iec
                a2(i,j) = missing_value
             do k=2,npz
                if (wz(i,j,k) >= 25000. ) continue ! nothing above 25 km
                if (a3(i,j,k) >= 18.5 ) then
                   a2(i,j) = wz(i,j,k)
                   exit
                endif
             enddo
             enddo
             enddo
             used=send_data(idiag%id_dbztop, a2, time)
          endif
          if (idiag%id_dbz_m10C > 0) then
             do j=jsc,jec
             do i=isc,iec
                a2(i,j) = missing_value
             do k=npz,1,-1
                if (wz(i,j,k) >= 25000. ) exit ! nothing above 25 km
                if (atm(n)%pt(i,j,k) <= 263.14 ) then
                   a2(i,j) = a3(i,j,k)
                   exit
                endif
             enddo
             enddo
             enddo
             used=send_data(idiag%id_dbz_m10C, a2, time)
          endif

          if (prt_minmax) then
             call mpp_max(allmax)
             if (master) write(*,*) 'max reflectivity = ', allmax, ' dBZ'
          endif

          deallocate(a3)
       endif
!-------------------------------------------------------
! Applying cubic-spline as the intepolator for (u,v,T,q)
!-------------------------------------------------------
       if(.not. allocated(a3)) allocate( a3(isc:iec,jsc:jec,nplev) )
! u-winds:
       idg(:) = idiag%id_u(:)

       do_cs_intp = .false.
       do i=1,nplev
          if ( idg(i)>0 ) then
               do_cs_intp = .true.
               exit
          endif
       enddo

       if ( do_cs_intp ) then
          call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%ua(isc:iec,jsc:jec,:), nplev,    &
                               pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), idg, a3, -1)
!                              plevs, Atm(n)%peln, idg, a3, -1)
          do i=1,nplev
             if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
          enddo
       endif

       if (idiag%id_u_plev>0) then
         id1(:) = 1
         call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%ua(isc:iec,jsc:jec,:), nplev,    &
                              pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), id1, a3, -1)
         used=send_data(idiag%id_u_plev, a3(isc:iec,jsc:jec,:), time)
       endif

! v-winds:
       idg(:) = idiag%id_v(:)

       do_cs_intp = .false.
       do i=1,nplev
          if ( idg(i)>0 ) then
               do_cs_intp = .true.
               exit
          endif
       enddo

       if ( do_cs_intp ) then
          call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%va(isc:iec,jsc:jec,:), nplev,    &
                               pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), idg, a3, -1)
!                              plevs, Atm(n)%peln, idg, a3, -1)
          do i=1,nplev
             if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
          enddo
       endif

       if (idiag%id_v_plev>0) then
         id1(:) = 1
         call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%va(isc:iec,jsc:jec,:), nplev,    &
                              pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), id1, a3, -1)
         used=send_data(idiag%id_v_plev, a3(isc:iec,jsc:jec,:), time)
       endif

! Specific humidity
       idg(:) = idiag%id_q(:)

       do_cs_intp = .false.
       do i=1,nplev
          if ( idg(i)>0 ) then
               do_cs_intp = .true.
               exit
          endif
       enddo

       if ( do_cs_intp ) then
          call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%q(isc:iec,jsc:jec,:,sphum), nplev, &
                               pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), idg, a3, 0)
!                              plevs, Atm(n)%peln, idg, a3, 0)
          do i=1,nplev
             if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
          enddo
       endif

       if (idiag%id_q_plev>0) then
         id1(:) = 1
         call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%q(isc:iec,jsc:jec,:,sphum), nplev, &
                              pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), id1, a3, 0)
         used=send_data(idiag%id_q_plev, a3(isc:iec,jsc:jec,:), time)
       endif

! Omega
       idg(:) = idiag%id_omg(:)

       do_cs_intp = .false.
       do i=1,nplev
          if ( idg(i)>0 ) then
               do_cs_intp = .true.
               exit
          endif
       enddo
       if ( do_cs_intp ) then
          call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%omga(isc:iec,jsc:jec,:), nplev,    &
                               pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), idg, a3, -1)
!                              plevs, Atm(n)%peln, idg, a3)
          do i=1,nplev
             if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
          enddo
       endif

       if (idiag%id_omg_plev>0) then
         id1(:) = 1
         call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%omga(isc:iec,jsc:jec,:), nplev,    &
                              pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), id1, a3, -1)
         used=send_data(idiag%id_omg_plev, a3(isc:iec,jsc:jec,:), time)
       endif

       if( allocated(a3) ) deallocate (a3)
! *** End cs_intp

       if ( idiag%id_sl12>0 ) then   ! 13th level wind speed (~ 222 mb for the 32L setup)
            do j=jsc,jec
               do i=isc,iec
                  a2(i,j) = sqrt(atm(n)%ua(i,j,12)**2 + atm(n)%va(i,j,12)**2)
               enddo
            enddo
            used=send_data(idiag%id_sl12, a2, time)
       endif
       if ( idiag%id_sl13>0 ) then   ! 13th level wind speed (~ 222 mb for the 32L setup)
            do j=jsc,jec
               do i=isc,iec
                  a2(i,j) = sqrt(atm(n)%ua(i,j,13)**2 + atm(n)%va(i,j,13)**2)
               enddo
            enddo
            used=send_data(idiag%id_sl13, a2, time)
       endif

       if ( (.not.atm(n)%flagstruct%hydrostatic) .and. idiag%id_w200>0 ) then
            call interpolate_vertical(isc, iec, jsc, jec, npz,   &
                                      200.e2, atm(n)%peln, atm(n)%w(isc:iec,jsc:jec,:), a2)
            used=send_data(idiag%id_w200, a2, time)
       endif
! 500-mb
       if ( (.not.atm(n)%flagstruct%hydrostatic) .and. idiag%id_w500>0 ) then
            call interpolate_vertical(isc, iec, jsc, jec, npz,   &
                                      500.e2, atm(n)%peln, atm(n)%w(isc:iec,jsc:jec,:), a2)
            used=send_data(idiag%id_w500, a2, time)
       endif
       if ( (.not.atm(n)%flagstruct%hydrostatic) .and. idiag%id_w700>0 ) then
            call interpolate_vertical(isc, iec, jsc, jec, npz,   &
                                      700.e2, atm(n)%peln, atm(n)%w(isc:iec,jsc:jec,:), a2)
            used=send_data(idiag%id_w700, a2, time)
       endif
       if ( (.not.atm(n)%flagstruct%hydrostatic) .and. idiag%id_w850>0 .or. idiag%id_x850>0) then
            call interpolate_vertical(isc, iec, jsc, jec, npz,   &
                                      850.e2, atm(n)%peln, atm(n)%w(isc:iec,jsc:jec,:), a2)
            used=send_data(idiag%id_w850, a2, time)

            if ( idiag%id_x850>0 .and. idiag%id_vort850>0 ) then
                 x850(:,:) = x850(:,:)*a2(:,:) 
                 used=send_data(idiag%id_x850, x850, time)
                 deallocate ( x850 )
            endif
       endif


       if ( .not.atm(n)%flagstruct%hydrostatic .and. idiag%id_w>0  ) then
          used=send_data(idiag%id_w, atm(n)%w(isc:iec,jsc:jec,:), time)
       endif
       if ( .not. atm(n)%flagstruct%hydrostatic .and. (idiag%id_wmaxup>0 .or. idiag%id_wmaxdn>0) ) then
          allocate(var2(isc:iec,jsc:jec))
          do j=jsc,jec
          do i=isc,iec
             a2(i,j) = 0.
             var2(i,j) = 0.
             do k=3,npz
                if (atm(n)%pe(i,k,j) <= 400.e2) continue
                a2(i,j) = max(a2(i,j),atm(n)%w(i,j,k))
                var2(i,j) = min(var2(i,j),atm(n)%w(i,j,k))
             enddo
          enddo
          enddo
          if (idiag%id_wmaxup > 0) then
             used=send_data(idiag%id_wmaxup, a2, time)
          endif
          if (idiag%id_wmaxdn > 0) then
             used=send_data(idiag%id_wmaxdn, var2, time)
          endif
          deallocate(var2)
       endif

       if(idiag%id_pt   > 0) used=send_data(idiag%id_pt  , atm(n)%pt  (isc:iec,jsc:jec,:), time)
       if(idiag%id_omga > 0) used=send_data(idiag%id_omga, atm(n)%omga(isc:iec,jsc:jec,:), time)
       if(idiag%id_diss > 0) used=send_data(idiag%id_diss, atm(n)%diss_est(isc:iec,jsc:jec,:), time)
       
       allocate( a3(isc:iec,jsc:jec,npz) )
       if(idiag%id_theta_e > 0 ) then

          if ( atm(n)%flagstruct%adiabatic .and. atm(n)%flagstruct%kord_tm>0 ) then
             do k=1,npz
                do j=jsc,jec
                   do i=isc,iec
                      a3(i,j,k) = atm(n)%pt(i,j,k)
                   enddo
                enddo
             enddo
          else
             call eqv_pot(a3, atm(n)%pt, atm(n)%delp, atm(n)%delz, atm(n)%peln, atm(n)%pkz, atm(n)%q(isd:ied,jsd:jed,1:npz,sphum),    &
                  isc, iec, jsc, jec, ngc, npz, atm(n)%flagstruct%hydrostatic, atm(n)%flagstruct%moist_phys)
          endif

          if (idiag%id_theta_e > 0) then
             if( prt_minmax ) call prt_maxmin('Theta_E', a3, isc, iec, jsc, jec, 0, npz, 1.)
             used=send_data(idiag%id_theta_e, a3, time)
          end if
          theta_d = get_tracer_index(model_atmos, 'theta_d')
          if ( theta_d>0 ) then
             if( prt_minmax ) then
                ! Check level-34 ~ 300 mb
                a2(:,:) = 0.
                do k=1,npz
                do j=jsc,jec
                do i=isc,iec
                   a2(i,j) = a2(i,j) + atm(n)%delp(i,j,k)*(atm(n)%q(i,j,k,theta_d)-a3(i,j,k))**2
                enddo
                enddo
                enddo
                call prt_mxm('PT_SUM', a2, isc, iec, jsc, jec, 0, 1, 1.e-5, atm(n)%gridstruct%area_64, atm(n)%domain)

                do k=1,npz
                do j=jsc,jec
                do i=isc,iec
                   a3(i,j,k) =  atm(n)%q(i,j,k,theta_d)/a3(i,j,k) - 1.
                enddo
                enddo
                enddo
                call prt_maxmin('Theta_Err (%)', a3, isc, iec, jsc, jec, 0, npz, 100.)
                !           if ( master ) write(*,*) 'PK0=', pk0, 'KAPPA=', kappa
             endif
          endif

       endif

       if(idiag%id_ppt> 0) then
! Potential temperature perturbation for gravity wave test_case
          allocate ( idiag%pt1(npz) )
          if( .not. allocated(a3) ) allocate ( a3(isc:iec,jsc:jec,npz) )
#ifdef TEST_GWAVES
          call gw_1d(npz, 1000.e2, atm(n)%ak, atm(n)%ak, atm(n)%ak(1), 10.e3, idiag%pt1)
#else
          idiag%pt1 = 0. 
#endif
          do k=1,npz
          do j=jsc,jec
             do i=isc,iec
!               wk(i,j,k) =  (Atm(n)%pt(i,j,k)-300.)/Atm(n)%pkz(i,j,k) * pk0
                wk(i,j,k) =  (atm(n)%pt(i,j,k)/atm(n)%pkz(i,j,k) - idiag%pt1(k)) * pk0
             enddo
          enddo
          enddo
          used=send_data(idiag%id_ppt, wk, time)

          if( prt_minmax ) then
              call prt_maxmin('PoTemp', wk, isc, iec, jsc, jec, 0, npz, 1.)
          endif

          if( allocated(a3) ) deallocate ( a3 )
          deallocate ( idiag%pt1 )
       endif


#ifndef SW_DYNAMICS
        do itrac=1, atm(n)%ncnst
          call get_tracer_names (model_atmos, itrac, tname)
          if (idiag%id_tracer(itrac) > 0 .and. itrac.gt.nq) then
            used = send_data(idiag%id_tracer(itrac), atm(n)%qdiag(isc:iec,jsc:jec,:,itrac), time )
          else
            used = send_data(idiag%id_tracer(itrac), atm(n)%q(isc:iec,jsc:jec,:,itrac), time )
          endif
          if (itrac .le. nq) then
            if( prt_minmax ) call prt_maxmin(trim(tname), atm(n)%q(:,:,1,itrac),     &
                              isc, iec, jsc, jec, ngc, npz, 1.)
          else
            if( prt_minmax ) call prt_maxmin(trim(tname), atm(n)%qdiag(:,:,1,itrac), &
                              isc, iec, jsc, jec, ngc, npz, 1.)
          endif
!-------------------------------
! ESM TRACER diagnostics output:
! jgj: per SJ email (jul 17 2008): q_dry = q_moist/(1-sphum)
! mass mixing ratio: q_dry = mass_tracer/mass_dryair = mass_tracer/(mass_air - mass_water) ~ q_moist/(1-sphum)
! co2_mmr = (wco2/wair) * co2_vmr
! Note: There is a check to ensure tracer number one is sphum

          if (idiag%id_tracer_dmmr(itrac) > 0 .or. idiag%id_tracer_dvmr(itrac) > 0) then
              if (itrac .gt. nq) then
                dmmr(:,:,:) = atm(n)%qdiag(isc:iec,jsc:jec,1:npz,itrac)  &
                              /(1.0-atm(n)%q(isc:iec,jsc:jec,1:npz,1))
              else
                dmmr(:,:,:) = atm(n)%q(isc:iec,jsc:jec,1:npz,itrac)  &
                              /(1.0-atm(n)%q(isc:iec,jsc:jec,1:npz,1))
              endif
              dvmr(:,:,:) = dmmr(isc:iec,jsc:jec,1:npz) * wtmair/idiag%w_mr(itrac)
              used = send_data(idiag%id_tracer_dmmr(itrac), dmmr, time )
              used = send_data(idiag%id_tracer_dvmr(itrac), dvmr, time )
              if( prt_minmax ) then
                 call prt_maxmin(trim(tname)//'_dmmr', dmmr, &
                    isc, iec, jsc, jec, 0, npz, 1.)
                 call prt_maxmin(trim(tname)//'_dvmr', dvmr, & 
                    isc, iec, jsc, jec, 0, npz, 1.)
            endif
          endif
        enddo

! Maximum overlap cloud fraction
      if ( .not. atm(n)%neststruct%nested )  then
        if ( cld_amt > 0 .and. prt_minmax ) then
          a2(:,:) = 0.
          do k=1,npz
             do j=jsc,jec
             do i=isc,iec
                a2(i,j) =  max(a2(i,j), atm(n)%q(i,j,k,cld_amt) )
             enddo
             enddo
          enddo
          call prt_gb_nh_sh('Max_cld GB_NH_SH_EQ',isc,iec, jsc,jec, a2, atm(n)%gridstruct%area_64(isc:iec,jsc:jec),   &
                            atm(n)%gridstruct%agrid_64(isc:iec,jsc:jec,2))
        endif
      endif

#endif

   ! enddo  ! end ntileMe do-loop

    deallocate ( a2 )
    deallocate ( u2 )
    deallocate ( v2 )
    deallocate ( wk )

    if (allocated(a3)) deallocate(a3)
    if (allocated(wz)) deallocate(wz)
    if (allocated(dmmr)) deallocate(dmmr)
    if (allocated(dvmr)) deallocate(dvmr)

    call nullify_domain()


 end subroutine fv_diag

 subroutine wind_max(isc, iec, jsc, jec ,isd, ied, jsd, jed, us, vs, ws_max, domain)
 integer isc, iec, jsc, jec
 integer isd, ied, jsd, jed
 real, intent(in), dimension(isc:iec,jsc:jec):: us, vs
 real, intent(out) :: ws_max(isc:iec,jsc:jec)
 type(domain2d), intent(INOUT) :: domain
! Local
 real :: wx(isc:iec,jsd:jed), ws(isd:ied,jsd:jed)
 integer:: i,j

 ws = 0.   ! fill corners with zeros
 do j=jsc,jec
    do i=isc,iec
       ws(i,j) = sqrt(us(i,j)**2 + vs(i,j)**2)
    enddo
 enddo

 call mpp_update_domains( ws, domain )

 do j=jsd,jed
    do i=isc,iec
       wx(i,j) = max(ws(i-3,j), ws(i-2,j), ws(i-1,j), ws(i,j), ws(i+1,j), ws(i+2,j), ws(i+3,j))
    enddo
 enddo

 do j=jsc,jec
    do i=isc,iec
       ws_max(i,j) = max(wx(i,j-3), wx(i,j-2), wx(i,j-1), wx(i,j), wx(i,j+1), wx(i,j+2), wx(i,j+3))
    enddo
 enddo

 end subroutine wind_max


 subroutine get_vorticity(isc, iec, jsc, jec ,isd, ied, jsd, jed, npz, u, v, vort, dx, dy, rarea)
 integer isd, ied, jsd, jed, npz
 integer isc, iec, jsc, jec
 real, intent(in)  :: u(isd:ied, jsd:jed+1, npz), v(isd:ied+1, jsd:jed, npz)
 real, intent(out) :: vort(isc:iec, jsc:jec, npz)
 real, intent(IN) :: dx(isd:ied,jsd:jed+1)
 real, intent(IN) :: dy(isd:ied+1,jsd:jed)
 real, intent(IN) :: rarea(isd:ied,jsd:jed)
! Local
 real :: utmp(isc:iec, jsc:jec+1), vtmp(isc:iec+1, jsc:jec)
 integer :: i,j,k

      do k=1,npz
         do j=jsc,jec+1
            do i=isc,iec
               utmp(i,j) = u(i,j,k)*dx(i,j)
            enddo
         enddo
         do j=jsc,jec
            do i=isc,iec+1
               vtmp(i,j) = v(i,j,k)*dy(i,j)
            enddo
         enddo

         do j=jsc,jec
            do i=isc,iec
               vort(i,j,k) = rarea(i,j)*(utmp(i,j)-utmp(i,j+1)-vtmp(i,j)+vtmp(i+1,j))
            enddo
         enddo
      enddo

 end subroutine get_vorticity


 subroutine get_height_field(is, ie, js, je, ng, km, hydrostatic, delz, wz, pt, q, peln, zvir)
  integer, intent(in):: is, ie, js, je, km, ng
  real, intent(in):: peln(is:ie,km+1,js:je)
  real, intent(in):: pt(is-ng:ie+ng,js-ng:je+ng,km)
  real, intent(in)::  q(is-ng:ie+ng,js-ng:je+ng,km,*) ! water vapor
  real, intent(in):: delz(is-ng:,js-ng:,1:)
  real, intent(in):: zvir
  logical, intent(in):: hydrostatic
  real, intent(out):: wz(is:ie,js:je,km+1)
!
  integer i,j,k
  real gg

      gg  = rdgas * ginv

      do j=js,je
         do i=is,ie
            wz(i,j,km+1) = idiag%zsurf(i,j)
         enddo
      if (hydrostatic ) then
         do k=km,1,-1
            do i=is,ie
#ifdef MULTI_GASES
               wz(i,j,k) = wz(i,j,k+1) + gg*pt(i,j,k)*virq(q(i,j,k,1:num_gas))  &
                          *(peln(i,k+1,j)-peln(i,k,j))
#else
               wz(i,j,k) = wz(i,j,k+1) + gg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))  &
                          *(peln(i,k+1,j)-peln(i,k,j))
#endif
            enddo
         enddo
      else
         do k=km,1,-1
            do i=is,ie
               wz(i,j,k) = wz(i,j,k+1)  - delz(i,j,k)
            enddo
         enddo
      endif
      enddo

 end subroutine get_height_field

 subroutine range_check(qname, q, is, ie, js, je, n_g, km, pos, q_low, q_hi, bad_range)
      character(len=*), intent(in)::  qname
      integer, intent(in):: is, ie, js, je
      integer, intent(in):: n_g, km
      real, intent(in)::    q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
      real, intent(in):: pos(is-n_g:ie+n_g, js-n_g:je+n_g,2)
      real, intent(in):: q_low, q_hi
      logical, optional, intent(out):: bad_range
!
      real qmin, qmax
      integer i,j,k

      if ( present(bad_range) ) bad_range = .false. 
      qmin = q(is,js,1)
      qmax = qmin

      do k=1,km
      do j=js,je
         do i=is,ie
            if( q(i,j,k) < qmin ) then
                qmin = q(i,j,k)
            elseif( q(i,j,k) > qmax ) then
                qmax = q(i,j,k)
            endif
          enddo
      enddo
      enddo

      call mp_reduce_min(qmin)
      call mp_reduce_max(qmax)

      if( qmin<q_low .or. qmax>q_hi ) then
          if(master) write(*,*) 'Range_check Warning:', qname, ' max = ', qmax, ' min = ', qmin
          if ( present(bad_range) ) then
               bad_range = .true. 
          endif
      endif

      if ( present(bad_range) ) then
! Print out where the bad value(s) is (are)
         if ( bad_range .EQV. .false. ) return
         do k=1,km
            do j=js,je
               do i=is,ie
                  if( q(i,j,k)<q_low .or. q(i,j,k)>q_hi ) then
                      write(*,*) 'Warn_K=',k,'(i,j)=',i,j, pos(i,j,1)*rad2deg, pos(i,j,2)*rad2deg, q(i,j,k)
                      if ( k/= 1 ) write(*,*) k-1, q(i,j,k-1)
                      if ( k/=km ) write(*,*) k+1, q(i,j,k+1)
                  endif
               enddo
            enddo
         enddo
         call mpp_error(note,'==> Error from range_check: data out of bound')
      endif

 end subroutine range_check

 subroutine prt_maxmin(qname, q, is, ie, js, je, n_g, km, fac)
      character(len=*), intent(in)::  qname
      integer, intent(in):: is, ie, js, je
      integer, intent(in):: n_g, km
      real, intent(in)::    q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
      real, intent(in)::    fac

      real qmin, qmax
      integer i,j,k
      !mpp_root_pe doesn't appear to recognize nested grid
      master = (mpp_pe()==mpp_root_pe()) .or. is_master()

      qmin = q(is,js,1)
      qmax = qmin

      do k=1,km
      do j=js,je
         do i=is,ie
!           qmin = min(qmin, q(i,j,k))
!           qmax = max(qmax, q(i,j,k))
            if( q(i,j,k) < qmin  ) then
                qmin = q(i,j,k)
            elseif( q(i,j,k) > qmax ) then
                qmax = q(i,j,k)
            endif
          enddo
      enddo
      enddo

      call mp_reduce_min(qmin)
      call mp_reduce_max(qmax)

      if(master) then
            write(*,*) qname//trim(gn), ' max = ', qmax*fac, ' min = ', qmin*fac
      endif

 end subroutine prt_maxmin

 subroutine prt_mxm(qname, q, is, ie, js, je, n_g, km, fac, area, domain)
      character(len=*), intent(in)::  qname
      integer, intent(in):: is, ie, js, je
      integer, intent(in):: n_g, km
      real, intent(in)::    q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
      real, intent(in)::    fac
! BUG !!!
!     real, intent(IN)::    area(is-n_g:ie+n_g, js-n_g:je+n_g, km)
      real(kind=R_GRID), intent(IN)::    area(is-3:ie+3, js-3:je+3)
      type(domain2d), intent(INOUT) :: domain
!
      real qmin, qmax, gmean
      integer i,j,k

      !mpp_root_pe doesn't appear to recognize nested grid
      master = (mpp_pe()==mpp_root_pe()) .or. is_master()
      qmin = q(is,js,1)
      qmax = qmin
      gmean = 0.

      do k=1,km
      do j=js,je
         do i=is,ie
!           qmin = min(qmin, q(i,j,k))
!           qmax = max(qmax, q(i,j,k))
            if( q(i,j,k) < qmin ) then
                qmin = q(i,j,k)
            elseif( q(i,j,k) > qmax ) then
                qmax = q(i,j,k)
            endif
          enddo
      enddo
      enddo

      call mp_reduce_min(qmin)
      call mp_reduce_max(qmax)

! SJL: BUG!!!
!     gmean = g_sum(domain, q(is,js,km), is, ie, js, je, 3, area, 1) 
      gmean = g_sum(domain, q(is:ie,js:je,km), is, ie, js, je, 3, area, 1) 

      if(master) write(6,*) qname, gn, qmax*fac, qmin*fac, gmean*fac

 end subroutine prt_mxm

 !Added nwat == 1 case for water vapor diagnostics
 subroutine prt_mass(km, nq, is, ie, js, je, n_g, nwat, ps, delp, q, area, domain)

 integer, intent(in):: is, ie, js, je
 integer, intent(in):: nq, n_g, km, nwat
 real, intent(in)::   ps(is-n_g:ie+n_g, js-n_g:je+n_g)
 real, intent(in):: delp(is-n_g:ie+n_g, js-n_g:je+n_g, km)
 real, intent(in)::  q(is-n_g:ie+n_g, js-n_g:je+n_g, km, nq)
 real(kind=R_GRID), intent(IN):: area(is-n_g:ie+n_g,js-n_g:je+n_g)
 type(domain2d), intent(INOUT) :: domain
! Local:
 real psq(is:ie,js:je,nwat), psqv(is:ie,js:je)
 real q_strat(is:ie,js:je)
 real qtot(nwat), qwat
 real psmo, totw, psdry
 integer k, n, kstrat

!Needed when calling prt_mass in fv_restart?
    sphum   = get_tracer_index(model_atmos, 'sphum')
    liq_wat = get_tracer_index(model_atmos, 'liq_wat')
    ice_wat = get_tracer_index(model_atmos, 'ice_wat')

    rainwat = get_tracer_index(model_atmos, 'rainwat')
    snowwat = get_tracer_index(model_atmos, 'snowwat')
    graupel = get_tracer_index(model_atmos, 'graupel')

 if ( nwat==0 ) then
      psmo = g_sum(domain, ps(is:ie,js:je), is, ie, js, je, n_g, area, 1) 
      if( master ) write(*,*) 'Total surface pressure (mb)', trim(gn), ' = ',  0.01*psmo
      call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,1  ), psqv(is,js)) 
      return
 endif

 psq(:,:,:) = 0.
 call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,sphum  ), psq(is,js,sphum  )) 

 if (liq_wat > 0)  &
      call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,liq_wat), psq(is,js,liq_wat))

 if (rainwat > 0)  &
      call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,rainwat), psq(is,js,rainwat))

!nwat == 4 => KESSLER, ice is probably garbage...
 if (ice_wat > 0)  &
      call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,ice_wat), psq(is,js,ice_wat))

 if (snowwat > 0) &
      call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,snowwat), psq(is,js,snowwat))
 if (graupel > 0) &
      call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,graupel), psq(is,js,graupel))


! Mean water vapor in the "stratosphere" (75 mb and above):
 if ( idiag%phalf(2)< 75. ) then
 kstrat = 1
 do k=1,km
    if ( idiag%phalf(k+1) > 75. ) exit
    kstrat = k
 enddo
 call z_sum(is, ie, js, je, kstrat, n_g, delp, q(is-n_g,js-n_g,1,sphum), q_strat(is,js)) 
 psmo = g_sum(domain, q_strat(is,js), is, ie, js, je, n_g, area, 1) * 1.e6           &
      / p_sum(is, ie, js, je, kstrat, n_g, delp, area, domain)
 if(master) write(*,*) 'Mean specific humidity (mg/kg) above 75 mb', trim(gn), '=', psmo
 endif


!-------------------
! Check global means
!-------------------
 psmo = g_sum(domain, ps(is:ie,js:je), is, ie, js, je, n_g, area, 1) 

 do n=1,nwat
    qtot(n) = g_sum(domain, psq(is,js,n), is, ie, js, je, n_g, area, 1) 
 enddo

 totw  = sum(qtot(1:nwat))
 psdry = psmo - totw

 if( master ) then
     write(*,*) 'Total surface pressure (mb)', trim(gn), ' = ',  0.01*psmo
     write(*,*) 'mean dry surface pressure', trim(gn), ' = ',    0.01*psdry
     write(*,*) 'Total Water Vapor (kg/m**2)', trim(gn), ' =',  qtot(sphum)*ginv
     if ( nwat> 2 ) then
          write(*,*) '--- Micro Phys water substances (kg/m**2) ---'
          write(*,*) 'Total cloud water', trim(gn), '=', qtot(liq_wat)*ginv
          if (rainwat > 0) &
               write(*,*) 'Total rain  water', trim(gn), '=', qtot(rainwat)*ginv
          if (ice_wat > 0) &
               write(*,*) 'Total cloud ice  ', trim(gn), '=', qtot(ice_wat)*ginv
          if (snowwat > 0) &
               write(*,*) 'Total snow       ', trim(gn), '=', qtot(snowwat)*ginv
          if (graupel > 0) &
               write(*,*) 'Total graupel    ', trim(gn), '=', qtot(graupel)*ginv
          write(*,*) '---------------------------------------------'
     elseif ( nwat==2 ) then
          write(*,*) 'GFS condensate (kg/m^2)', trim(gn), '=', qtot(liq_wat)*ginv
     endif
  endif

 end subroutine prt_mass


 subroutine z_sum(is, ie, js, je, km, n_g, delp, q, sum2)
 integer, intent(in):: is, ie, js, je,  n_g, km
 real, intent(in):: delp(is-n_g:ie+n_g, js-n_g:je+n_g, km)
 real, intent(in)::    q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
 real, intent(out):: sum2(is:ie,js:je)

 integer i,j,k

 do j=js,je
    do i=is,ie
       sum2(i,j) = delp(i,j,1)*q(i,j,1)
    enddo
    do k=2,km
       do i=is,ie
          sum2(i,j) = sum2(i,j) + delp(i,j,k)*q(i,j,k)
       enddo
    enddo
 enddo

 end subroutine z_sum


 real function p_sum(is, ie, js, je, km, n_g, delp, area, domain)
 integer, intent(in):: is, ie, js, je,  n_g, km
 real, intent(in):: delp(is-n_g:ie+n_g, js-n_g:je+n_g, km)
 real(kind=R_GRID), intent(IN) :: area(is-n_g:ie+n_g, js-n_g:je+n_g)
 real :: sum2(is:ie,js:je)
 integer i,j,k
 type(domain2d), intent(INOUT) :: domain

!$OMP parallel do default(none) shared(is,ie,js,je,km,sum2,delp)
 do j=js,je
    do i=is,ie
       sum2(i,j) = delp(i,j,1)
    enddo
    do k=2,km
       do i=is,ie
          sum2(i,j) = sum2(i,j) + delp(i,j,k)
       enddo
    enddo
 enddo
 p_sum = g_sum(domain, sum2, is, ie, js, je, n_g, area, 1)

 end function p_sum



 subroutine get_pressure_given_height(is, ie, js, je, ng, km, wz, kd, height,   &
                                      ts, peln, a2, fac)

 integer,  intent(in):: is, ie, js, je, km, ng
 integer,  intent(in):: kd
 real, intent(in):: wz(is:ie,js:je,km+1)
 real, intent(in):: ts(is-ng:ie+ng,js-ng:je+ng)
 real, intent(in):: peln(is:ie,km+1,js:je)
 real, intent(in):: height(kd)
 real, intent(out):: a2(is:ie,js:je,kd)
 real, optional, intent(in):: fac

! local:
 integer n, i,j,k
 real ptmp, tm


 do n=1,kd

!$OMP parallel do default(none) shared(is,ie,js,je,n,height,wz,km,peln,a2,ginv,ts,fac) &
!$OMP                          private(ptmp, tm)
    do j=js,je

       do 1000 i=is,ie

         if ( height(n) >= wz(i,j,km+1) ) then
!---------------------
! Search from top down
!---------------------
          do k=1,km
             if( height(n) < wz(i,j,k) .and. height(n) >= wz(i,j,k+1) ) then
! Found it!
                 ptmp = peln(i,k,j) + (peln(i,k+1,j)-peln(i,k,j)) *   &
                       (wz(i,j,k)-height(n)) / (wz(i,j,k)-wz(i,j,k+1))
                 a2(i,j,n) = exp(ptmp)
                 go to 500
             endif
          enddo

         else
!-----------------------------------------
! xtrapolation: mean laspe rate 6.5 deg/km
!-----------------------------------------
                tm = rdgas*ginv*(ts(i,j) + 3.25e-3*(wz(i,j,km)-height(n)))
          a2(i,j,n) = exp( peln(i,km+1,j) + (wz(i,j,km+1) - height(n))/tm )
         endif
500      if ( present(fac) ) a2(i,j,n) = fac * a2(i,j,n)
1000   continue
    enddo
 enddo

 end subroutine get_pressure_given_height


 subroutine get_height_given_pressure(is, ie, js, je, km, wz, kd, id, log_p, peln, a2)
 integer,  intent(in):: is, ie, js, je, km
 integer,  intent(in):: kd
 integer,  intent(in):: id(kd)
 real, intent(in):: log_p(kd)
 real, intent(in):: wz(is:ie,js:je,km+1)
 real, intent(in):: peln(is:ie,km+1,js:je)
 real, intent(out):: a2(is:ie,js:je,kd)
! local:
 real, dimension(2*km+1):: pn, gz
 integer n,i,j,k, k1, k2, l

 k2 = max(12, km/2+1)

!$OMP parallel do default(none) shared(k2,is,ie,js,je,km,kd,id,log_p,peln,a2,wz)   &
!$OMP             private(i,j,n,k,k1,l,pn,gz)
 do j=js,je
    do i=is,ie
!---------------
! Mirror method:
!---------------
       do k=1,km+1
          pn(k) = peln(i,k,j)
          gz(k) = wz(i,j,k)
       enddo
       do k=km+2, km+k2
          l = 2*(km+1) - k
          gz(k) = 2.*gz(km+1) - gz(l)
          pn(k) = 2.*pn(km+1) - pn(l)
       enddo
       k1 = 1
       do 1000 n=1,kd
          if( id(n)<0 ) goto 1000
          do k=k1,km+k2-1
             if( log_p(n) <= pn(k+1) .and. log_p(n) >= pn(k) ) then
                 a2(i,j,n) = gz(k) + (gz(k+1)-gz(k))*(log_p(n)-pn(k))/(pn(k+1)-pn(k))
                 k1 = k
                 go to 1000
             endif
          enddo
1000   continue
    enddo
 enddo

 end subroutine get_height_given_pressure

 subroutine prt_height(qname, is, ie, js, je, ng, km, press, phis, delz, peln, area, lat)
 character(len=*), intent(in)::  qname
 integer,  intent(in):: is, ie, js, je, ng, km
 real, intent(in):: press
 real, intent(in):: peln(is:ie,km+1,js:je)
 real, intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
 real, intent(in):: delz(is-ng:ie+ng,js-ng:je+ng,km)
 real(kind=R_GRID), intent(in), dimension(is:ie, js:je):: area, lat
! local:
 real:: a2(is:ie,js:je)
 real(kind=R_GRID), dimension(2*km+1):: pn, gz
 real(kind=R_GRID):: log_p
 integer i,j,k, k2, l

 log_p = log(press)
 k2 = max(12, km/2+1)

!$OMP parallel do default(none) shared(k2,is,ie,js,je,km,log_p,peln,phis,delz,a2)   &
!$OMP             private(i,j,k,l,pn,gz)
 do j=js,je
    do 1000 i=is,ie
!---------------
! Mirror method:
!---------------
       do k=1,km+1
          pn(k) = peln(i,k,j)
       enddo
       gz(km+1) = phis(i,j)/grav
       do k=km,1,-1
          gz(k) = gz(k+1) - delz(i,j,k)
       enddo
       do k=km+2, km+k2
          l = 2*(km+1) - k
          gz(k) = 2.*gz(km+1) - gz(l)
          pn(k) = 2.*pn(km+1) - pn(l)
       enddo

       do k=1,km+k2-1
          if( log_p <= pn(k+1) .and. log_p >= pn(k) ) then
              a2(i,j) = gz(k) + (gz(k+1)-gz(k))*(log_p-pn(k))/(pn(k+1)-pn(k))
              go to 1000
          endif
       enddo
1000   continue
 enddo
 call prt_gb_nh_sh(qname, is,ie, js,je, a2, area, lat)

 end subroutine prt_height

 subroutine prt_gb_nh_sh(qname, is,ie, js,je, a2, area, lat)
  character(len=*), intent(in)::  qname
  integer, intent(in):: is, ie, js, je
  real, intent(in), dimension(is:ie, js:je):: a2
  real(kind=R_GRID), intent(in), dimension(is:ie, js:je):: area, lat
! Local:
  real(R_GRID), parameter:: rad2deg = 180./pi
  real(R_GRID):: slat
  real:: t_eq, t_nh, t_sh, t_gb
  real:: area_eq, area_nh, area_sh, area_gb
  integer:: i,j

     t_eq = 0.   ;    t_nh = 0.;    t_sh = 0.;    t_gb = 0.
     area_eq = 0.; area_nh = 0.; area_sh = 0.; area_gb = 0.
     do j=js,je
        do i=is,ie
           slat = lat(i,j)*rad2deg
           area_gb = area_gb + area(i,j)
           t_gb = t_gb + a2(i,j)*area(i,j)
           if( (slat>-20. .and. slat<20.) ) then
                area_eq = area_eq + area(i,j)
                t_eq = t_eq + a2(i,j)*area(i,j)
           elseif( slat>=20. .and. slat<80. ) then
                area_nh = area_nh + area(i,j)
                t_nh = t_nh + a2(i,j)*area(i,j)
           elseif( slat<=-20. .and. slat>-80. ) then
                area_sh = area_sh + area(i,j)
                t_sh = t_sh + a2(i,j)*area(i,j)
           endif
        enddo
     enddo
     call mp_reduce_sum(area_gb)
     call mp_reduce_sum(   t_gb)
     call mp_reduce_sum(area_nh)
     call mp_reduce_sum(   t_nh)
     call mp_reduce_sum(area_sh)
     call mp_reduce_sum(   t_sh)
     call mp_reduce_sum(area_eq)
     call mp_reduce_sum(   t_eq)
     !Bugfix for non-global domains
     if (area_gb <= 1.) area_gb = -1.0
     if (area_nh <= 1.) area_nh = -1.0
     if (area_sh <= 1.) area_sh = -1.0
     if (area_eq <= 1.) area_eq = -1.0
     if (is_master()) write(*,*) qname, t_gb/area_gb, t_nh/area_nh, t_sh/area_sh, t_eq/area_eq

 end subroutine prt_gb_nh_sh

 subroutine cs3_interpolator(is, ie, js, je, km, qin, kd, pout, wz, pe, id, qout, iv)
! iv =-1: winds
! iv = 0: positive definite scalars
! iv = 1: temperature
 integer,  intent(in):: is, ie, js, je, km, iv
 integer,  intent(in):: kd
 integer,  intent(in):: id(kd)
 real, intent(in):: pout(kd)    ! must be monotonically increasing with increasing k
 real, intent(in):: pe(is:ie,km+1,js:je)
 real, intent(in):: qin(is:ie,js:je,km)
 real, intent(in):: wz(is:ie,js:je,km+1)
 real, intent(out):: qout(is:ie,js:je,kd)
! local:
 real, parameter:: gcp = grav / cp_air
 real:: qe(is:ie,km+1)
 real, dimension(is:ie,km):: q2, dp
 real:: s0, a6
 integer:: i,j,k, n, k1

!$OMP parallel do default(none) shared(iv,id,is,ie,js,je,km,kd,pout,qin,qout,pe,wz) & 
!$OMP             private(k1,s0,a6,q2,dp,qe)
 do j=js,je

   do i=is,ie
      do k=1,km
         dp(i,k) = pe(i,k+1,j) - pe(i,k,j)
         q2(i,k) = qin(i,j,k)
      enddo
   enddo

   call cs_prof(q2, dp, qe, km, is, ie, iv)

   do i=is,ie
     k1 = 1
     do n=1,kd
        if ( id(n) > 0 ) then
          if( pout(n) <= pe(i,1,j) ) then
! Higher than the top:
              qout(i,j,n) = qe(i,1)
          elseif ( pout(n) >= pe(i,km+1,j) ) then
! lower than the bottom surface:
            if ( iv==1 ) then    ! Temperature
! lower than the bottom surface:
! mean (hydro) potential temp based on lowest 2-3 layers (NCEP method)
! temp = ptm * p**cappa = ptm * exp(cappa*log(pout))
               qout(i,j,n) = gcp*exp(kappa*pout(n)) * (wz(i,j,km-2) - wz(i,j,km))  /   &
                               ( exp(kappa*pe(i,km,j)) - exp(kappa*pe(i,km-2,j)) )
            else
               qout(i,j,n) = qe(i,km+1)
            endif
          else 
            do k=k1,km
               if ( pout(n)>=pe(i,k,j) .and. pout(n) <= pe(i,k+1,j) ) then
! PPM distribution: f(s) = AL + s*[(AR-AL) + A6*(1-s)]         ( 0 <= s <= 1 )
                   a6 = 3.*(2.*q2(i,k) - (qe(i,k)+qe(i,k+1)))
                   s0 = (pout(n)-pe(i,k,j)) / dp(i,k)
                   qout(i,j,n) = qe(i,k) + s0*(qe(i,k+1)-qe(i,k)+a6*(1.-s0))
                   k1 = k     ! next level
                   go to 500
               endif
            enddo
          endif
500       if ( iv==0 ) qout(i,j,n) = max(0., qout(i,j,n))
        endif
     enddo
   enddo
 enddo

! Send_data here

 end subroutine cs3_interpolator

 subroutine cs_interpolator(is, ie, js, je, km, qin, zout, wz, qout, qmin)
 integer,  intent(in):: is, ie, js, je, km
 real, intent(in):: zout, qmin
 real, intent(in):: qin(is:ie,js:je,km)
 real, intent(in):: wz(is:ie,js:je,km+1)
 real, intent(out):: qout(is:ie,js:je)
! local:
 real:: qe(is:ie,km+1)
 real, dimension(is:ie,km):: q2, dz
 real:: s0, a6
 integer:: i,j,k

!$OMP parallel do default(none) shared(qmin,is,ie,js,je,km,zout,qin,qout,wz) & 
!$OMP             private(s0,a6,q2,dz,qe)
 do j=js,je

   do i=is,ie
      do k=1,km
         dz(i,k) = wz(i,j,k) - wz(i,j,k+1)
         q2(i,k) = qin(i,j,k)
      enddo
   enddo

   call cs_prof(q2, dz, qe, km, is, ie, 1)

   do i=is,ie
      if( zout >= wz(i,j,1) ) then
! Higher than the top:
          qout(i,j) = qe(i,1)
      elseif ( zout <= wz(i,j,km+1) ) then
          qout(i,j) = qe(i,km+1)
      else 
          do k=1,km
             if ( zout<=wz(i,j,k) .and. zout >= wz(i,j,k+1) ) then
! PPM distribution: f(s) = AL + s*[(AR-AL) + A6*(1-s)]         ( 0 <= s <= 1 )
                  a6 = 3.*(2.*q2(i,k) - (qe(i,k)+qe(i,k+1)))
                  s0 = (wz(i,j,k)-zout) / dz(i,k)
                  qout(i,j) = qe(i,k) + s0*(qe(i,k+1)-qe(i,k)+a6*(1.-s0))
                  go to 500
             endif
          enddo
      endif
500   qout(i,j) = max(qmin, qout(i,j))
   enddo
 enddo

! Send_data here

 end subroutine cs_interpolator

 subroutine cs_prof(q2, delp, q, km, i1, i2, iv)
! Latest: Dec 2015 S.-J. Lin, NOAA/GFDL
 integer, intent(in):: i1, i2, km
 integer, intent(in):: iv
 real, intent(in)   :: q2(i1:i2,km)
 real, intent(in)   :: delp(i1:i2,km)     ! <layer pressure thickness
 real, intent(out):: q(i1:i2,km+1)
!-----------------------------------------------------------------------
 real  gam(i1:i2,km)
 real   d4(i1:i2)
 real   bet, a_bot, grat
 integer i, k

  do i=i1,i2
         grat = delp(i,2) / delp(i,1)   ! grid ratio
          bet = grat*(grat+0.5)
       q(i,1) = ( (grat+grat)*(grat+1.)*q2(i,1) + q2(i,2) ) / bet
     gam(i,1) = ( 1. + grat*(grat+1.5) ) / bet
  enddo

  do k=2,km
     do i=i1,i2
           d4(i) = delp(i,k-1) / delp(i,k)
             bet =  2. + d4(i) + d4(i) - gam(i,k-1)
          q(i,k) = ( 3.*(q2(i,k-1)+d4(i)*q2(i,k)) - q(i,k-1) )/bet
        gam(i,k) = d4(i) / bet
     enddo
  enddo
 
  do i=i1,i2
         a_bot = 1. + d4(i)*(d4(i)+1.5)
     q(i,km+1) = (2.*d4(i)*(d4(i)+1.)*q2(i,km)+q2(i,km-1)-a_bot*q(i,km))  &
               / ( d4(i)*(d4(i)+0.5) - a_bot*gam(i,km) )
  enddo

  do k=km,1,-1
     do i=i1,i2
        q(i,k) = q(i,k) - gam(i,k)*q(i,k+1)
     enddo
  enddo

! Apply *large-scale* constraints 
  do i=i1,i2
     q(i,2) = min( q(i,2), max(q2(i,1), q2(i,2)) )
     q(i,2) = max( q(i,2), min(q2(i,1), q2(i,2)) )
  enddo

  do k=2,km
     do i=i1,i2
        gam(i,k) = q2(i,k) - q2(i,k-1)
     enddo
  enddo

! Interior:
  do k=3,km-1
     do i=i1,i2
        if ( gam(i,k-1)*gam(i,k+1)>0. ) then
! Apply large-scale constraint to ALL fields if not local max/min
             q(i,k) = min( q(i,k), max(q2(i,k-1),q2(i,k)) )
             q(i,k) = max( q(i,k), min(q2(i,k-1),q2(i,k)) )
        else
          if ( gam(i,k-1) > 0. ) then
! There exists a local max
               q(i,k) = max(q(i,k), min(q2(i,k-1),q2(i,k)))
          else
! There exists a local min
               q(i,k) = min(q(i,k), max(q2(i,k-1),q2(i,k)))
               if ( iv==0 ) q(i,k) = max(0., q(i,k))
          endif
        endif
     enddo
  enddo

! Bottom:
  do i=i1,i2
     q(i,km) = min( q(i,km), max(q2(i,km-1), q2(i,km)) )
     q(i,km) = max( q(i,km), min(q2(i,km-1), q2(i,km)) )
  enddo
  
 end subroutine cs_prof


 subroutine interpolate_vertical(is, ie, js, je, km, plev, peln, a3, a2)

 integer,  intent(in):: is, ie, js, je, km
 real, intent(in):: peln(is:ie,km+1,js:je)
 real, intent(in):: a3(is:ie,js:je,km)
 real, intent(in):: plev
 real, intent(out):: a2(is:ie,js:je)
! local:
 real pm(km)
 real logp
 integer i,j,k

 logp = log(plev)

!$OMP parallel do default(none) shared(is,ie,js,je,km,peln,logp,a2,a3) & 
!$OMP                          private(pm)
 do j=js,je
    do 1000 i=is,ie

       do k=1,km
          pm(k) = 0.5*(peln(i,k,j)+peln(i,k+1,j))
       enddo

       if( logp <= pm(1) ) then
           a2(i,j) = a3(i,j,1)
       elseif ( logp >= pm(km) ) then
           a2(i,j) = a3(i,j,km)
       else 
           do k=1,km-1
              if( logp <= pm(k+1) .and. logp >= pm(k) ) then
                  a2(i,j) = a3(i,j,k) + (a3(i,j,k+1)-a3(i,j,k))*(logp-pm(k))/(pm(k+1)-pm(k))
                  go to 1000
              endif
           enddo
       endif
1000   continue
 enddo

 end subroutine interpolate_vertical


 subroutine interpolate_z(is, ie, js, je, km, zl, hght, a3, a2)

 integer,  intent(in):: is, ie, js, je, km
 real, intent(in):: hght(is:ie,js:je,km+1)
 real, intent(in):: a3(is:ie,js:je,km)
 real, intent(in):: zl
 real, intent(out):: a2(is:ie,js:je)
! local:
 real zm(km)
 integer i,j,k


!$OMP parallel do default(none) shared(is,ie,js,je,km,hght,zl,a2,a3) private(zm)
 do j=js,je
    do 1000 i=is,ie
       do k=1,km
          zm(k) = 0.5*(hght(i,j,k)+hght(i,j,k+1))
       enddo
       if( zl >= zm(1) ) then
           a2(i,j) = a3(i,j,1)
       elseif ( zl <= zm(km) ) then
           a2(i,j) = a3(i,j,km)
       else 
           do k=1,km-1
              if( zl <= zm(k) .and. zl >= zm(k+1) ) then
                  a2(i,j) = a3(i,j,k) + (a3(i,j,k+1)-a3(i,j,k))*(zm(k)-zl)/(zm(k)-zm(k+1))
                  go to 1000
              endif
           enddo
       endif
1000   continue
 enddo

 end subroutine interpolate_z

 subroutine helicity_relative(is, ie, js, je, ng, km, zvir, sphum, srh,   &
                              ua, va, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
! !INPUT PARAMETERS:
   integer, intent(in):: is, ie, js, je, ng, km, sphum
   real, intent(in):: grav, zvir, z_bot, z_top
   real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, ua, va
   real, intent(in):: delz(is-ng:ie+ng,js-ng:je+ng,km)
   real, intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
   real, intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
   real, intent(in):: peln(is:ie,km+1,js:je) 
   logical, intent(in):: hydrostatic
   real, intent(out):: srh(is:ie,js:je)   ! unit: (m/s)**2
!   real, parameter:: z_crit = 3.e3   ! lowest 3-km
!---------------------------------------------------------------------------------
! SRH = 150-299 ... supercells possible with weak tornadoes
! SRH = 300-449 ... very favourable to supercells development and strong tornadoes
! SRH > 450 ... violent tornadoes
!---------------------------------------------------------------------------------
! if z_top = 1E3, the threshold for supercells is 100 (m/s)**2
! Coded by S.-J. Lin for CONUS regional climate simulations
!
   real:: rdg
   real, dimension(is:ie):: zh, uc, vc, dz, zh0
   integer i, j, k, k0, k1
   logical below(is:ie)

   rdg = rdgas / grav

!$OMP parallel do default(none) shared(is,ie,js,je,km,hydrostatic,rdg,pt,zvir,sphum, &
#ifdef MULTI_GASES
!$OMP                                  num_gas,                         &
#endif
!$OMP                                  peln,delz,ua,va,srh,z_bot,z_top) &
!$OMP                          private(zh,uc,vc,dz,k0,k1,zh0,below)
   do j=js,je

      do i=is,ie
         uc(i) = 0.
         vc(i) = 0.
         zh(i) = 0.
         srh(i,j) = 0.
         below(i) = .true.
         zh0(i) = 0.

!        if ( phis(i,j)/grav < 1.E3 ) then
         do k=km,1,-1
            if ( hydrostatic ) then
#ifdef MULTI_GASES
                 dz(i) = rdg*pt(i,j,k)*virq(q(i,j,k,1:num_gas))*(peln(i,k+1,j)-peln(i,k,j))
#else
                 dz(i) = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
#endif
            else
                 dz(i) = - delz(i,j,k)
            endif
            zh(i) = zh(i) + dz(i)
            if (zh(i) <= z_bot ) continue
            if (zh(i) > z_bot .and. below(i)) then
               uc(i) = ua(i,j,k)*dz(i)
               vc(i) = va(i,j,k)*dz(i)
               zh0(i) = zh(i) - dz(i)
               k1 = k
               below(i) = .false.
! Compute mean winds below z_top
            elseif ( zh(i) < z_top ) then
                uc(i) = uc(i) + ua(i,j,k)*dz(i)
                vc(i) = vc(i) + va(i,j,k)*dz(i)
                k0 = k
            else
                uc(i) = uc(i) / (zh(i)-dz(i) - zh0(i))
                vc(i) = vc(i) / (zh(i)-dz(i) - zh0(i))
                goto 123
            endif
         enddo 
123      continue

! Lowest layer wind shear computed betw top edge and mid-layer
         k = k1
         srh(i,j) = 0.5*(va(i,j,k1)-vc(i))*(ua(i,j,k1-1)-ua(i,j,k1))  -  &
                    0.5*(ua(i,j,k1)-uc(i))*(va(i,j,k1-1)-va(i,j,k1))
         do k=k0, k1-1
            srh(i,j) = srh(i,j) + 0.5*(va(i,j,k)-vc(i))*(ua(i,j,k-1)-ua(i,j,k+1)) -  &
                                  0.5*(ua(i,j,k)-uc(i))*(va(i,j,k-1)-va(i,j,k+1))
         enddo
!        endif
      enddo  ! i-loop
   enddo   ! j-loop

 end subroutine helicity_relative

 subroutine helicity_relative_caps(is, ie, js, je, ng, km, zvir, sphum, srh, uc, vc,  &
                              ua, va, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
! !INPUT PARAMETERS:
   integer, intent(in):: is, ie, js, je, ng, km, sphum
   real, intent(in):: grav, zvir, z_bot, z_top
   real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, ua, va
   real, intent(in):: delz(is-ng:ie+ng,js-ng:je+ng,km)
   real, intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
   real, intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
   real, intent(in):: peln(is:ie,km+1,js:je) 
   real, intent(in):: uc(is:ie,js:je), vc(is:ie,js:je)
   logical, intent(in):: hydrostatic
   real, intent(out):: srh(is:ie,js:je)   ! unit: (m/s)**2
!---------------------------------------------------------------------------------
! SRH = 150-299 ... supercells possible with weak tornadoes
! SRH = 300-449 ... very favourable to supercells development and strong tornadoes
! SRH > 450 ... violent tornadoes
!---------------------------------------------------------------------------------
! if z_crit = 1E3, the threshold for supercells is 100 (m/s)**2
! Coded by S.-J. Lin for CONUS regional climate simulations
!
   real:: rdg
   real, dimension(is:ie):: zh, dz, zh0
   integer i, j, k, k0, k1, n
   logical below

   rdg = rdgas / grav

!$OMP parallel do default(none) shared(is,ie,js,je,km,hydrostatic,rdg,pt,zvir,sphum, &
#ifdef MULTI_GASES
!$OMP                                  num_gas,                               &
#endif
!$OMP                                  peln,delz,ua,va,srh,uc,vc,z_bot,z_top) &
!$OMP                          private(zh,dz,k0,k1,zh0,below)
   do j=js,je

      do i=is,ie
         srh(i,j) = 0.
         zh(i) = 0.
         zh0 = 0.
         below = .true.

  k_loop:do k=km,1,-1
            if ( hydrostatic ) then
#ifdef MULTI_GASES
                 dz(i) = rdg*pt(i,j,k)*virq(q(i,j,k,1:num_gas))*(peln(i,k+1,j)-peln(i,k,j))
#else
                 dz(i) = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
#endif
            else
                 dz(i) = -delz(i,j,k)
            endif

            zh(i) = zh(i) + dz(i)
            if (zh(i) <= z_bot ) continue
            if (zh(i) > z_bot .and. below) then
               zh0(i) = zh(i) - dz(i)
               k1 = k
               below = .false.
! Compute mean winds below z_top
            elseif ( zh(i) < z_top ) then
                k0 = k
            else
               EXIT k_loop
            endif

         enddo k_loop

! Lowest layer wind shear computed betw top edge and mid-layer
         k = k1
         srh(i,j) = 0.5*(va(i,j,k1)-vc(i,j))*(ua(i,j,k1-1)-ua(i,j,k1))  -  &
                    0.5*(ua(i,j,k1)-uc(i,j))*(va(i,j,k1-1)-va(i,j,k1))
         do k=k0, k1-1
            srh(i,j) = srh(i,j) + 0.5*(va(i,j,k)-vc(i,j))*(ua(i,j,k-1)-ua(i,j,k+1)) -  &
                                  0.5*(ua(i,j,k)-uc(i,j))*(va(i,j,k-1)-va(i,j,k+1))
         enddo
      enddo  ! i-loop
   enddo   ! j-loop

 end subroutine helicity_relative_caps


 subroutine bunkers_vector(is, ie, js, je, ng, km, zvir, sphum, uc, vc,  &
                           ua, va, delz, q, hydrostatic, pt, peln, phis, grav)

   integer, intent(in):: is, ie, js, je, ng, km, sphum
   real, intent(in):: grav, zvir
   real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, ua, va
   real, intent(in):: delz(is-ng:ie+ng,js-ng:je+ng,km)
   real, intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
   real, intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
   real, intent(in):: peln(is:ie,km+1,js:je) 
   logical, intent(in):: hydrostatic
   real, intent(out):: uc(is:ie,js:je), vc(is:ie,js:je)

   real:: rdg
   real :: zh, dz, usfc, vsfc, u6km, v6km, umn, vmn
   real :: ushr, vshr, shrmag
   integer i, j, k, n
   real, parameter :: bunkers_d = 7.5 ! Empirically derived parameter
   logical :: has_sfc, has_6km

   rdg = rdgas / grav

!$OMP parallel do default(none) shared(is,ie,js,je,ng,km,hydrostatic,rdg,pt,zvir,sphum, &
#ifdef MULTI_GASES
!$OMP                                  num_gas,               &
#endif
!$OMP                                  peln,delz,ua,va,uc,vc) &
!$OMP                           private(zh,dz,usfc,vsfc,u6km,v6km,umn,vmn, &
!$OMP                                  ushr,vshr,shrmag)
   do j=js,je
      do i=is,ie
         zh = 0.
         usfc = 0.
         vsfc = 0.
         u6km = 0.
         v6km = 0.
         umn = 0.
         vmn = 0.

         usfc = ua(i,j,km)
         vsfc = va(i,j,km)

  k_loop:do k=km,1,-1
            if ( hydrostatic ) then
#ifdef MULTI_GASES
                 dz = rdg*pt(i,j,k)*virq(q(i,j,k,1:num_gas))*(peln(i,k+1,j)-peln(i,k,j))
#else
                 dz = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
#endif
            else
                 dz = -delz(i,j,k)
            endif
            zh = zh + dz

            if (zh < 6000) then
                u6km = ua(i,j,k)
                v6km = va(i,j,k)

                umn = umn + ua(i,j,k)*dz
                vmn = vmn + va(i,j,k)*dz
            else
               EXIT k_loop
            endif

         enddo k_loop

         u6km = u6km + (ua(i,j,k) - u6km) / dz * (6000. - (zh - dz))
         v6km = v6km + (va(i,j,k) - v6km) / dz * (6000. - (zh - dz))

         umn = umn / (zh - dz)
         vmn = vmn / (zh - dz)

         ushr = u6km - usfc
         vshr = v6km - vsfc
         shrmag = sqrt(ushr * ushr + vshr * vshr)
         uc(i,j) = umn + bunkers_d * vshr / shrmag
         vc(i,j) = vmn - bunkers_d * ushr / shrmag

      enddo  ! i-loop
   enddo   ! j-loop

 end subroutine bunkers_vector


 subroutine updraft_helicity(is, ie, js, je, ng, km, zvir, sphum, uh,   &
                             w, vort, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
! !INPUT PARAMETERS:
   integer, intent(in):: is, ie, js, je, ng, km, sphum
   real, intent(in):: grav, zvir, z_bot, z_top
   real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, w
   real, intent(in), dimension(is:ie,js:je,km):: vort
   real, intent(in):: delz(is-ng:ie+ng,js-ng:je+ng,km)
   real, intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
   real, intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
   real, intent(in):: peln(is:ie,km+1,js:je) 
   logical, intent(in):: hydrostatic
   real, intent(out):: uh(is:ie,js:je)   ! unit: (m/s)**2
! Coded by S.-J. Lin for CONUS regional climate simulations
! Modified for UH by LMH
!
   real:: rdg
   real, dimension(is:ie):: zh, dz, zh0
   integer i, j, k, n
   logical below(is:ie)

   rdg = rdgas / grav

!$OMP parallel do default(none) shared(is,ie,js,je,ng,km,hydrostatic,rdg,pt,zvir,sphum, &
#ifdef MULTI_GASES
!$OMP                                  num_gas,                         &
#endif
!$OMP                                  peln,delz,w,vort,uh,z_bot,z_top) &
!$OMP                          private(zh,dz,zh0,below)
   do j=js,je

      do i=is,ie
         zh(i) = 0.
         uh(i,j) = 0.
         below(i) = .true.
         zh0(i) = 0.

!        if ( phis(i,j)/grav < 1.E3 ) then
         do k=km,1,-1
            if ( hydrostatic ) then
#ifdef MULTI_GASES
                 dz(i) = rdg*pt(i,j,k)*virq(q(i,j,k,1:num_gas))*(peln(i,k+1,j)-peln(i,k,j))
#else
                 dz(i) = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
#endif
            else
                 dz(i) = - delz(i,j,k)
            endif
            zh(i) = zh(i) + dz(i)
            if (zh(i) <= z_bot ) continue
            if (zh(i) > z_bot .and. below(i)) then
               uh(i,j) = vort(i,j,k)*w(i,j,k)*(zh(i) - z_bot)
               below(i) = .false.
! Compute mean winds below z_top
            elseif ( zh(i) < z_top ) then
               uh(i,j) = uh(i,j) + vort(i,j,k)*w(i,j,k)*dz(i)
            else
               uh(i,j) = uh(i,j) + vort(i,j,k)*w(i,j,k)*(z_top - (zh(i)-dz(i)) ) 
               goto 123
            endif
         enddo
123      continue

      enddo  ! i-loop
   enddo   ! j-loop

 end subroutine updraft_helicity

 subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav)

! !INPUT PARAMETERS:
   integer, intent(in)::  is, ie, js, je, ng, km
   real, intent(in):: grav
   real, intent(in):: pt(is-ng:ie+ng,js-ng:je+ng,km) 
   real, intent(in):: pkz(is:ie,js:je,km) 
   real, intent(in):: delp(is-ng:ie+ng,js-ng:je+ng,km)
   real, intent(in):: f_d(is-ng:ie+ng,js-ng:je+ng) 

! vort is relative vorticity as input. Becomes PV on output
      real, intent(inout):: vort(is:ie,js:je,km)

! !DESCRIPTION:
!        EPV = 1/r * (vort+f_d) * d(S)/dz; where S is a conservative scalar
!        r the fluid density, and S is chosen to be the entropy here: S = log(pt)
!        pt == potential temperature.
! Computation are performed assuming the input is on "x-y-z" Cartesian coordinate.
! The approximation made here is that the relative vorticity computed on constant
! z-surface is not that different from the hybrid sigma-p coordinate.
! See page 39, Pedlosky 1979: Geophysical Fluid Dynamics
!
! The follwoing simplified form is strictly correct only if vort is computed on 
! constant z surfaces. In addition hydrostatic approximation is made.
!        EPV = - GRAV * (vort+f_d) / del(p) * del(pt) / pt 
! where del() is the vertical difference operator.
!
! programmer: S.-J. Lin, shian-jiann.lin@noaa.gov
!
!EOP
!---------------------------------------------------------------------
!BOC
      real w3d(is:ie,js:je,km)
      real te(is:ie,js:je,km+1), t2(is:ie,km), delp2(is:ie,km)
      real te2(is:ie,km+1)
      integer i, j, k

#ifdef SW_DYNAMICS
!$OMP parallel do default(none) shared(is,ie,js,je,vort,grav,f_d,delp)
        do j=js,je
          do i=is,ie
            vort(i,j,1) = grav * (vort(i,j,1)+f_d(i,j)) / delp(i,j,1)
          enddo
        enddo
#else
! Compute PT at layer edges.
!$OMP parallel do default(none) shared(is,ie,js,je,km,pt,pkz,w3d,delp,te2,te) &
!$OMP                          private(t2, delp2) 
     do j=js,je
        do k=1,km
          do i=is,ie
               t2(i,k) = pt(i,j,k) / pkz(i,j,k)
              w3d(i,j,k) = t2(i,k)
            delp2(i,k) = delp(i,j,k)
          enddo
        enddo

        call ppme(t2, te2, delp2, ie-is+1, km)

        do k=1,km+1
           do i=is,ie
              te(i,j,k) = te2(i,k)
           enddo
        enddo
     enddo

!$OMP parallel do default(none) shared(is,ie,js,je,km,vort,f_d,te,w3d,delp,grav)
     do k=1,km
        do j=js,je
          do i=is,ie
! Entropy is the thermodynamic variable in the following form
            vort(i,j,k) = (vort(i,j,k)+f_d(i,j)) * ( te(i,j,k)-te(i,j,k+1) )  &
                          / ( w3d(i,j,k)*delp(i,j,k) )  * grav
          enddo
        enddo
     enddo
#endif

 end subroutine pv_entropy


 subroutine ppme(p,qe,delp,im,km)

  integer, intent(in):: im, km
  real, intent(in)::    p(im,km)
  real, intent(in):: delp(im,km)
  real, intent(out)::qe(im,km+1)

! local arrays.
      real dc(im,km),delq(im,km), a6(im,km)
      real c1, c2, c3, tmp, qmax, qmin
      real a1, a2, s1, s2, s3, s4, ss3, s32, s34, s42
      real a3, b2, sc, dm, d1, d2, f1, f2, f3, f4
      real qm, dq
      integer i, k, km1

      km1 = km - 1

      do 500 k=2,km
      do 500 i=1,im
500   a6(i,k) = delp(i,k-1) + delp(i,k)

      do 1000 k=1,km1
      do 1000 i=1,im
      delq(i,k) = p(i,k+1) - p(i,k)
1000  continue

      do 1220 k=2,km1
      do 1220 i=1,im
      c1 = (delp(i,k-1)+0.5*delp(i,k))/a6(i,k+1)
      c2 = (delp(i,k+1)+0.5*delp(i,k))/a6(i,k)
      tmp = delp(i,k)*(c1*delq(i,k) + c2*delq(i,k-1)) /    &
                                    (a6(i,k)+delp(i,k+1))
      qmax = max(p(i,k-1),p(i,k),p(i,k+1)) - p(i,k)
      qmin = p(i,k) - min(p(i,k-1),p(i,k),p(i,k+1))
      dc(i,k) = sign(min(abs(tmp),qmax,qmin), tmp)
1220  continue

!****6***0*********0*********0*********0*********0*********0**********72
! 4th order interpolation of the provisional cell edge value
!****6***0*********0*********0*********0*********0*********0**********72

   do k=3,km1
      do i=1,im
         c1 = delq(i,k-1)*delp(i,k-1) / a6(i,k)
         a1 = a6(i,k-1) / (a6(i,k) + delp(i,k-1))
         a2 = a6(i,k+1) / (a6(i,k) + delp(i,k))
         qe(i,k) = p(i,k-1) + c1 + 2./(a6(i,k-1)+a6(i,k+1)) *        &
                   ( delp(i,k)*(c1*(a1 - a2)+a2*dc(i,k-1)) -         &
                                delp(i,k-1)*a1*dc(i,k  ) )
      enddo
   enddo

! three-cell parabolic subgrid distribution at model top

   do i=1,im
! three-cell PP-distribution
! Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp
! a3 = a / 3
! b2 = b / 2
      s1 = delp(i,1)
      s2 = delp(i,2) + s1
!
      s3 = delp(i,2) + delp(i,3)
      s4 = s3 + delp(i,4)
      ss3 =  s3 + s1
      s32 = s3*s3
      s42 = s4*s4
      s34 = s3*s4
! model top
      a3 = (delq(i,2) - delq(i,1)*s3/s2) / (s3*ss3)
!
      if(abs(a3) .gt. 1.e-14) then
         b2 =  delq(i,1)/s2 - a3*(s1+s2)
         sc = -b2/(3.*a3)
         if(sc .lt. 0. .or. sc .gt. s1) then
             qe(i,1) = p(i,1) - s1*(a3*s1 + b2)
         else
             qe(i,1) = p(i,1) - delq(i,1)*s1/s2
         endif
      else
! Linear
         qe(i,1) = p(i,1) - delq(i,1)*s1/s2
      endif
      dc(i,1) = p(i,1) - qe(i,1)
! compute coef. for the off-centered area preserving cubic poly.
      dm = delp(i,1) / (s34*ss3*(delp(i,2)+s3)*(s4+delp(i,1)))
      f1 = delp(i,2)*s34 / ( s2*ss3*(s4+delp(i,1)) )
      f2 = (delp(i,2)+s3) * (ss3*(delp(i,2)*s3+s34+delp(i,2)*s4)   &
            + s42*(delp(i,2)+s3+s32/s2))
      f3 = -delp(i,2)*( ss3*(s32*(s3+s4)/(s4-delp(i,2))            &
            + (delp(i,2)*s3+s34+delp(i,2)*s4))                     &
            + s42*(delp(i,2)+s3) )
      f4 = ss3*delp(i,2)*s32*(delp(i,2)+s3) / (s4-delp(i,2))
      qe(i,2) = f1*p(i,1)+(f2*p(i,2)+f3*p(i,3)+f4*p(i,4))*dm
   enddo

! Bottom
! Area preserving cubic with 2nd deriv. = 0 at the surface
   do i=1,im
      d1 = delp(i,km)
      d2 = delp(i,km1)
      qm = (d2*p(i,km)+d1*p(i,km1)) / (d1+d2)
      dq = 2.*(p(i,km1)-p(i,km)) / (d1+d2)
      c1 = (qe(i,km1)-qm-d2*dq) / (d2*(2.*d2*d2+d1*(d2+3.*d1)))
      c3 = dq - 2.0*c1*(d2*(5.*d1+d2)-3.*d1**2)
      qe(i,km  ) = qm - c1*d1*d2*(d2+3.*d1)
      qe(i,km+1) = d1*(8.*c1*d1**2-c3) + qe(i,km)
   enddo

 end subroutine ppme

subroutine rh_calc (pfull, t, qv, rh, do_cmip)
  
   real, intent (in),  dimension(:,:) :: pfull, t, qv
   real, intent (out), dimension(:,:) :: rh
   real, dimension(size(t,1),size(t,2)) :: esat
   logical, intent(in), optional :: do_cmip

   real, parameter :: d622 = rdgas/rvgas
   real, parameter :: d378 = 1.-d622

   logical :: do_simple = .false.

! because Betts-Miller uses a simplified scheme for calculating the relative humidity

     if (do_simple) then
        call lookup_es (t, esat)
        rh(:,:) = pfull(:,:)
        rh(:,:) = max(rh(:,:),esat(:,:))  !limit where pfull ~ esat
        rh(:,:)=100.*qv(:,:)/(d622*esat(:,:)/rh(:,:))
     else
        if (present(do_cmip)) then
         call compute_qs (t, pfull, rh, q=qv, es_over_liq_and_ice = .true.)
         rh(:,:)=100.*qv(:,:)/rh(:,:)
        else
        call compute_qs (t, pfull, rh, q=qv)
        rh(:,:)=100.*qv(:,:)/rh(:,:)
     endif
     endif

end subroutine rh_calc

#ifdef SIMPLIFIED_THETA_E

#ifdef MULTI_GASES
subroutine eqv_pot(theta_e, pt, delp, delz, peln, pkz, qi, is, ie, js, je, ng, npz, &
                   hydrostatic, moist)
#else
subroutine eqv_pot(theta_e, pt, delp, delz, peln, pkz, q, is, ie, js, je, ng, npz, &
                   hydrostatic, moist)
#endif
    integer, intent(in):: is,ie,js,je,ng,npz
#ifdef MULTI_GASES
    real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,npz,*):: qi
#endif
    real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,npz):: pt, delp, q
    real, intent(in), dimension(is-ng:     ,js-ng:     ,1: ):: delz
    real, intent(in), dimension(is:ie,npz+1,js:je):: peln
    real, intent(in):: pkz(is:ie,js:je,npz) 
    logical, intent(in):: hydrostatic, moist
! Output:
    real, dimension(is:ie,js:je,npz), intent(out) :: theta_e
! local
    real, parameter:: tice = 273.16
    real, parameter:: c_liq = 4190.
#ifdef SIM_NGGPS
    real, parameter:: dc_vap = 0.
#else
    real, parameter:: dc_vap = cp_vapor - c_liq     
#endif
    real(kind=R_GRID), dimension(is:ie):: pd, rq
    real(kind=R_GRID) :: wfac
    integer :: i,j,k

    if ( moist ) then
         wfac = 1.
    else
         wfac = 0.
    endif
#ifdef MULLTI_GASES
    do j=js,je
       do i=is,ie
          do k=1,npz
             q(i,j,k) = qi(i,j,k,1)
          enddo
       enddo
    enddo
#endif

!$OMP parallel do default(none) shared(pk0,wfac,moist,pkz,is,ie,js,je,npz,pt,q,delp,peln,delz,theta_e,hydrostatic)  &
!$OMP  private(pd, rq)
    do k = 1,npz
       do j = js,je

        if ( hydrostatic ) then
            do i=is,ie
               rq(i) = max(0., wfac*q(i,j,k))
               pd(i) = (1.-rq(i))*delp(i,j,k) / (peln(i,k+1,j) - peln(i,k,j))
            enddo
        else
! Dry pressure: p = r R T
            do i=is,ie
               rq(i) = max(0., wfac*q(i,j,k))
               pd(i) = -rdgas*pt(i,j,k)*(1.-rq(i))*delp(i,j,k)/(grav*delz(i,j,k))
            enddo
        endif

        if ( moist ) then
            do i=is,ie
               rq(i) = max(0., q(i,j,k))
!              rh(i) = max(1.e-12, rq(i)/wqs1(pt(i,j,k),den(i)))   ! relative humidity
!              theta_e(i,j,k) = exp(rq(i)/cp_air*((hlv+dc_vap*(pt(i,j,k)-tice))/pt(i,j,k) -   &
!                                   rvgas*log(rh(i))) + kappa*log(1.e5/pd(i))) * pt(i,j,k)
! Simplified form: (ignoring the RH term)
#ifdef SIM_NGGPS
#ifdef MULTI_GASES
               theta_e(i,j,k) = pt(i,j,k)*exp(kappa * (virqd(qi(i,j,k,:))/vicpqd(qi(i,j,k,:)))*log(1.e5/pd(i))) *  &
                                          exp(rq(i)*hlv/(cp_air*vicpqd(qi(i,j,k,:))*pt(i,j,k)))
#else
               theta_e(i,j,k) = pt(i,j,k)*exp(kappa*log(1.e5/pd(i))) *  &
                                          exp(rq(i)*hlv/(cp_air*pt(i,j,k)))
#endif
#else
#ifdef MULTI_GASES
               theta_e(i,j,k) = pt(i,j,k)*exp( rq(i)/(cp_air*vicpqd(q(i,j,k,:))*pt(i,j,k))*(hlv+dc_vap*(pt(i,j,k)-tice)) &
                                             + rdgas*virqd(qi(i,j,k,:)) / (cp_air*vicpqd(qi(i,j,k,:)))*log(1.e5/pd(i)) )
#else
               theta_e(i,j,k) = pt(i,j,k)*exp( rq(i)/(cp_air*pt(i,j,k))*(hlv+dc_vap*(pt(i,j,k)-tice)) &
                                             + kappa*log(1.e5/pd(i)) )
#endif
#endif
            enddo
        else
          if ( hydrostatic ) then
             do i=is,ie
                theta_e(i,j,k) = pt(i,j,k)*pk0/pkz(i,j,k)
             enddo
          else
             do i=is,ie
!               theta_e(i,j,k) = pt(i,j,k)*(1.e5/pd(i))**kappa
#ifdef MULTI_GASES
                theta_e(i,j,k) = pt(i,j,k)*exp( kappa * (virqd(qi(i,j,k,:))/vicpqd(qi(i,j,k,:)))*log(1.e5/pd(i)) )
#else
                theta_e(i,j,k) = pt(i,j,k)*exp( kappa*log(1.e5/pd(i)) )
#endif
             enddo
          endif
        endif
      enddo ! j-loop
    enddo   ! k-loop

end subroutine eqv_pot

#else

#ifdef MULTI_GASES
subroutine eqv_pot(theta_e, pt, delp, delz, peln, pkz, qi, is, ie, js, je, ng, npz, &
#else
subroutine eqv_pot(theta_e, pt, delp, delz, peln, pkz, q, is, ie, js, je, ng, npz, &
#endif
                   hydrostatic, moist)
! Modified by SJL
    integer, intent(in):: is,ie,js,je,ng,npz
#ifdef MULTI_GASES
    real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,npz,*):: qi
#else
    real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,npz):: q
#endif
    real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,npz):: pt, delp
    real, intent(in), dimension(is-ng:     ,js-ng:     ,1: ):: delz
    real, intent(in), dimension(is:ie,npz+1,js:je):: peln
    real, intent(in):: pkz(is:ie,js:je,npz) 
    logical, intent(in):: hydrostatic, moist
! Output:
    real, dimension(is:ie,js:je,npz), intent(out) :: theta_e
! local
#ifdef MULTI_GASES
    real, dimension(is-ng:ie+ng,js-ng:je+ng,npz):: q
#endif
    real, parameter:: cv_vap = cp_vapor - rvgas  ! 1384.5
    real, parameter:: cappa_b = 0.2854
    real(kind=R_GRID):: cv_air, cappa, zvir
    real(kind=R_GRID):: p_mb(is:ie)
    real(kind=R_GRID) :: r, e, t_l, rdg, capa
    integer :: i,j,k, n

#ifdef MULTI_GASES
    q(:,:,:) = qi(:,:,:,1)
#endif
    cv_air =  cp_air - rdgas
    rdg = -rdgas/grav
    if ( moist ) then
         zvir = rvgas/rdgas - 1.
    else
         zvir = 0.
    endif

!$OMP parallel do default(none) shared(moist,pk0,pkz,cv_air,zvir,rdg,is,ie,js,je,ng,npz, &
#ifdef MULTI_GASES
!$OMP      qi,num_gas,                                                                   &
#endif
!$OMP      pt,q,delp,peln,delz,theta_e,hydrostatic)  &
!$OMP      private(cappa,p_mb, r, e, t_l, capa)
    do k = 1,npz
       cappa = cappa_b
      do j = js,je
! get pressure in mb
        if ( hydrostatic ) then
            do i=is,ie
               p_mb(i) = 0.01*delp(i,j,k) / (peln(i,k+1,j) - peln(i,k,j))
            enddo
        else
            do i=is,ie
#ifdef MULTI_GASES
               p_mb(i) = 0.01*rdg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)*virq(qi(i,j,k,1:num_gas))
#else
               p_mb(i) = 0.01*rdg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)*(1.+zvir*q(i,j,k))
#endif
            enddo
        endif
        if ( moist ) then
          do i = is,ie
#ifdef MULTI_GASES
          cappa = rdgas/(rdgas+((1.-q(i,j,k))*cv_air+q(i,j,k)*cv_vap)/virq(qi(i,j,k,1:num_gas)))
#else
          cappa = rdgas/(rdgas+((1.-q(i,j,k))*cv_air+q(i,j,k)*cv_vap)/(1.+zvir*q(i,j,k)))
#endif
! get "dry" mixing ratio of m_vapor/m_tot in g/kg
          r = q(i,j,k)/(1.-q(i,j,k))*1000.
          r = max(1.e-10, r)
! get water vapor pressure
          e = p_mb(i)*r/(622.+r)
! get temperature at the lifting condensation level
! eq. 21 of Bolton 1980
          t_l = 2840./(3.5*log(pt(i,j,k))-log(e)-4.805)+55.
! get the equivalent potential temperature
!         theta_e(i,j,k) = pt(i,j,k)*exp( (cappa*(1.-0.28e-3*r)*log(1000./p_mb(i))) * &
!                          exp( (3.376/t_l-0.00254)*r*(1.+0.81e-3*r) )
           capa = cappa*(1. - r*0.28e-3)
           theta_e(i,j,k) = exp( (3.376/t_l-0.00254)*r*(1.+r*0.81e-3) )*pt(i,j,k)*(1000./p_mb(i))**capa
          enddo
        else
          if ( hydrostatic ) then
             do i = is,ie
                theta_e(i,j,k) = pt(i,j,k)*pk0/pkz(i,j,k)
             enddo
          else
             do i = is,ie
#ifdef MULTI_GASES
                theta_e(i,j,k) = pt(i,j,k)*exp( kappa * virqd(qi(i,j,k,1:num_gas))/vicpqd(qi(i,j,k,1:num_gas)) *log(1000./p_mb(i)) )
#else
                theta_e(i,j,k) = pt(i,j,k)*exp( kappa*log(1000./p_mb(i)) )
#endif
             enddo
          endif
        endif
      enddo ! j-loop
    enddo   ! k-loop

end subroutine eqv_pot

#endif

 subroutine nh_total_energy(is, ie, js, je, isd, ied, jsd, jed, km,  &
                            w, delz, pt, delp, q, hs, area, domain,  &
                            sphum, liq_wat, rainwat, ice_wat,        &
                            snowwat, graupel, nwat, ua, va, moist_phys, te)
! INPUT PARAMETERS:
   integer,  intent(in):: km, is, ie, js, je, isd, ied, jsd, jed
   integer,  intent(in):: nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel
   real, intent(in), dimension(isd:ied,jsd:jed,km):: ua, va, pt, delp, w, delz
   real, intent(in), dimension(isd:ied,jsd:jed,km,nwat):: q
   real, intent(in):: hs(isd:ied,jsd:jed)
   real, intent(in):: area(isd:ied, jsd:jed)
   logical, intent(in):: moist_phys
   type(domain2d), intent(INOUT) :: domain
   real, intent(out):: te(is:ie,js:je)
! Local
   real, parameter:: c_liq = 4190.
   real(kind=R_Grid) ::    area_l(isd:ied, jsd:jed)
   real, parameter:: cv_vap = cp_vapor - rvgas  
   real  phiz(is:ie,km+1)
   real, dimension(is:ie):: cvm, qc
   real cv_air, psm
   integer i, j, k

   area_l = area
   cv_air =  cp_air - rdgas

!$OMP parallel do default(none) shared(te,nwat,is,ie,js,je,isd,ied,jsd,jed,km,ua,va,   &
#ifdef MULTI_GASES
!$OMP          num_gas,                                                                &
#endif
!$OMP          w,q,pt,delp,delz,hs,cv_air,moist_phys,sphum,liq_wat,rainwat,ice_wat,snowwat,graupel) &
!$OMP          private(phiz,cvm, qc)
  do j=js,je

     do i=is,ie
        te(i,j) = 0.
        phiz(i,km+1) = hs(i,j)
     enddo

     do i=is,ie
        do k=km,1,-1
           phiz(i,k) = phiz(i,k+1) - grav*delz(i,j,k)
        enddo
     enddo

     if ( moist_phys ) then
        do k=1,km
           call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat,    &
                         ice_wat, snowwat, graupel, q, qc, cvm)
           do i=is,ie
              te(i,j) = te(i,j) + delp(i,j,k)*( cvm(i)*pt(i,j,k) + hlv*q(i,j,k,sphum) +  &
                      0.5*(phiz(i,k)+phiz(i,k+1)+ua(i,j,k)**2+va(i,j,k)**2+w(i,j,k)**2) )
           enddo
        enddo
     else
       do k=1,km
          do i=is,ie
#ifdef MULTI_GASES
             te(i,j) = te(i,j) + delp(i,j,k)*( cv_air*vicvqd(q(i,j,k,1:num_gas))*pt(i,j,k) +  &
                     0.5*(phiz(i,k)+phiz(i,k+1)+ua(i,j,k)**2+va(i,j,k)**2+w(i,j,k)**2) )
#else
             te(i,j) = te(i,j) + delp(i,j,k)*( cv_air*pt(i,j,k) +  &
                     0.5*(phiz(i,k)+phiz(i,k+1)+ua(i,j,k)**2+va(i,j,k)**2+w(i,j,k)**2) )
#endif
          enddo
       enddo
     endif
! Unit: kg*(m/s)^2/m^2 = Joule/m^2
     do i=is,ie
        te(i,j) = te(i,j)/grav
     enddo
  enddo

  psm = g_sum(domain, te, is, ie, js, je, 3, area_l, 1) 
  if( master ) write(*,*) 'TE ( Joule/m^2 * E9) =',  psm * 1.e-9

  end subroutine nh_total_energy

 subroutine dbzcalc(q, pt, delp, peln, delz, &
      dbz, maxdbz, allmax, bd, npz, ncnst, &
      hydrostatic, zvir, in0r, in0s, in0g, iliqskin)

! Code from Mark Stoelinga's dbzcalc.f from the RIP package. 
! Currently just using values taken directly from that code, which is
! consistent for the MM5 Reisner-2 microphysics. From that file:

!     This routine computes equivalent reflectivity factor (in dBZ) at
!     each model grid point.  In calculating Ze, the RIP algorithm makes
!     assumptions consistent with those made in an early version
!     (ca. 1996) of the bulk mixed-phase microphysical scheme in the MM5
!     model (i.e., the scheme known as "Resiner-2").  For each species:
!
!     1. Particles are assumed to be spheres of constant density.  The
!     densities of rain drops, snow particles, and graupel particles are
!     taken to be rho_r = rho_l = 1000 kg m^-3, rho_s = 100 kg m^-3, and
!     rho_g = 400 kg m^-3, respectively. (l refers to the density of
!     liquid water.)
!
!     2. The size distribution (in terms of the actual diameter of the
!     particles, rather than the melted diameter or the equivalent solid
!     ice sphere diameter) is assumed to follow an exponential
!     distribution of the form N(D) = N_0 * exp( lambda*D ).
!
!     3. If in0X=0, the intercept parameter is assumed constant (as in
!     early Reisner-2), with values of 8x10^6, 2x10^7, and 4x10^6 m^-4,
!     for rain, snow, and graupel, respectively.  Various choices of
!     in0X are available (or can be added).  Currently, in0X=1 gives the
!     variable intercept for each species that is consistent with
!     Thompson, Rasmussen, and Manning (2004, Monthly Weather Review,
!     Vol. 132, No. 2, pp. 519-542.)
!
!     4. If iliqskin=1, frozen particles that are at a temperature above
!     freezing are assumed to scatter as a liquid particle.
!
!     More information on the derivation of simulated reflectivity in RIP
!     can be found in Stoelinga (2005, unpublished write-up).  Contact
!     Mark Stoelinga (stoeling@atmos.washington.edu) for a copy.  

! 22sep16: Modifying to use the GFDL MP parameters. If doing so remember
!   that the GFDL MP assumes a constant intercept (in0X = .false.)
!   Ferrier-Aligo has an option for fixed slope (rather than fixed intercept).
!   Thompson presumably is an extension of Reisner MP.

   type(fv_grid_bounds_type), intent(IN) :: bd
   integer, intent(IN) :: npz, ncnst
   real,    intent(IN),  dimension(bd%isd:bd%ied, bd%jsd:bd%jed, npz) :: pt, delp, delz
   real,    intent(IN),  dimension(bd%isd:bd%ied, bd%jsd:bd%jed, npz, ncnst) :: q
   real,    intent(IN),  dimension(bd%is :bd%ie,  npz+1, bd%js:bd%je) :: peln
   real,    intent(OUT), dimension(bd%is :bd%ie,  bd%js :bd%je , npz) :: dbz
   real,    intent(OUT), dimension(bd%is :bd%ie,  bd%js :bd%je)      :: maxdbz
   logical, intent(IN) :: hydrostatic, in0r, in0s, in0g, iliqskin
   real,    intent(IN) :: zvir
   real,    intent(OUT) :: allmax

   !Parameters for constant intercepts (in0[rsg] = .false.)
   !Using GFDL MP values
   real(kind=R_GRID), parameter:: rn0_r = 8.e6 ! m^-4
   real(kind=R_GRID), parameter:: rn0_s = 3.e6 ! m^-4
   real(kind=R_GRID), parameter:: rn0_g = 4.e6 ! m^-4
   real(kind=R_GRID), parameter:: vconr = 2503.23638966667
   real(kind=R_GRID), parameter:: vcong =   87.2382675
   real(kind=R_GRID), parameter:: vcons =    6.6280504
   real(kind=R_GRID), parameter:: normr = 25132741228.7183
   real(kind=R_GRID), parameter:: normg =  5026548245.74367
   real(kind=R_GRID), parameter:: norms =   942477796.076938

   !Constants for variable intercepts
   !Will need to be changed based on MP scheme
   real, parameter :: r1=1.e-15
   real, parameter :: ron=8.e6
   real, parameter :: ron2=1.e10
   real, parameter :: son=2.e7
   real, parameter :: gon=5.e7
   real, parameter :: ron_min = 8.e6
   real, parameter :: ron_qr0 = 0.00010
   real, parameter :: ron_delqr0 = 0.25*ron_qr0
   real, parameter :: ron_const1r = (ron2-ron_min)*0.5
   real, parameter :: ron_const2r = (ron2+ron_min)*0.5

   !Other constants
   real, parameter :: gamma_seven = 720.
   !The following values are also used in GFDL MP
   real, parameter :: rho_r  = 1.0e3  ! LFO83
   real, parameter :: rho_s  = 100.   ! kg m^-3 
   real, parameter :: rho_g0 = 400.   ! kg m^-3
   real, parameter :: rho_g  = 500.   ! graupel-hail mix
!  real, parameter :: rho_g  = 900.   ! hail/frozen rain
   real, parameter :: alpha  = 0.224
   real(kind=R_GRID), parameter :: factor_r = gamma_seven * 1.e18 * (1./(pi*rho_r))**1.75
   real(kind=R_GRID), parameter :: factor_s = gamma_seven * 1.e18 * (1./(pi*rho_s))**1.75 &
        * (rho_s/rho_r)**2 * alpha
   real(kind=R_GRID), parameter :: factor_g = gamma_seven * 1.e18 * (1./(pi*rho_g))**1.75 &
        * (rho_g/rho_r)**2 * alpha
   real, parameter :: qmin = 1.e-12
   real, parameter :: tice = 273.16

! Double precision
   real(kind=R_GRID):: rhoair(bd%is:bd%ie)
   real(kind=R_GRID):: qr1, qs1, qg1, t1, t2, t3, rwat, denfac, vtr, vtg, vts
   real(kind=R_GRID):: factorb_s, factorb_g
   real(kind=R_GRID):: temp_c, pres, sonv, gonv, ronv, z_e

   integer :: i,j,k
   integer :: is, ie, js, je

   is = bd%is
   ie = bd%ie
   js = bd%js
   je = bd%je
   if (rainwat < 1) return

   dbz(:,:,1:mp_top) = -20.
   maxdbz(:,:) = -20. !Minimum value
   allmax = -20.


!$OMP parallel do default(shared) private(rhoair,t1,t2,t3,denfac,vtr,vtg,vts,z_e)
   do k=mp_top+1, npz
   do j=js, je
      if (hydrostatic) then
         do i=is, ie
#ifdef MULTI_GASES
            rhoair(i) = delp(i,j,k)/( (peln(i,k+1,j)-peln(i,k,j)) * rdgas * pt(i,j,k) * virq(q(i,j,k,1:num_gas)) )
#else
            rhoair(i) = delp(i,j,k)/( (peln(i,k+1,j)-peln(i,k,j)) * rdgas * pt(i,j,k) * ( 1. + zvir*q(i,j,k,sphum) ) )
#endif
         enddo
      else
         do i=is, ie
            rhoair(i) = -delp(i,j,k)/(grav*delz(i,j,k)) ! moist air density
         enddo
      endif
      do i=is, ie
! The following form vectorizes better & more consistent with GFDL_MP
! SJL notes: Marshall-Palmer, dBZ = 200*precip**1.6, precip = 3.6e6*t1/rho_r*vtr  ! [mm/hr]
! GFDL_MP terminal fall speeds are used
! Date modified 20170701
! Account for excessively high cloud water -> autoconvert (diag only) excess cloud water
         t1 = rhoair(i)*max(qmin, q(i,j,k,rainwat)+dim(q(i,j,k,liq_wat), 1.0e-3))
         t2 = rhoair(i)*max(qmin, q(i,j,k,snowwat))
         if (graupel > 0) then
           t3 = rhoair(i)*max(qmin, q(i,j,k,graupel))
         else
           t3 = rhoair(i)*qmin
         endif
         denfac = sqrt(min(10., 1.2/rhoair(i)))
         vtr = max(1.e-3, vconr*denfac*exp(0.2   *log(t1/normr)))
         vtg = max(1.e-3, vcong*denfac*exp(0.125 *log(t3/normg)))
!     vts = max(1.e-3, vcons*denfac*exp(0.0625*log(t2/norms)))
         z_e = 200.*(exp(1.6*log(3.6e6*t1/rho_r*vtr)) + exp(1.6*log(3.6e6*t3/rho_g0*vtg))) + (factor_s/alpha)*t2*exp(0.75*log(t2/rn0_s))
!     z_e = 200.*(exp(1.6*log(3.6e6*t1/rho_r*vtr)) + exp(1.6*log(3.6e6*t3/rho_g*vtg)) + exp(1.6*log(3.6e6*t2/rho_s*vts)))
         dbz(i,j,k) = 10.*log10( max(0.01, z_e) )
      enddo
   enddo
   enddo

!$OMP parallel do default(shared)
   do j=js, je
      do k=mp_top+1, npz
         do i=is, ie
            maxdbz(i,j) = max(dbz(i,j,k), maxdbz(i,j))
         enddo
      enddo
   enddo

   do j=js, je
      do i=is, ie
         allmax = max(maxdbz(i,j), allmax)
      enddo
   enddo

 end subroutine dbzcalc
!
 subroutine max_vorticity_hy1(is, ie, js, je, km, vort, maxvorthy1)
   integer, intent(in):: is, ie, js, je, km
   real, intent(in), dimension(is:ie,js:je,km):: vort
   real, intent(inout), dimension(is:ie,js:je):: maxvorthy1
   integer i, j, k

   do j=js,je
      do i=is,ie
         maxvorthy1(i,j)=max(maxvorthy1(i,j),vort(i,j,km))
      enddo  ! i-loop
   enddo   ! j-loop
 end subroutine max_vorticity_hy1
!
! subroutine max_vorticity(is, ie, js, je, ng, km, zvir, sphum, delz, q, hydrostatic, &
!                          pt, peln, phis, grav, vort, maxvorthy1, maxvort, z_bot, z_top)
 subroutine max_vorticity(is, ie, js, je, ng, km, zvir, sphum, delz, q, hydrostatic, &
                          pt, peln, phis, grav, vort, maxvort, z_bot, z_top)
   integer, intent(in):: is, ie, js, je, ng, km, sphum
   real, intent(in):: grav, zvir, z_bot, z_top
   real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt
   real, intent(in), dimension(is:ie,js:je,km):: vort
   real, intent(in):: delz(is-ng:ie+ng,js-ng:je+ng,km)
   real, intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
   real, intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
   real, intent(in):: peln(is:ie,km+1,js:je)
   logical, intent(in):: hydrostatic
!   real, intent(inout), dimension(is:ie,js:je):: maxvorthy1,maxvort
   real, intent(inout), dimension(is:ie,js:je):: maxvort

   real:: rdg
   real, dimension(is:ie):: zh, dz, zh0
   integer i, j, k,klevel
   logical below(is:ie)

   rdg = rdgas / grav

   do j=js,je

      do i=is,ie
         zh(i) = 0.
         below(i) = .true.
         zh0(i) = 0.

     k_loop:do k=km,1,-1
            if ( hydrostatic ) then
                 dz(i) = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
            else
                 dz(i) = - delz(i,j,k)
            endif
            zh(i) = zh(i) + dz(i)
            if (zh(i) <= z_bot ) continue
            if (zh(i) > z_bot .and. below(i)) then
               maxvort(i,j) = max(maxvort(i,j),vort(i,j,k))
               below(i) = .false.
            elseif ( zh(i) < z_top ) then
               maxvort(i,j) = max(maxvort(i,j),vort(i,j,k))
            else
               maxvort(i,j) = max(maxvort(i,j),vort(i,j,k))
               EXIT k_loop 
            endif
         enddo k_loop
!      maxvorthy1(i,j)=max(maxvorthy1(i,j),vort(i,j,km))
      enddo  ! i-loop
   enddo   ! j-loop


 end subroutine max_vorticity

 subroutine max_uh(is, ie, js, je, ng, km, zvir, sphum, uphmax,uphmin,   &
                  w, vort, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
! !INPUT PARAMETERS:
   integer, intent(in):: is, ie, js, je, ng, km, sphum
   real, intent(in):: grav, zvir, z_bot, z_top
   real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, w
   real, intent(in), dimension(is:ie,js:je,km):: vort
   real, intent(in):: delz(is-ng:ie+ng,js-ng:je+ng,km)
   real, intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
   real, intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
   real, intent(in):: peln(is:ie,km+1,js:je)
   logical, intent(in):: hydrostatic
   real :: uh(is:ie,js:je)   ! unit: (m/s)**2
   real, intent(inout), dimension(is:ie,js:je):: uphmax,uphmin
! Coded by S.-J. Lin for CONUS regional climate simulations
! Modified for UH by LMH
!
   real:: rdg
   real, dimension(is:ie):: zh, dz, zh0
   integer i, j, k
   logical below(is:ie)

   rdg = rdgas / grav
   do j=js,je

      do i=is,ie
         zh(i) = 0.
         uh(i,j) = 0.
         below(i) = .true.
         zh0(i) = 0.

!        if ( phis(i,j)/grav < 1.E3 ) then
  k_loop:do k=km,1,-1
            if ( hydrostatic ) then
                 dz(i) = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
            else
                 dz(i) = - delz(i,j,k)
            endif
            zh(i) = zh(i) + dz(i)
            if (zh(i) <= z_bot ) continue
            if (zh(i) > z_bot .and. below(i)) then
               if(w(i,j,k).lt.0)then
                  uh(i,j) = 0.
                  EXIT k_loop
               endif 
               uh(i,j) = vort(i,j,k)*w(i,j,k)*(zh(i) - z_bot)
               below(i) = .false.
! Compute mean winds below z_top
            elseif ( zh(i) < z_top ) then
               if(w(i,j,k).lt.0)then
                  uh(i,j) = 0.
                  EXIT k_loop
               endif
               uh(i,j) = uh(i,j) + vort(i,j,k)*w(i,j,k)*dz(i)
            else
               if(w(i,j,k).lt.0)then
                  uh(i,j) = 0.
                  EXIT k_loop
               endif
               uh(i,j) = uh(i,j) + vort(i,j,k)*w(i,j,k)*(z_top - (zh(i)-dz(i)) )
               EXIT k_loop
            endif
         enddo k_loop
        if (uh(i,j) > uphmax(i,j)) then 
           uphmax(i,j) = uh(i,j)
        elseif (uh(i,j) < uphmin(i,j)) then 
           uphmin(i,j) = uh(i,j)
        endif 
      enddo  ! i-loop
   enddo   ! j-loop

 end subroutine max_uh 
 subroutine max_vv(is,ie,js,je,npz,ng,up2,dn2,pe,w)
! !INPUT PARAMETERS:
   integer, intent(in):: is, ie, js, je, ng, npz
   integer :: i,j,k
   real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,npz):: w
   real, intent(in):: pe(is-1:ie+1,npz+1,js-1:je+1)
   real, intent(inout), dimension(is:ie,js:je):: up2,dn2
      do j=js,je
         do i=is,ie
            do k=3,npz
              if (pe(i,k,j) >= 100.e2) then 
                  up2(i,j) = max(up2(i,j),w(i,j,k))
                  dn2(i,j) = min(dn2(i,j),w(i,j,k))
              endif
             enddo
         enddo
      enddo
 end subroutine max_vv 

 subroutine fv_diag_init_gn(Atm)
   type(fv_atmos_type), intent(inout), target :: atm
   
   if (atm%grid_Number > 1) then
      write(gn,"(A2,I1)") " g", atm%grid_number
   else
      gn = ""
   end if
   
 end subroutine fv_diag_init_gn

    subroutine getcape( nk , p , t , dz, q, the, cape , cin, source_in )
    implicit none

    integer, intent(in) :: nk
    real, dimension(nk), intent(in) :: p,t,dz,q,the
    real, intent(out) :: cape,cin
    integer, intent(IN), OPTIONAL :: source_in

!-----------------------------------------------------------------------
!
!  getcape - a fortran90 subroutine to calculate Convective Available
!            Potential Energy (CAPE) from a sounding.
!
!  Version 1.02                           Last modified:  10 October 2008
!
!  Author:  George H. Bryan
!           Mesoscale and Microscale Meteorology Division
!           National Center for Atmospheric Research
!           Boulder, Colorado, USA
!           gbryan@ucar.edu
!
!  Disclaimer:  This code is made available WITHOUT WARRANTY.
!
!  References:  Bolton (1980, MWR, p. 1046) (constants and definitions)
!               Bryan and Fritsch (2004, MWR, p. 2421) (ice processes)
!
!-----------------------------------------------------------------------
!
!  Input:     nk - number of levels in the sounding (integer)
!
!           p - one-dimensional array of pressure (Pa) (real)
!
!           t - one-dimensional array of temperature (K) (real)
!
!           dz - one-dimensional array of height thicknesses (m) (real)
!
!           q - one-dimensional array of specific humidity (kg/kg) (real)
!
!          source - source parcel:
!                   1 = surface (default)
!                   2 = most unstable (max theta-e)
!                   3 = mixed-layer (specify ml_depth)
!
!  Output:  cape - Convective Available Potential Energy (J/kg) (real)
!
!            cin - Convective Inhibition (J/kg) (real)
!
!-----------------------------------------------------------------------
!  User options:

    real, parameter :: pinc = 10000.0
                                        ! (smaller number yields more accurate
                                        !  results,larger number makes code 
                                        !  go faster)


    real, parameter :: ml_depth =  200.0

    integer, parameter :: adiabat = 1

!-----------------------------------------------------------------------
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----------------------------------------------------------------------
!            No need to modify anything below here:
!-----------------------------------------------------------------------

    integer :: source = 1
    logical :: doit,ice,cloud,not_converged
    integer :: k,kmin,n,nloop,i,orec
    real, dimension(nk) :: pi,th,thv,z,pt,pb,pc,pn,ptv

    real :: maxthe,parea,narea,lfc
    real :: th1,p1,t1,qv1,ql1,qi1,b1,pi1,thv1,qt,dp,frac
    real :: th2,p2,t2,qv2,ql2,qi2,b2,pi2,thv2
    real :: thlast,fliq,fice,tbar,qvbar,qlbar,qibar,lhv,lhs,lhf,rm,cpm
    real*8 :: avgth,avgqv

!-----------------------------------------------------------------------

    real, parameter :: g     = 9.81
    real, parameter :: p00   = 100000.0
    real, parameter :: cp    = 1005.7
    real, parameter :: rd    = 287.04
    real, parameter :: rv    = 461.5
    real, parameter :: xlv   = 2501000.0
    real, parameter :: xls   = 2836017.0
    real, parameter :: t0    = 273.15
    real, parameter :: cpv   = 1875.0
    real, parameter :: cpl   = 4190.0
    real, parameter :: cpi   = 2118.636
    real, parameter :: lv1   = xlv+(cpl-cpv)*t0
    real, parameter :: lv2   = cpl-cpv
    real, parameter :: ls1   = xls+(cpi-cpv)*t0
    real, parameter :: ls2   = cpi-cpv

    real, parameter :: rp00  = 1.0/p00
    real, parameter :: eps   = rd/rv
    real, parameter :: reps  = rv/rd
    real, parameter :: rddcp = rd/cp
    real, parameter :: cpdrd = cp/rd
    real, parameter :: cpdg  = cp/g

    real, parameter :: converge = 0.1

    integer, parameter :: debug_level =   0

    if (present(source_in)) source = source_in

!-----------------------------------------------------------------------

!---- convert p,t to mks units; get pi,th,thv ----!

    do k=1,nk
       pi(k) = (p(k)*rp00)**rddcp
       th(k) = t(k)/pi(k)
      thv(k) = th(k)*(1.0+reps*q(k))/(1.0+q(k))
    enddo

!---- get height using the hydrostatic equation ----!

    z(nk) = 0.5*dz(nk)
    do k=nk-1,1,-1
      z(k) = z(k+1) + 0.5*(dz(k+1)+dz(k))
    enddo

!---- find source parcel ----!

  IF(source.eq.1)THEN
    ! use surface parcel
    kmin = nk

  ELSEIF(source.eq.2)THEN
    ! use most unstable parcel (max theta-e)

    IF(p(1).lt.50000.0)THEN
      ! first report is above 500 mb ... just use the first level reported
      kmin = nk
      maxthe = the(nk)
    ELSE
      ! find max thetae below 500 mb
      maxthe = 0.0
      do k=nk,1,-1
        if(p(k).ge.50000.0)then
          if( the(nk).gt.maxthe )then
            maxthe = the(nk)
            kmin = k
          endif
        endif
      enddo
    ENDIF
    if(debug_level.ge.100) print *,'  kmin,maxthe = ',kmin,maxthe

!!$  ELSEIF(source.eq.3)THEN
!!$    ! use mixed layer
!!$
!!$    IF( dz(nk).gt.ml_depth )THEN
!!$      ! the second level is above the mixed-layer depth:  just use the
!!$      ! lowest level
!!$
!!$      avgth = th(nk)
!!$      avgqv = q(nk)
!!$      kmin = nk
!!$
!!$    ELSEIF( z(1).lt.ml_depth )THEN
!!$      ! the top-most level is within the mixed layer:  just use the
!!$      ! upper-most level (not 
!!$
!!$      avgth = th(1)
!!$      avgqv = q(1)
!!$      kmin = 1
!!$
!!$    ELSE
!!$      ! calculate the mixed-layer properties:
!!$
!!$      avgth = 0.0
!!$      avgqv = 0.0
!!$      k = nk-1
!!$      if(debug_level.ge.100) print *,'  ml_depth = ',ml_depth
!!$      if(debug_level.ge.100) print *,'  k,z,th,q:'
!!$      if(debug_level.ge.100) print *,nk,z(nk),th(nk),q(nk)
!!$
!!$      do while( (z(k).le.ml_depth) .and. (k.ge.1) )
!!$
!!$        if(debug_level.ge.100) print *,k,z(k),th(k),q(k)
!!$
!!$        avgth = avgth + dz(k)*th(k)
!!$        avgqv = avgqv + dz(k)*q(k)
!!$
!!$        k = k - 1
!!$
!!$      enddo
!!$
!!$      th2 = th(k+1)+(th(k)-th(k+1))*(ml_depth-z(k-1))/dz(k)
!!$      qv2 =  q(k+1)+( q(k)- q(k+1))*(ml_depth-z(k-1))/dz(k)
!!$
!!$      if(debug_level.ge.100) print *,999,ml_depth,th2,qv2
!!$
!!$      avgth = avgth + 0.5*(ml_depth-z(k-1))*(th2+th(k-1))
!!$      avgqv = avgqv + 0.5*(ml_depth-z(k-1))*(qv2+q(k-1))
!!$
!!$      if(debug_level.ge.100) print *,k,z(k),th(k),q(k)
!!$
!!$      avgth = avgth/ml_depth
!!$      avgqv = avgqv/ml_depth
!!$
!!$      kmin = nk
!!$
!!$    ENDIF
!!$
!!$    if(debug_level.ge.100) print *,avgth,avgqv

  ELSE

    print *
    print *,'  Unknown value for source'
    print *
    print *,'  source = ',source
    print *
    call mpp_error(fatal, " Unknown CAPE source")

  ENDIF

!---- define parcel properties at initial location ----!
    narea = 0.0

  if( (source.eq.1).or.(source.eq.2) )then
    k    = kmin
    th2  = th(kmin)
    pi2  = pi(kmin)
    p2   = p(kmin)
    t2   = t(kmin)
    thv2 = thv(kmin)
    qv2  = q(kmin)
    b2   = 0.0
  elseif( source.eq.3 )then
    k    = kmin
    th2  = avgth
    qv2  = avgqv
    thv2 = th2*(1.0+reps*qv2)/(1.0+qv2)
    pi2  = pi(kmin)
    p2   = p(kmin)
    t2   = th2*pi2
    b2   = g*( thv2-thv(kmin) )/thv(kmin)
  endif

    ql2 = 0.0
    qi2 = 0.0
    qt  = qv2

    cape = 0.0
    cin  = 0.0
    lfc  = 0.0

    doit = .true.
    cloud = .false.
    if(adiabat.eq.1.or.adiabat.eq.2)then
      ice = .false.
    else
      ice = .true.
    endif

!      the = getthe(p2,t2,t2,qv2)
!      if(debug_level.ge.100) print *,'  the = ',the

!---- begin ascent of parcel ----!

      if(debug_level.ge.100)then
        print *,'  Start loop:'
        print *,'  p2,th2,qv2 = ',p2,th2,qv2
      endif

    do while( doit .and. (k.gt.1) )

        k = k-1
       b1 =  b2

       dp = p(k)-p(k-1)

      if( dp.lt.pinc )then
        nloop = 1
      else
        nloop = 1 + int( dp/pinc )
        dp = dp/float(nloop)
      endif

      do n=1,nloop

         p1 =  p2
         t1 =  t2
        pi1 = pi2
        th1 = th2
        qv1 = qv2
        ql1 = ql2
        qi1 = qi2
        thv1 = thv2

        p2 = p2 - dp
        pi2 = (p2*rp00)**rddcp

        thlast = th1
        i = 0
        not_converged = .true.

        do while( not_converged )
          i = i + 1
          t2 = thlast*pi2
          if(ice)then
            fliq = max(min((t2-233.15)/(273.15-233.15),1.0),0.0)
            fice = 1.0-fliq
          else
            fliq = 1.0
            fice = 0.0
          endif
          qv2 = min( qt , fliq*getqvs(p2,t2) + fice*getqvi(p2,t2) )
          qi2 = max( fice*(qt-qv2) , 0.0 )
          ql2 = max( qt-qv2-qi2 , 0.0 )

          tbar  = 0.5*(t1+t2)
          qvbar = 0.5*(qv1+qv2)
          qlbar = 0.5*(ql1+ql2)
          qibar = 0.5*(qi1+qi2)

          lhv = lv1-lv2*tbar
          lhs = ls1-ls2*tbar
          lhf = lhs-lhv

          rm=rd+rv*qvbar
          cpm=cp+cpv*qvbar+cpl*qlbar+cpi*qibar
          th2=th1*exp(  lhv*(ql2-ql1)/(cpm*tbar)     &
                       +lhs*(qi2-qi1)/(cpm*tbar)     &
                       +(rm/cpm-rd/cp)*alog(p2/p1) )

          if(i.gt.90) print *,i,th2,thlast,th2-thlast
          if(i.gt.100)then
            print *
            print *,'  Error:  lack of convergence'
            print *
            print *,'  ... stopping iteration '
            print *
            stop 1001
          endif
          if( abs(th2-thlast).gt.converge )then
            thlast=thlast+0.3*(th2-thlast)
          else
            not_converged = .false.
          endif
        enddo

        ! Latest pressure increment is complete.  Calculate some
        ! important stuff:

        if( ql2.ge.1.0e-10 ) cloud = .true.

        IF(adiabat.eq.1.or.adiabat.eq.3)THEN
          ! pseudoadiabat
          qt  = qv2
          ql2 = 0.0
          qi2 = 0.0
        ELSEIF(adiabat.le.0.or.adiabat.ge.5)THEN
          print *
          print *,'  Undefined adiabat'
          print *
          stop 10000
        ENDIF

      enddo

      thv2 = th2*(1.0+reps*qv2)/(1.0+qv2+ql2+qi2)
        b2 = g*( thv2-thv(k) )/thv(k)

!      the = getthe(p2,t2,t2,qv2)

      ! Get contributions to CAPE and CIN:

      if( (b2.ge.0.0) .and. (b1.lt.0.0) )then
        ! first trip into positive area
        !ps = p(k-1)+(p(k)-p(k-1))*(0.0-b1)/(b2-b1)
        frac = b2/(b2-b1)
        parea =  0.5*b2*dz(k)*frac
        narea = narea-0.5*b1*dz(k)*(1.0-frac)
        if(debug_level.ge.200)then
          print *,'      b1,b2 = ',b1,b2
          !print *,'      p1,ps,p2 = ',p(k-1),ps,p(k)
          print *,'      frac = ',frac
          print *,'      parea = ',parea
          print *,'      narea = ',narea
        endif
        cin  = cin  + narea
        narea = 0.0
      elseif( (b2.lt.0.0) .and. (b1.gt.0.0) )then
        ! first trip into neg area
        !ps = p(k-1)+(p(k)-p(k-1))*(0.0-b1)/(b2-b1)
        frac = b1/(b1-b2)
        parea =  0.5*b1*dz(k)*frac
        narea = -0.5*b2*dz(k)*(1.0-frac)
        if(debug_level.ge.200)then
          print *,'      b1,b2 = ',b1,b2
          !print *,'      p1,ps,p2 = ',p(k-1),ps,p(k)
          print *,'      frac = ',frac
          print *,'      parea = ',parea
          print *,'      narea = ',narea
        endif
      elseif( b2.lt.0.0 )then
        ! still collecting negative buoyancy
        parea =  0.0
        narea = narea-0.5*dz(k)*(b1+b2)
      else
        ! still collecting positive buoyancy
        parea =  0.5*dz(k)*(b1+b2)
        narea =  0.0
      endif

      cape = cape + max(0.0,parea)

      if(debug_level.ge.200)then
        write(6,102) p2,b1,b2,cape,cin,cloud
102     format(5(f13.4),2x,l1)
      endif

      if( (p(k).le.10000.0).and.(b2.lt.0.0) )then
        ! stop if b < 0 and p < 100 mb
        doit = .false.
      endif

    enddo

!---- All done ----!

    return
  end subroutine getcape
!-----------------------------------------------------------------------


  real function getqvs(p,t)
    implicit none

    real :: p,t,es

    real, parameter :: eps = 287.04/461.5

    es = 611.2*exp(17.67*(t-273.15)/(t-29.65))
    getqvs = eps*es/(p-es)

    return
  end function getqvs
!-----------------------------------------------------------------------


  real function getqvi(p,t)
    implicit none

    real :: p,t,es

    real, parameter :: eps = 287.04/461.5

    es = 611.2*exp(21.8745584*(t-273.15)/(t-7.66))
    getqvi = eps*es/(p-es)

    return
  end function getqvi
!-----------------------------------------------------------------------

end module fv_diagnostics_mod