module spectral_dynamics_mod

#ifdef INTERNAL_FILE_NML
use mpp_mod, only: input_nml_file
#else
use fms_mod, only: open_namelist_file
#endif

use                fms_mod, only: mpp_pe, mpp_root_pe, error_mesg, NOTE, FATAL, write_version_number, stdlog, &
                                  close_file, open_restart_file, file_exist, set_domain,                      &
                                  read_data, write_data, check_nml_error, lowercase, uppercase, mpp_npes,     &
                                  field_size

use          constants_mod, only: rdgas, rvgas, grav, cp_air, omega, radius, pi

use       time_manager_mod, only: time_type, get_time, set_time, get_calendar_type, NO_CALENDAR, &
                                  get_date, interval_alarm, operator( - ), operator( + )

use      field_manager_mod, only: MODEL_ATMOS, parse

use     tracer_manager_mod, only: get_number_tracers, query_method, get_tracer_index, get_tracer_names, NO_TRACER

use       diag_manager_mod, only: diag_axis_init, register_diag_field, register_static_field, send_data

use         transforms_mod, only: transforms_init,         transforms_end,            &
                                  get_grid_boundaries,     area_weighted_global_mean, &
                                  compute_gradient_cos,    trans_spherical_to_grid,   &
                                  trans_grid_to_spherical, divide_by_cos,             &
                                  triangular_truncation,   rhomboidal_truncation,     &
                                  compute_laplacian,       get_eigen_laplacian,       &
                                  get_spherical_wave,      get_sin_lat,               &
                                  vor_div_from_uv_grid,    uv_grid_from_vor_div,      &
                                  horizontal_advection,    get_grid_domain,           &
                                  get_spec_domain,         grid_domain,               &
                                  spectral_domain,         get_deg_lon, get_deg_lat

use     vert_advection_mod, only: vert_advection, SECOND_CENTERED, FOURTH_CENTERED, VAN_LEER_LINEAR, FINITE_VOLUME_PARABOLIC, &
                                  ADVECTIVE_FORM

use           implicit_mod, only: implicit_init, implicit_end, implicit_correction

use   press_and_geopot_mod, only: press_and_geopot_init, press_and_geopot_end, &
                                  pressure_variables, compute_geopotential, compute_pressures_and_heights

use   spectral_damping_mod, only: spectral_damping_init, spectral_damping_end, compute_spectral_damping, &
                                  compute_spectral_damping_vor, compute_spectral_damping_div

use           leapfrog_mod, only: leapfrog, leapfrog_2level_A, leapfrog_2level_B

use       fv_advection_mod, only: fv_advection_init, fv_advection_end, a_grid_horiz_advection
                                  
use    water_borrowing_mod, only: water_borrowing

use    global_integral_mod, only: mass_weighted_global_integral

use spectral_init_cond_mod, only: spectral_init_cond

use        tracer_type_mod, only: tracer_type, tracer_type_version, tracer_type_tagname

use every_step_diagnostics_mod, only: every_step_diagnostics_init, every_step_diagnostics, every_step_diagnostics_end

use        mpp_domains_mod, only: mpp_global_field

use       polvani_2004_mod, only: polvani_2004
use       polvani_2007_mod, only: polvani_2007, polvani_2007_tracer_init, get_polvani_2007_tracers
use   jablonowski_2006_mod, only: jablonowski_2006

!===============================================================================================
implicit none
private
!===============================================================================================

public :: spectral_dynamics_init, spectral_dynamics, spectral_dynamics_end, get_num_levels
public :: get_use_virtual_temperature, get_reference_sea_level_press, get_surf_geopotential
public :: get_pk_bk, complete_robert_filter, complete_update_of_future
public :: get_axis_id, spectral_diagnostics, get_initial_fields

!===============================================================================================

character(len=128), parameter :: version = '$Id: spectral_dynamics.F90,v 19.0 2012/01/06 20:35:46 fms Exp $'
character(len=128), parameter :: tagname = '$Name: tikal $'

!===============================================================================================
! variables needed for diagnostics
integer :: id_ps, id_u, id_v, id_t, id_vor, id_div, id_omega, id_wspd, id_slp
integer :: id_pres_full, id_pres_half, id_zfull, id_zhalf, id_vort_norm, id_EKE
integer :: id_uu, id_vv, id_tt, id_omega_omega, id_uv, id_omega_t
integer, allocatable, dimension(:) :: id_tr
real :: gamma, expf, expf_inverse
character(len=8) :: mod_name = 'dynamics'
integer, dimension(4) :: axis_id
!===============================================================================================

integer, parameter :: num_time_levels = 2

logical :: module_is_initialized = .false.
logical :: dry_model
logical :: robert_complete_for_tracers=.true., robert_complete_for_fields=.true. ! Needed only for error checks during code development

type(time_type) :: Time_step, Alarm_time, Alarm_interval ! Used to determine when it is time to print global integrals.

real,    allocatable, dimension(:) :: sin_lat, coriolis
real,    allocatable, dimension(:) :: pk, bk, dpk, dbk

complex, allocatable, dimension(:,:,:,:)   :: vors, divs, ts ! last dimension is for time level
complex, allocatable, dimension(:,:,:  )   :: ln_ps          ! last dimension is for time level
complex, allocatable, dimension(:,:,:,:,:) :: spec_tracers   ! 4'th dimension is for time level, last dimension is for tracer number

real, allocatable, dimension(:,:,:    ) :: psg               ! last dimension is for time level
real, allocatable, dimension(:,:,:,:  ) :: ug, vg, tg        ! last dimension is for time level
real, allocatable, dimension(:,:,:,:,:) :: grid_tracers      ! 4'th dimension is for time level, last dimension is for tracer number
real, allocatable, dimension(:,:      ) :: surf_geopotential
real, allocatable, dimension(:,:,:    ) :: vorg, divg        ! no time levels needed

integer, allocatable, dimension(:) :: tracer_vert_advect_scheme

real    :: virtual_factor, dt_real
integer :: pe, npes, num_tracers, nhum, t_vert_advect_scheme, uv_vert_advect_scheme, step_number
real    :: mean_energy_previous, mean_water_previous, mean_surf_press_previous
integer :: ms, me, ns, ne, is, ie, js, je
integer :: previous, current, future

character(len=32), parameter :: default_representation = 'spectral'
character(len=32), parameter :: default_advect_vert    = 'second_centered'
character(len=32), parameter :: default_hole_filling   = 'off'

!===============================================================================================
! namelist variables

logical :: do_mass_correction     = .true. , &
           do_water_correction    = .true. , &
           do_energy_correction   = .true. , &
           use_virtual_temperature= .false., &
           use_implicit           = .true.,  &
           triang_trunc           = .true.

integer :: damping_order       = 2, &
           damping_order_vor   =-1, &
           damping_order_div   =-1, &
           lon_max             = 128, & ! T42
           lat_max             = 64,  & ! T42
           num_fourier         = 42,  & ! T42
           num_spherical       = 43,  & ! T42
           fourier_inc         = 1,   &
           num_levels          = 18,  &
           num_steps           = 1

integer, dimension(2) ::  print_interval=(/1,0/)

character(len=64) :: vert_coord_option      = 'even_sigma',   &
                     damping_option         = 'resolution_dependent', &
                     vert_advect_uv         = default_advect_vert,    &
                     vert_advect_t          = default_advect_vert,    &
                     vert_difference_option = 'simmons_and_burridge', &
                     initial_state_option   = 'quiescent'

real    :: damping_coeff       = 1.15740741e-4, & ! (one tenth day)**-1
           damping_coeff_vor   = -1., &
           damping_coeff_div   = -1., &
           eddy_sponge_coeff   =  0., &
           zmu_sponge_coeff    =  0., &
           zmv_sponge_coeff    =  0., &
           robert_coeff        = .04, &
           alpha_implicit      = .5,  &
           longitude_origin    =  0., &
           scale_heights       =  4., &
           surf_res            = .1,  &
           p_press             = .1,  &
           p_sigma             = .3,  &
           exponent            = 2.5, &
         ocean_topog_smoothing = .93, &
           initial_sphum       = 0.0, &
     reference_sea_level_press =  101325.
!===============================================================================================

real, dimension(2) :: valid_range_t = (/100.,500./)

namelist /spectral_dynamics_nml/ use_virtual_temperature, damping_option,                            &
                                 damping_order, damping_coeff, damping_order_vor, damping_coeff_vor, &
                                 damping_order_div, damping_coeff_div, do_mass_correction,           &
                                 do_water_correction, do_energy_correction, vert_advect_uv,          &
                                 vert_advect_t, use_implicit, longitude_origin, robert_coeff,        &
                                 alpha_implicit, vert_difference_option,                             &
                                 reference_sea_level_press, lon_max, lat_max, num_levels,            &
                                 num_fourier, num_spherical, fourier_inc, triang_trunc,              &
                                 vert_coord_option, scale_heights, surf_res,      &
                                 p_press, p_sigma, exponent, ocean_topog_smoothing, initial_sphum,   &
                                 valid_range_t, eddy_sponge_coeff, zmu_sponge_coeff, zmv_sponge_coeff, &
                                 print_interval, num_steps, initial_state_option
                                 
contains

!===============================================================================================

subroutine spectral_dynamics_init(Time, Time_step_in, tracer_attributes, dry_model_out, nhum_out, ocean_mask)

type(time_type), intent(in) :: Time, Time_step_in
type(tracer_type), intent(inout), dimension(:) :: tracer_attributes
logical, intent(out) :: dry_model_out
integer, intent(out) :: nhum_out
logical, optional, intent(in), dimension(:,:) :: ocean_mask

integer :: num_total_wavenumbers, unit, k, seconds, days, ierr, io, ntr, nsphum, nmix_rat
logical :: south_to_north = .true.
real    :: ref_surf_p_implicit, robert_coeff_tracers

real,    allocatable, dimension(:,:) :: eigen
integer, allocatable, dimension(:,:) :: wavenumber
real,    allocatable, dimension(:)   :: glon_bnd, glat_bnd, ref_temperature_implicit
character(len=32) :: scheme, params
character(len=128) :: tname, longname, units

! < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < >

if(module_is_initialized) return

#ifdef INTERNAL_FILE_NML
    read (input_nml_file, nml=spectral_dynamics_nml, iostat=io)
    ierr = check_nml_error(io, 'spectral_dynamics_nml')
#else
    unit = open_namelist_file()
    ierr=1
    do while (ierr /= 0)
      read(unit, nml=spectral_dynamics_nml, iostat=io, end=20)
      ierr = check_nml_error (io, 'spectral_dynamics_nml')
    enddo
20  call close_file (unit)
#endif

call write_version_number(version, tagname)
if(mpp_pe() == mpp_root_pe()) write (stdlog(), nml=spectral_dynamics_nml)
call write_version_number(version, tagname)

Time_step  = Time_step_in
Alarm_interval = set_time(print_interval(2), print_interval(1))
Alarm_time = Time + Alarm_interval

if(damping_order_vor == -1 ) damping_order_vor = damping_order
if(damping_order_div == -1 ) damping_order_div = damping_order
if(damping_coeff_vor == -1.) damping_coeff_vor = damping_coeff
if(damping_coeff_div == -1.) damping_coeff_div = damping_coeff

call check_dynamics_nml

if(uppercase(trim(vert_advect_uv)) == 'SECOND_CENTERED') then
  uv_vert_advect_scheme = SECOND_CENTERED
else if(uppercase(trim(vert_advect_uv)) == 'FOURTH_CENTERED') then
  uv_vert_advect_scheme = FOURTH_CENTERED
else if(uppercase(trim(vert_advect_uv)) == 'VAN_LEER_LINEAR') then
  uv_vert_advect_scheme = VAN_LEER_LINEAR
else if(uppercase(trim(vert_advect_uv)) == 'FINITE_VOLUME_PARABOLIC') then
  uv_vert_advect_scheme = FINITE_VOLUME_PARABOLIC
else
  call error_mesg('spectral_dynamics_init','"'//trim(vert_advect_uv)//'"'//' is not a valid value for vert_advect_uv.', FATAL)
endif

if(uppercase(trim(vert_advect_t)) == 'SECOND_CENTERED') then
  t_vert_advect_scheme = SECOND_CENTERED
else if(uppercase(trim(vert_advect_t)) == 'FOURTH_CENTERED') then
  t_vert_advect_scheme = FOURTH_CENTERED
else if(uppercase(trim(vert_advect_t)) == 'VAN_LEER_LINEAR') then
  t_vert_advect_scheme = VAN_LEER_LINEAR
else if(uppercase(trim(vert_advect_t)) == 'FINITE_VOLUME_PARABOLIC') then
  t_vert_advect_scheme = FINITE_VOLUME_PARABOLIC
else
  call error_mesg('spectral_dynamics_init','"'//trim(vert_advect_t)//'"'//' is not a valid value for vert_advect_t.', FATAL)
endif

! transforms_init must be called before the remaining 
! restart data can be read or the fields can be allocated.
! This is because transforms_init calls spec_mpp_init,
! which is where domains are determined.

call transforms_init(radius, lat_max, lon_max, num_fourier, fourier_inc, num_spherical, south_to_north=south_to_north, &
                     triang_trunc=triang_trunc, longitude_origin=longitude_origin)

call get_grid_domain(is, ie, js, je)
call get_spec_domain(ms, me, ns, ne)
call get_number_tracers(MODEL_ATMOS, num_prog=num_tracers)
call allocate_fields

do ntr=1,num_tracers

  call get_tracer_names(MODEL_ATMOS, ntr, tname, longname, units)
  tracer_attributes(ntr)%name = lowercase(tname)

  if(query_method('numerical_representation', MODEL_ATMOS, ntr, scheme)) then
    tracer_attributes(ntr)%numerical_representation = scheme
  else
    tracer_attributes(ntr)%numerical_representation = default_representation
  endif

  if(query_method('advect_vert', MODEL_ATMOS, ntr, scheme)) then
    tracer_attributes(ntr)%advect_vert = scheme
  else
    tracer_attributes(ntr)%advect_vert = default_advect_vert
  endif

  if(query_method('hole_filling', MODEL_ATMOS, ntr, scheme)) then
    tracer_attributes(ntr)%hole_filling = scheme
  else
    tracer_attributes(ntr)%hole_filling = default_hole_filling
  endif

  if(query_method('robert_filter', MODEL_ATMOS, ntr, scheme, params)) then
    if(uppercase(scheme) == 'OFF') then
      tracer_attributes(ntr)%robert_coeff = 0.0
    else
      if(parse(params,'robert_coeff', robert_coeff_tracers) == 1) then
        tracer_attributes(ntr)%robert_coeff = robert_coeff_tracers
      else
        tracer_attributes(ntr)%robert_coeff = robert_coeff
      endif
    endif
  else
    tracer_attributes(ntr)%robert_coeff = robert_coeff
  endif

  if(trim(tracer_attributes(ntr)%numerical_representation) == 'spectral') then
    tracer_attributes(ntr)%advect_horiz = 'spectral'
  else if(trim(tracer_attributes(ntr)%numerical_representation) == 'grid') then
    tracer_attributes(ntr)%advect_horiz = 'van_leer'
  else
    call error_mesg('spectral_dynamics_init',trim(tracer_attributes(ntr)%numerical_representation)// &
           ' is an invalid numerical_representation', FATAL)
  endif

  if(trim(tracer_attributes(ntr)%numerical_representation) == 'grid') then
    if(trim(tracer_attributes(ntr)%hole_filling) == 'on') then
      call error_mesg('spectral_dynamics_init','Warning: hole_filling scheme = on will be ignored for grid tracer '// &
      tracer_attributes(ntr)%name,NOTE)
    endif
  endif

enddo

nsphum   = get_tracer_index(MODEL_ATMOS, 'sphum')
nmix_rat = get_tracer_index(MODEL_ATMOS, 'mix_rat')

if(nsphum == NO_TRACER) then
  if(nmix_rat == NO_TRACER) then
    nhum = 0
    dry_model = .true.
  else
    nhum = nmix_rat 
    dry_model = .false.
  endif
else
  if(nmix_rat == NO_TRACER) then
    nhum = nsphum
    dry_model = .false.
  else
    call error_mesg('spectral_dynamics_init','sphum and mix_rat cannot both be specified as tracers at the same time', FATAL)
  endif
endif
dry_model_out = dry_model
nhum_out = nhum

allocate(tracer_vert_advect_scheme(num_tracers))

do ntr=1,num_tracers
  if(uppercase(trim(tracer_attributes(ntr)%advect_vert)) == 'SECOND_CENTERED') then
    tracer_vert_advect_scheme(ntr) = SECOND_CENTERED
  else if(uppercase(trim(tracer_attributes(ntr)%advect_vert)) == 'FOURTH_CENTERED') then
    tracer_vert_advect_scheme(ntr) = FOURTH_CENTERED
  else if(uppercase(trim(tracer_attributes(ntr)%advect_vert)) == 'VAN_LEER_LINEAR') then
    tracer_vert_advect_scheme(ntr) = VAN_LEER_LINEAR
  else if(uppercase(trim(tracer_attributes(ntr)%advect_vert)) == 'FINITE_VOLUME_PARABOLIC') then
    tracer_vert_advect_scheme(ntr) = FINITE_VOLUME_PARABOLIC
  else
    call error_mesg('spectral_dynamics_init',trim(tracer_attributes(ntr)%advect_vert)// &
         ' is not a valid vertical advection scheme for tracers. Check your field_table.', FATAL)
  endif
enddo

call read_restart_or_do_coldstart(tracer_attributes, ocean_mask)

call press_and_geopot_init(pk, bk, use_virtual_temperature, vert_difference_option)

call spectral_diagnostics_init(Time)

call every_step_diagnostics_init(Time, lon_max, lat_max, num_levels, reference_sea_level_press)

if(do_water_correction .and. .not.do_mass_correction) then
  call error_mesg('spectral_dynamics_init', 'water_correction requires mass_correction', FATAL)
endif

if(do_energy_correction .and. .not.do_mass_correction) then
  call error_mesg('spectral_dynamics_init', 'energy_correction requires mass_correction', FATAL)
endif

pe   = mpp_pe()
npes = mpp_npes()

if(triang_trunc) then
  num_total_wavenumbers = num_spherical - 1
else
  num_total_wavenumbers = num_spherical - 1 + fourier_inc*num_fourier
end if

if(use_virtual_temperature) then
  virtual_factor = (rvgas/rdgas) - 1.0
end if

allocate (sin_lat (js:je))
allocate (coriolis(js:je))

call get_sin_lat (sin_lat)

coriolis = 2*omega*sin_lat

allocate (glon_bnd (lon_max + 1))
allocate (glat_bnd (lat_max + 1))

call get_grid_boundaries(glon_bnd, glat_bnd, global=.true.)
call fv_advection_init  (lon_max, lat_max, glat_bnd, 360./fourier_inc)

deallocate (glat_bnd)
deallocate (glon_bnd)

allocate(dpk(size(pk,1)-1))
allocate(dbk(size(bk,1)-1))

do k=1,size(pk,1)-1
  dpk(k) = pk(k+1) - pk(k)
  dbk(k) = bk(k+1) - bk(k)
enddo

call spectral_damping_init(damping_coeff, damping_order, damping_option, num_fourier, num_spherical, &
                           num_levels, eddy_sponge_coeff, zmu_sponge_coeff, zmv_sponge_coeff,        &
                           damping_coeff_vor=damping_coeff_vor, damping_order_vor=damping_order_vor, &
                           damping_coeff_div=damping_coeff_div, damping_order_div=damping_order_div)

if(use_implicit) then
  allocate(wavenumber(0:num_fourier,0:num_spherical))
  allocate(     eigen(0:num_fourier,0:num_spherical))
  allocate (ref_temperature_implicit(num_levels))
  ref_temperature_implicit(:) = 300.
  ref_surf_p_implicit = reference_sea_level_press

  call get_spherical_wave(wavenumber)
  call get_eigen_laplacian(eigen)

  call implicit_init(pk, bk, ref_temperature_implicit, ref_surf_p_implicit, num_total_wavenumbers, &
                     eigen, wavenumber, alpha_implicit, vert_difference_option)

  deallocate(eigen, wavenumber, ref_temperature_implicit)
endif

call set_domain(grid_domain)

call get_time(Time_step, seconds, days)
dt_real = 86400*days + seconds

module_is_initialized = .true.
return
end subroutine spectral_dynamics_init

!===============================================================================================
subroutine read_restart_or_do_coldstart(tracer_attributes, ocean_mask)

! For backward compatibility, this routine has the capability
! to read native data restart files written by inchon code.

type(tracer_type), intent(inout), dimension(:) :: tracer_attributes
logical, optional, intent(in), dimension(:,:) :: ocean_mask

integer :: m, n, k, nt, ntr
integer, dimension(4) :: siz
real, dimension(ms:me, ns:ne, num_levels) :: real_part, imag_part
character(len=64) :: file, tr_name
character(len=4) :: ch1,ch2,ch3,ch4,ch5,ch6

file = 'INPUT/spectral_dynamics.res.nc'
if(file_exist(trim(file))) then
  call field_size(trim(file), 'vors_real', siz)
  if(num_fourier /= siz(1)-1 .or. num_spherical /= siz(2)-1 .or. num_levels /= siz(3)) then
    write(ch1,'(i4)') siz(1)-1
    write(ch2,'(i4)') siz(2)-1
    write(ch3,'(i4)') siz(3)
    write(ch4,'(i4)') num_fourier
    write(ch5,'(i4)') num_spherical
    write(ch6,'(i4)') num_levels
    call error_mesg('spectral_dynamics_init','Resolution of restart data does not match resolution specified on namelist.'// &
    ' Restart data: num_fourier='//ch1//', num_spherical='//ch2//', num_levels='//ch3// &
       '  Namelist: num_fourier='//ch4//', num_spherical='//ch5//', num_levels='//ch6, FATAL)
  endif
  call field_size(trim(file), 'ug', siz)
  if(lon_max /= siz(1) .or. lat_max /= siz(2)) then
    write(ch1,'(i4)') siz(1)
    write(ch2,'(i4)') siz(2)
    write(ch3,'(i4)') lon_max
    write(ch4,'(i4)') lat_max
    call error_mesg('spectral_dynamics_init','Resolution of restart data does not match resolution specified on namelist.'// &
    ' Restart data: lon_max='//ch1//', lat_max='//ch2//'  Namelist: lon_max='//ch3//', lat_max='//ch4, FATAL)
  endif
  call read_data(trim(file), 'previous', previous, no_domain=.true.)
  call read_data(trim(file), 'current',  current,  no_domain=.true.)
  call read_data(trim(file), 'pk', pk, no_domain=.true.)
  call read_data(trim(file), 'bk', bk, no_domain=.true.)
  do nt=1,num_time_levels
    call read_data(trim(file), 'vors_real',  real_part, spectral_domain, timelevel=nt)
    call read_data(trim(file), 'vors_imag',  imag_part, spectral_domain, timelevel=nt)
    do k=1,num_levels; do n=ns,ne; do m=ms,me
      vors(m,n,k,nt) = cmplx(real_part(m,n,k),imag_part(m,n,k))
    enddo; enddo; enddo
    call read_data(trim(file), 'divs_real',  real_part, spectral_domain, timelevel=nt)
    call read_data(trim(file), 'divs_imag',  imag_part, spectral_domain, timelevel=nt)
    do k=1,num_levels; do n=ns,ne; do m=ms,me
      divs(m,n,k,nt) = cmplx(real_part(m,n,k),imag_part(m,n,k))
    enddo; enddo; enddo
    call read_data(trim(file), 'ts_real',  real_part, spectral_domain, timelevel=nt)
    call read_data(trim(file), 'ts_imag',  imag_part, spectral_domain, timelevel=nt)
    do k=1,num_levels; do n=ns,ne; do m=ms,me
      ts(m,n,k,nt) = cmplx(real_part(m,n,k),imag_part(m,n,k))
    enddo; enddo; enddo
    call read_data(trim(file), 'ln_ps_real', real_part(:,:,1), spectral_domain, timelevel=nt)
    call read_data(trim(file), 'ln_ps_imag', imag_part(:,:,1), spectral_domain, timelevel=nt)
    do n=ns,ne; do m=ms,me
      ln_ps(m,n,nt) = cmplx(real_part(m,n,1),imag_part(m,n,1))
    enddo; enddo
    call read_data(trim(file), 'ug',   ug(:,:,:,nt), grid_domain, timelevel=nt)
    call read_data(trim(file), 'vg',   vg(:,:,:,nt), grid_domain, timelevel=nt)
    call read_data(trim(file), 'tg',   tg(:,:,:,nt), grid_domain, timelevel=nt)
    call read_data(trim(file), 'psg', psg(:,:,  nt), grid_domain, timelevel=nt)
    do ntr = 1,num_tracers
      tr_name = trim(tracer_attributes(ntr)%name)
      call read_data(trim(file), trim(tr_name), grid_tracers(:,:,:,nt,ntr), grid_domain, timelevel=nt)
      if(uppercase(trim(tracer_attributes(ntr)%numerical_representation)) == 'SPECTRAL') then
        call read_data(trim(file), trim(tr_name)//'_real', real_part, spectral_domain, timelevel=nt)
        call read_data(trim(file), trim(tr_name)//'_imag', imag_part, spectral_domain, timelevel=nt)
        do k=1,num_levels; do n=ns,ne; do m=ms,me
          spec_tracers(m,n,k,nt,ntr) = cmplx(real_part(m,n,k),imag_part(m,n,k))
        enddo; enddo; enddo
      endif
    enddo ! loop over tracers
  enddo ! loop over time levels
  call read_data(trim(file), 'vorg', vorg, grid_domain)
  call read_data(trim(file), 'divg', divg, grid_domain)
  call read_data(trim(file), 'surf_geopotential', surf_geopotential, grid_domain)
else if(file_exist('INPUT/spectral_dynamics.res')) then
  call error_mesg('spectral_dynamics_init', &
              'Binary restart file, INPUT/spectral_dynamics.res, is not supported by this version of spectral_dynamics.f90',FATAL)
else
  do ntr = 1,num_tracers
    if(trim(tracer_attributes(ntr)%name) == 'sphum') then
      grid_tracers(:,:,:,:,ntr) = initial_sphum
    else if(trim(tracer_attributes(ntr)%name) == 'mix_rat') then
      grid_tracers(:,:,:,:,ntr) = 0.
    else
      grid_tracers(:,:,:,:,ntr) = 0.
    endif   
  enddo

  previous = 1
  current  = 1
  if(initial_state_option == 'quiescent' .or. initial_state_option == 'input') then
    call spectral_init_cond(initial_state_option, tracer_attributes, reference_sea_level_press, triang_trunc, use_virtual_temperature,&
                            vert_coord_option, vert_difference_option, scale_heights, surf_res, p_press, p_sigma,  &
                            exponent, ocean_topog_smoothing, pk, bk,                                               &
                            vors(:,:,:,1), divs(:,:,:,1), ts(:,:,:,1), ln_ps(:,:,1), ug(:,:,:,1), vg(:,:,:,1),     &
                            tg(:,:,:,1), psg(:,:,1), vorg, divg, grid_tracers(:,:,:,1,:), surf_geopotential, ocean_mask)
  else if(initial_state_option == 'polvani_2004') then
    call polvani_2004(reference_sea_level_press, triang_trunc, vert_difference_option, pk, bk, &
                 vors(:,:,:,1), divs(:,:,:,1), ts(:,:,:,1), ln_ps(:,:,1), ug(:,:,:,1),  vg(:,:,:,1), &
                 tg(:,:,:,1), psg(:,:,1), vorg, divg, surf_geopotential)
  else if(initial_state_option == 'polvani_2007') then
    call polvani_2007(triang_trunc, vert_difference_option, pk, bk, &
                 vors(:,:,:,1), divs(:,:,:,1), ts(:,:,:,1), ln_ps(:,:,1), ug(:,:,:,1),  vg(:,:,:,1), &
                 tg(:,:,:,1), psg(:,:,1), vorg, divg, surf_geopotential)
    call get_polvani_2007_tracers(grid_tracers(:,:,:,1,:))
  else if(initial_state_option == 'jablonowski_2006') then
    call jablonowski_2006(reference_sea_level_press, triang_trunc, vert_coord_option, vert_difference_option, &
                 scale_heights, surf_res, p_press, p_sigma, exponent, pk, bk, &
                 vors(:,:,:,1), divs(:,:,:,1), ts(:,:,:,1), ln_ps(:,:,1), ug(:,:,:,1),  vg(:,:,:,1), &
                 tg(:,:,:,1), psg(:,:,1), vorg, divg, surf_geopotential)
  else
    call error_mesg('spectral_dynamics_init', trim(initial_state_option)//' is not a valid value of initial_state_option', FATAL)
  end if

  vors (:,:,:,2) = vors (:,:,:,1)
  divs (:,:,:,2) = divs (:,:,:,1)
  ts   (:,:,:,2) = ts   (:,:,:,1)
  ln_ps(:,:,  2) = ln_ps(:,:,  1)
  ug   (:,:,:,2) = ug   (:,:,:,1)
  vg   (:,:,:,2) = vg   (:,:,:,1)
  tg   (:,:,:,2) = tg   (:,:,:,1)
  psg  (:,:,  2) = psg  (:,:,  1)
  grid_tracers(:,:,:,2,:) = grid_tracers(:,:,:,1,:)

  do ntr = 1,num_tracers
    call trans_grid_to_spherical(grid_tracers(:,:,:,1,ntr), spec_tracers(:,:,:,1,ntr))
    spec_tracers(:,:,:,2,ntr) = spec_tracers(:,:,:,1,ntr)
  enddo
endif

return
end subroutine read_restart_or_do_coldstart
!===============================================================================================

subroutine allocate_fields

allocate (psg    (is:ie, js:je,             num_time_levels))
allocate (ug     (is:ie, js:je, num_levels, num_time_levels))
allocate (vg     (is:ie, js:je, num_levels, num_time_levels))
allocate (tg     (is:ie, js:je, num_levels, num_time_levels))

allocate (ln_ps(ms:me, ns:ne,             num_time_levels))
allocate ( vors(ms:me, ns:ne, num_levels, num_time_levels))
allocate ( divs(ms:me, ns:ne, num_levels, num_time_levels))
allocate (   ts(ms:me, ns:ne, num_levels, num_time_levels))

allocate (pk(num_levels+1), bk(num_levels+1))

allocate (vorg(is:ie, js:je, num_levels))
allocate (divg(is:ie, js:je, num_levels))

allocate (surf_geopotential(is:ie, js:je))

allocate (grid_tracers(is:ie, js:je, num_levels, num_time_levels, num_tracers))
allocate (spec_tracers(ms:me, ns:ne, num_levels, num_time_levels, num_tracers))

! Filling allocatable arrays with zeros immediately after allocation facilitates code debugging
psg=0.; ug=0.; vg=0.; tg=0.
ln_ps=cmplx(0.,0.); vors=cmplx(0.,0.); divs=cmplx(0.,0.); ts=cmplx(0.,0.); spec_tracers=cmplx(0.,0.)
pk=0.; bk=0.; vorg=0.; divg=0.; surf_geopotential=0.; grid_tracers=0.

return
end subroutine allocate_fields
!===============================================================================================
subroutine check_dynamics_nml

character(len=8)  :: ch_tmp1
character(len=16) :: ch_tmp2
integer :: itmp

!  Check for invalid values and incompatible combinations of namelist variables

if(num_fourier <= 0) then
  write(ch_tmp1,'(i8)') num_fourier
  call error_mesg('check_dynamics_nml',ch_tmp1//'is an invalid value for num_fourier.', FATAL)
endif

if(num_spherical <= 0) then
  write(ch_tmp1,'(i8)') num_spherical
  call error_mesg('check_dynamics_nml',ch_tmp1//'is an invalid value for num_spherical.', FATAL)
endif

if(fourier_inc <= 0) then
  write(ch_tmp1,'(i8)') fourier_inc
  call error_mesg('check_dynamics_nml',ch_tmp1//'is an invalid value for fourier_inc.', FATAL)
endif

if(num_levels <= 0) then
  write(ch_tmp1,'(i8)') num_levels
  call error_mesg('check_dynamics_nml',ch_tmp1//'is an invalid value for num_levels.', FATAL)
endif

if(lon_max < 3*num_fourier+1) then
  write(ch_tmp1,'(i8)') lon_max
  write(ch_tmp2,'(i8)') num_fourier
  call error_mesg('check_dynamics_nml','number of longitude points is too small for number of fourier waves.&
                  & lon_max='//ch_tmp1//'  num_fourier='//ch_tmp2, FATAL)
endif

if (triang_trunc) then
  itmp = 3
else
  itmp = 5
endif
if(2*lat_max < itmp*(num_spherical-1)+1) then
  write(ch_tmp1,'(i8)') lat_max
  write(ch_tmp2,'(i8)') num_spherical
  call error_mesg('check_dynamics_nml','number of latitude points is too small for number of meridional waves.&
                  &  lat_max='//ch_tmp1//'  num_spherical='//ch_tmp2, FATAL)
endif

if(damping_order < 0) then
  write(ch_tmp1,'(i8)') damping_order
  call error_mesg('check_dynamics_nml',ch_tmp1//' is an invalid value for damping_order.', FATAL)
endif

if(damping_order_vor < 0) then
  write(ch_tmp1,'(i8)') damping_order_vor
  call error_mesg('check_dynamics_nml',ch_tmp1//' is an invalid value for damping_order_vor.', FATAL)
endif

if(damping_order_div < 0) then
  write(ch_tmp1,'(i8)') damping_order_div
  call error_mesg('check_dynamics_nml',ch_tmp1//' is an invalid value for damping_order_div.', FATAL)
endif

if(damping_coeff < 0.) then
  write(ch_tmp2,'(e16.8)') damping_coeff
  call error_mesg('check_dynamics_nml',ch_tmp2//' is an invalid value for damping_coeff.', FATAL)
endif

if(damping_coeff_vor < 0.) then
  write(ch_tmp2,'(e16.8)') damping_coeff_vor
  call error_mesg('check_dynamics_nml',ch_tmp2//' is an invalid value for damping_coeff_vor.', FATAL)
endif

if(damping_coeff_div < 0.) then
  write(ch_tmp2,'(e16.8)') damping_coeff_div
  call error_mesg('check_dynamics_nml',ch_tmp2//' is an invalid value for damping_coeff_div.', FATAL)
endif

if(robert_coeff < 0. .or. robert_coeff > 1.) then
  write(ch_tmp2,'(1pe16.8)') robert_coeff
  call error_mesg('check_dynamics_nml',ch_tmp2//' is an invalid value for robert_coeff.', FATAL)
endif

if((do_energy_correction .or. do_water_correction) .and. .not.do_mass_correction) then
  call error_mesg('check_dynamics_nml','.not.do_mass_correction must be .true. when either & 
           &do_energy_correction or do_water_correction is .true.', FATAL)
endif

return

end subroutine check_dynamics_nml
!===============================================================================================
subroutine get_initial_fields(ug_out, vg_out, tg_out, psg_out, grid_tracers_out)
real, intent(out), dimension(:,:,:)   :: ug_out, vg_out, tg_out
real, intent(out), dimension(:,:)     :: psg_out
real, intent(out), dimension(:,:,:,:) :: grid_tracers_out

if(.not.module_is_initialized) then
  call error_mesg('get_initial_fields','dynamics has not been initialized',FATAL)
endif

if(previous /= 1 .or. current /= 1) then
  call error_mesg('get_initial_fields','This routine may be called only to get the&
                  & initial values after a cold_start',FATAL)
endif

ug_out  =  ug(:,:,:,1)
vg_out  =  vg(:,:,:,1)
tg_out  =  tg(:,:,:,1)
psg_out = psg(:,:,  1)
grid_tracers_out = grid_tracers(:,:,:,1,:)

end subroutine get_initial_fields
!===============================================================================================

subroutine spectral_dynamics(Time, psg_final, ug_final, vg_final, tg_final, tracer_attributes, grid_tracers_final, &
                             time_level_out, dt_psg, dt_ug, dt_vg, dt_tg, dt_tracers, wg_full, p_full, p_half, z_full)

type(time_type),  intent(in) :: Time
real, intent(out), dimension(is:, js:      ) :: psg_final
real, intent(out), dimension(is:, js:, :   ) :: ug_final, vg_final, tg_final
real, intent(out), dimension(is:, js:, :,:,:) :: grid_tracers_final
type(tracer_type),intent(inout), dimension(:) :: tracer_attributes
integer, intent(in)                           :: time_level_out

real, intent(inout), dimension(is:, js:      ) :: dt_psg
real, intent(inout), dimension(is:, js:, :   ) :: dt_ug, dt_vg, dt_tg
real, intent(inout), dimension(is:, js:, :, :) :: dt_tracers
real, intent(out),   dimension(is:, js:, :   ) :: wg_full, p_full
real, intent(out),   dimension(is:, js:, :   ) :: p_half
real, intent(in),    dimension(is:, js:, :   ) :: z_full

! < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < >

type(time_type) :: Time_diag

complex, dimension(ms:me, ns:ne              ) :: dt_ln_ps
complex, dimension(ms:me, ns:ne, num_levels  ) :: dt_vors, dt_divs, dt_ts, phis_plus_ke
real   , dimension(is:ie, js:je, num_levels  ) :: virtual_t, dp, dt_grid_tmp
real   , dimension(is:ie, js:je              ) :: dx_psg, dy_psg, ln_psg, dt_ln_psg
real   , dimension(is:ie, js:je, num_levels  ) :: phig_full, ln_p_full, phig_full_plus_ke
real   , dimension(is:ie, js:je, num_levels+1) :: phig_half, ln_p_half, wg

real, dimension(is:ie, js:je                         ) :: dt_psg_tmp
real, dimension(is:ie, js:je, num_levels             ) :: dt_ug_tmp, dt_vg_tmp, dt_tg_tmp
real, dimension(is:ie, js:je, num_levels, num_tracers) :: dt_tracers_tmp

integer :: j, k, time_level, seconds, days
real    :: delta_t

! < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < > < >

if(.not.module_is_initialized) then
  call error_mesg('spectral_dynamics','dynamics has not been initialized ', FATAL)
endif

step_loop: do step_number=1,num_steps

if(previous == current) then
  delta_t = dt_real/num_steps
else
  delta_t = 2*dt_real/num_steps
endif
if(num_time_levels == 2) then
  future = 3 - current
else
  call error_mesg('spectral_dynamics','Do not know how to set time pointers when num_time_levels does not equal 2',FATAL)
endif

dt_psg_tmp = dt_psg
dt_ug_tmp  = dt_ug
dt_vg_tmp  = dt_vg
dt_tg_tmp  = dt_tg
dt_tracers_tmp = dt_tracers

call initialize_corrections(dt_ug, dt_vg, dt_tg, dt_tracers, delta_t)

call pressure_variables (p_half, ln_p_half, p_full, ln_p_full, psg(:,:,current))

call compute_pressure_gradient  (ln_ps(:,:,current), psg(:,:,current), dx_psg, dy_psg)

if (use_virtual_temperature .and. .not.dry_model) then
  virtual_t = tg(:,:,:,current)*(1.0 + virtual_factor*grid_tracers(:,:,:,current,nhum))
else
  virtual_t = tg(:,:,:,current)
endif

call four_in_one (divg, ug(:,:,:,current), vg(:,:,:,current), virtual_t, psg(:,:,current), &
   ln_p_half, ln_p_full, p_full, dx_psg, dy_psg, dt_psg_tmp, wg, wg_full, dt_tg_tmp, dt_ug_tmp, dt_vg_tmp)

if(dry_model) then
  call compute_geopotential(tg(:,:,:,current), ln_p_half, ln_p_full, surf_geopotential, phig_full, phig_half)
else
  call compute_geopotential(tg(:,:,:,current), ln_p_half, ln_p_full, surf_geopotential, phig_full, phig_half, &
                            grid_tracers(:,:,:,current,nhum))
endif

dt_ln_psg = dt_psg_tmp/psg(:,:,current)
call trans_grid_to_spherical(dt_ln_psg, dt_ln_ps)

dp = p_half(:,:,2:num_levels+1) - p_half(:,:,1:num_levels)

if(uv_vert_advect_scheme == SECOND_CENTERED .or.  uv_vert_advect_scheme == FOURTH_CENTERED)         time_level=current
if(uv_vert_advect_scheme == VAN_LEER_LINEAR .or.  uv_vert_advect_scheme == FINITE_VOLUME_PARABOLIC) time_level=previous
call vert_advection(delta_t, wg, dp, ug(:,:,:,time_level), dt_grid_tmp, scheme=uv_vert_advect_scheme, form=ADVECTIVE_FORM)
dt_ug_tmp = dt_ug_tmp + dt_grid_tmp
call vert_advection(delta_t, wg, dp, vg(:,:,:,time_level), dt_grid_tmp, scheme=uv_vert_advect_scheme, form=ADVECTIVE_FORM)
dt_vg_tmp = dt_vg_tmp + dt_grid_tmp

if(t_vert_advect_scheme == SECOND_CENTERED .or.  t_vert_advect_scheme == FOURTH_CENTERED)         time_level=current
if(t_vert_advect_scheme == VAN_LEER_LINEAR .or.  t_vert_advect_scheme == FINITE_VOLUME_PARABOLIC) time_level=previous
call vert_advection(delta_t, wg, dp, tg(:,:,:,time_level),  dt_grid_tmp, scheme=t_vert_advect_scheme, form=ADVECTIVE_FORM)
dt_tg_tmp = dt_tg_tmp + dt_grid_tmp

call horizontal_advection(ts(:,:,:,current), ug(:,:,:,current), vg(:,:,:,current), dt_tg_tmp)
call trans_grid_to_spherical(dt_tg_tmp, dt_ts)

do k=1,num_levels
  do j = js,je
    dt_ug_tmp(:,j,k) = dt_ug_tmp(:,j,k) + (vorg(:,j,k) + coriolis(j))*vg(:,j,k,current)
    dt_vg_tmp(:,j,k) = dt_vg_tmp(:,j,k) - (vorg(:,j,k) + coriolis(j))*ug(:,j,k,current)
  enddo
enddo

call vor_div_from_uv_grid(dt_ug_tmp, dt_vg_tmp, dt_vors, dt_divs, triang = triang_trunc)

phig_full_plus_ke = phig_full + .5*(ug(:,:,:,current)**2 + vg(:,:,:,current)**2)
call trans_grid_to_spherical(phig_full_plus_ke, phis_plus_ke)
dt_divs = dt_divs - compute_laplacian(phis_plus_ke)

if(use_implicit) call implicit_correction (dt_divs, dt_ts, dt_ln_ps, divs, ts, ln_ps, delta_t, previous, current)

call compute_spectral_damping_vor (vors(:,:,:,previous), dt_vors, delta_t)
call compute_spectral_damping_div (divs(:,:,:,previous), dt_divs, delta_t)
call compute_spectral_damping     (ts  (:,:,:,previous), dt_ts,   delta_t)

if(.not.robert_complete_for_fields) then
  call error_mesg('spectral_dynamics','robert_complete_for_fields should be .true.',FATAL)
endif
if(.not.robert_complete_for_tracers) then
  call error_mesg('spectral_dynamics','robert_complete_for_tracers should be .true.',FATAL)
endif

if(step_number == num_steps) then
  call leapfrog_2level_A(ln_ps, dt_ln_ps, previous, current, future, delta_t, robert_coeff)
  call leapfrog_2level_A(vors,  dt_vors,  previous, current, future, delta_t, robert_coeff)
  call leapfrog_2level_A(divs,  dt_divs,  previous, current, future, delta_t, robert_coeff)
  call leapfrog_2level_A(ts,    dt_ts,    previous, current, future, delta_t, robert_coeff)
  robert_complete_for_fields = .false.
else
  call leapfrog         (ln_ps, dt_ln_ps, previous, current, future, delta_t, robert_coeff)
  call leapfrog         (vors , dt_vors , previous, current, future, delta_t, robert_coeff)
  call leapfrog         (divs , dt_divs , previous, current, future, delta_t, robert_coeff)
  call leapfrog         (ts   , dt_ts   , previous, current, future, delta_t, robert_coeff)
  robert_complete_for_fields = .true.
endif

call trans_spherical_to_grid(divs(:,:,:,future), divg)
call trans_spherical_to_grid(vors(:,:,:,future), vorg)
call uv_grid_from_vor_div(vors(:,:,:,future), divs(:,:,:,future), ug(:,:,:,future), vg(:,:,:,future))
call trans_spherical_to_grid(ts   (:,:,:,future), tg(:,:,:,future))
call trans_spherical_to_grid(ln_ps(:,:,  future), ln_psg)
psg(:,:,future) = exp(ln_psg)

if(minval(tg(:,:,:,future)) < valid_range_t(1) .or. maxval(tg(:,:,:,future)) > valid_range_t(2)) then
  call error_mesg('spectral_dynamics','temperatures out of valid range', FATAL)
endif

call update_tracers(tracer_attributes, dt_tracers_tmp, wg, p_half, delta_t)

call compute_corrections(delta_t, tracer_attributes)

previous = current
current  = future

call get_time(Time, seconds, days)
seconds = seconds + step_number*int(dt_real/2)
Time_diag = set_time(seconds, days)
call every_step_diagnostics( &
     Time_diag, psg(:,:,current), ug(:,:,:,current), vg(:,:,:,current), tg(:,:,:,current), grid_tracers(:,:,:,:,:), current)

enddo step_loop

psg_final = psg(:,:,  current)
ug_final  =  ug(:,:,:,current)
vg_final  =  vg(:,:,:,current)
tg_final  =  tg(:,:,:,current)
grid_tracers_final(:,:,:,time_level_out,:) = grid_tracers(:,:,:,current,:)

return
end subroutine spectral_dynamics

!================================================================================

subroutine four_in_one(divg, u_grid, v_grid, t_grid, p_surf, ln_p_half, ln_p_full, p_full, &
                       dx_psg, dy_psg, dt_psg, wg, wg_full, dt_tg, dt_ug, dt_vg)

real, intent(in),    dimension(:,:,:) :: divg, u_grid, v_grid, t_grid, ln_p_full, p_full
real, intent(in),    dimension(:,:  ) :: p_surf, dx_psg, dy_psg
real, intent(inout), dimension(:,:  ) :: dt_psg
real, intent(out),   dimension(:,:,:) :: wg_full
real, intent(out),   dimension(:,:,:) :: wg
real, intent(in),    dimension(:,:,:) :: ln_p_half
real, intent(inout), dimension(:,:,:) :: dt_tg, dt_ug, dt_vg

!  wg is dimensioned (is:ie, js:je, num_levels+1)
!  wg(:,:,k) = downward mass flux/per unit area across the K+1/2
!  cell boundary. This is the "vertical velocity" in the hybrid coordinate system.
!  When vertical coordinate is pure sigma: wg = psg*d(sigma)/dt

real, dimension(is:ie, js:je) :: dp, dp_inv, dlog_1, dlog_2, dlog_3, dmean, dmean_tot
real, dimension(is:ie, js:je) :: x1, x2, x3, x4, x5, p_surf_inv

real :: kappa
integer :: k

kappa = rdgas/cp_air
   
dmean_tot = 0.

if(vert_difference_option == 'simmons_and_burridge') then
  do k = 1,num_levels
    dp = dpk(k) + dbk(k)*p_surf
    dp_inv = 1/dp
    dlog_1 = ln_p_half(:,:,k+1) - ln_p_full(:,:,k)
    dlog_2 = ln_p_full(:,:,k)   - ln_p_half(:,:,k)
    dlog_3 = ln_p_half(:,:,k+1) - ln_p_half(:,:,k)
    x1 = (bk(k+1)*dlog_1 + bk(k)*dlog_2)*dp_inv
    x2 = x1*dx_psg
    x3 = x1*dy_psg
    dt_ug(:,:,k) = dt_ug(:,:,k) - rdgas*t_grid(:,:,k)*x2
    dt_vg(:,:,k) = dt_vg(:,:,k) - rdgas*t_grid(:,:,k)*x3
    dmean = divg(:,:,k)*dp + dbk(k)*(u_grid(:,:,k)*dx_psg + v_grid(:,:,k)*dy_psg)
    x4 = (dmean_tot*dlog_3 + dmean*dlog_1)*dp_inv
    x5 = x4 - u_grid(:,:,k)*x2 - v_grid(:,:,k)*x3
    dt_tg(:,:,k) = dt_tg(:,:,k) - kappa*t_grid(:,:,k) * x5
    wg_full(:,:,k) = -x5*p_full(:,:,k)
    dmean_tot = dmean_tot + dmean
    wg(:,:,k+1) = - dmean_tot
  enddo
else if(vert_difference_option == 'mcm') then
  p_surf_inv = 1.0/p_surf
  do k = 1,num_levels
    dp = dpk(k) + dbk(k)*p_surf
    x2 = dx_psg*p_surf_inv
    x3 = dy_psg*p_surf_inv
    dt_ug(:,:,k) = dt_ug(:,:,k) - rdgas*t_grid(:,:,k)*x2
    dt_vg(:,:,k) = dt_vg(:,:,k) - rdgas*t_grid(:,:,k)*x3
    dmean = divg(:,:,k)*dp + dbk(k)*(u_grid(:,:,k)*dx_psg + v_grid(:,:,k)*dy_psg)
    x4 = (dmean_tot + 0.5*dmean)/p_full(:,:,k)
    x5 = x4 - u_grid(:,:,k)*x2 - v_grid(:,:,k)*x3
    dt_tg(:,:,k) = dt_tg(:,:,k) - kappa*t_grid(:,:,k) * x5
    wg_full(:,:,k) = -x5*p_full(:,:,k)
    dmean_tot = dmean_tot + dmean
    wg(:,:,k+1) = - dmean_tot
  enddo
endif

dt_psg = dt_psg - dmean_tot

do k = 2,num_levels
  wg(:,:,k) = wg(:,:,k) + dmean_tot*bk(k)
enddo

wg(:,:,1           ) = 0.0
wg(:,:,num_levels+1) = 0.0

return
end subroutine four_in_one

!================================================================================

subroutine update_tracers(tracer_attributes, dt_tr, wg, p_half, delta_t)

type(tracer_type), intent(inout), dimension(:) :: tracer_attributes
real   , intent(inout), dimension(:,:,:,:) :: dt_tr
real   , intent(in   ), dimension(:,:,:  ) :: wg, p_half
real   , intent(in   )  :: delta_t

complex, dimension(ms:me, ns:ne, num_levels) :: dt_trs
real,    dimension(is:ie, js:je, num_levels) :: dp, dt_tmp, tr_future
integer :: ntr, time_level

dp = p_half(:,:,2:num_levels+1) - p_half(:,:,1:num_levels)

do ntr = 1, num_tracers
  if(trim(tracer_attributes(ntr)%numerical_representation) == 'spectral') then
    call horizontal_advection(spec_tracers(:,:,:,current,ntr), ug(:,:,:,current), vg(:,:,:,current), dt_tr(:,:,:,ntr))
    if(tracer_vert_advect_scheme(ntr) == SECOND_CENTERED .or. &
       tracer_vert_advect_scheme(ntr) == FOURTH_CENTERED)         time_level=current
    if(tracer_vert_advect_scheme(ntr) == VAN_LEER_LINEAR .or. &
       tracer_vert_advect_scheme(ntr) == FINITE_VOLUME_PARABOLIC) time_level=previous
    call vert_advection(delta_t, wg, dp, grid_tracers(:,:,:,time_level,ntr), dt_tmp, &
                     scheme=tracer_vert_advect_scheme(ntr), form=ADVECTIVE_FORM)
    dt_tr(:,:,:,ntr) = dt_tr(:,:,:,ntr) + dt_tmp
    if(trim(tracer_attributes(ntr)%hole_filling) == 'on') then
      call water_borrowing (dt_tr(:,:,:,ntr), grid_tracers(:,:,:,previous,ntr), current, p_half, delta_t)
    endif
    call trans_grid_to_spherical  (dt_tr(:,:,:,ntr), dt_trs)
    call compute_spectral_damping (spec_tracers(:,:,:,previous,ntr), dt_trs, delta_t)
    if(step_number == num_steps) then
      call leapfrog_2level_A(spec_tracers(:,:,:,:,ntr),dt_trs,previous,current,future,delta_t,tracer_attributes(ntr)%robert_coeff)
      robert_complete_for_tracers = .false.
    else
      call leapfrog(spec_tracers(:,:,:,:,ntr),dt_trs,previous,current,future,delta_t,tracer_attributes(ntr)%robert_coeff)
      robert_complete_for_tracers = .true.
    endif
    call trans_spherical_to_grid  (spec_tracers(:,:,:,future,ntr), grid_tracers(:,:,:,future,ntr))
  else if(trim(tracer_attributes(ntr)%numerical_representation) == 'grid') then
    tr_future = grid_tracers(:,:,:,previous,ntr) + delta_t*dt_tr(:,:,:,ntr)
    dt_tr(:,:,:,ntr) = 0.0
    call a_grid_horiz_advection (ug(:,:,:,current), vg(:,:,:,current), tr_future, delta_t, dt_tr(:,:,:,ntr))
    tr_future = tr_future + delta_t*dt_tr(:,:,:,ntr)
    dp = p_half(:,:,2:num_levels+1) - p_half(:,:,1:num_levels)
    call vert_advection(delta_t, wg, dp, tr_future, dt_tmp, scheme=tracer_vert_advect_scheme(ntr), form=ADVECTIVE_FORM)
    tr_future = tr_future + delta_t*dt_tmp
    if(step_number == num_steps) then
      grid_tracers(:,:,:,current,ntr) = grid_tracers(:,:,:,current,ntr) + &
      tracer_attributes(ntr)%robert_coeff*(grid_tracers(:,:,:,previous,ntr) - 2.0*grid_tracers(:,:,:,current,ntr))
      robert_complete_for_tracers = .false.
    else
      grid_tracers(:,:,:,current,ntr) = grid_tracers(:,:,:,current,ntr) + &
      tracer_attributes(ntr)%robert_coeff*(grid_tracers(:,:,:,previous,ntr) - 2.0*grid_tracers(:,:,:,current,ntr) + tr_future)
      robert_complete_for_tracers = .true.
    endif
    grid_tracers(:,:,:,future,ntr) = tr_future
  else
    call error_mesg('update_tracers',trim(tracer_attributes(ntr)%numerical_representation)// &
           ' is an invalid numerical_representation', FATAL)
  endif
enddo

return
end subroutine update_tracers

!=================================================================================================

subroutine compute_pressure_gradient(ln_ps, psg, dx_psg, dy_psg)

complex, intent(in ), dimension(:,:) :: ln_ps
real   , intent(in ), dimension(:,:) :: psg
real   , intent(out), dimension(:,:) :: dx_psg, dy_psg

complex, dimension(ms:me, ns:ne) :: dx_ln_ps, dy_ln_ps

call compute_gradient_cos(ln_ps, dx_ln_ps, dy_ln_ps)
call trans_spherical_to_grid(dx_ln_ps, dx_psg)
call trans_spherical_to_grid(dy_ln_ps, dy_psg)
dx_psg = psg*dx_psg
dy_psg = psg*dy_psg
call divide_by_cos(dx_psg)
call divide_by_cos(dy_psg)

return
end subroutine compute_pressure_gradient 

!===================================================================================

subroutine compute_corrections(delta_t, tracer_attributes)

real,              intent(in   )               :: delta_t
type(tracer_type), intent(inout), dimension(:) :: tracer_attributes

real :: mass_correction_factor, temperature_correction, water_correction_factor
real :: mean_surf_press_tmp,    mean_energy_tmp,        mean_water_tmp

if(do_mass_correction) then
  mean_surf_press_tmp = area_weighted_global_mean(psg(:,:,future))
  mass_correction_factor = mean_surf_press_previous/mean_surf_press_tmp
  psg(:,:,future) = mass_correction_factor*psg(:,:,future)
  if(ms == 0 .and. ns == 0) then
    ln_ps(0,0,future) = ln_ps(0,0,future) + sqrt(2.)*log(mass_correction_factor)
  endif
endif

if(do_energy_correction) then
  mean_energy_tmp = mass_weighted_global_integral( &
          0.5*(ug(:,:,:,future)**2 + vg(:,:,:,future)**2) + cp_air*tg(:,:,:,future), psg(:,:,future))
  temperature_correction = grav*(mean_energy_previous - mean_energy_tmp)/(cp_air*mean_surf_press_previous)
  tg(:,:,:,future) = tg(:,:,:,future) + temperature_correction
  if(ms == 0 .and. ns == 0) then
    ts(0,0,:,future) = ts(0,0,:,future) + sqrt(2.)*temperature_correction
  endif
endif

if(do_water_correction) then
  if(dry_model) then
    call error_mesg('compute_corrections','do_water_correction must be .false. in a dry model (default is .true.)', FATAL)
  else
    mean_water_tmp  = mass_weighted_global_integral(grid_tracers(:,:,:,future,nhum), psg(:,:,future))
    if(mean_water_tmp > 0.) then
      water_correction_factor = mean_water_previous/mean_water_tmp
      grid_tracers(:,:,:,future,nhum) = water_correction_factor*grid_tracers(:,:,:,future,nhum)
      if(tracer_attributes(nhum)%numerical_representation == 'spectral') then
        spec_tracers(:,:,:,future,nhum) = water_correction_factor*spec_tracers(:,:,:,future,nhum)
      endif
    endif
  endif
endif

return
end subroutine compute_corrections 

!===================================================================================

subroutine initialize_corrections(dt_ug, dt_vg, dt_tg, dt_tracers, delta_t)

real,    intent(in), dimension(:,:,:)   :: dt_ug, dt_vg, dt_tg
real,    intent(in), dimension(:,:,:,:) :: dt_tracers
real,    intent(in) :: delta_t

real, dimension(is:ie, js:je, num_levels) :: energy

if(do_mass_correction) then
  mean_surf_press_previous = area_weighted_global_mean(psg(:,:,previous))
endif

if(do_energy_correction) then
   energy  =  0.5*((ug(:,:,:,previous)+dt_ug*delta_t)**2 + (vg(:,:,:,previous)+dt_vg*delta_t)**2)   &
               +cp_air*(tg(:,:,:,previous)+dt_tg*delta_t)
   mean_energy_previous = mass_weighted_global_integral( energy, psg(:,:,previous))

   energy  =  0.5*((ug(:,:,:,current)+dt_ug*delta_t)**2 + (vg(:,:,:,current)+dt_vg*delta_t)**2)   &
               +cp_air*(tg(:,:,:,current)+dt_tg*delta_t)
endif

if(do_water_correction) then
  if(dry_model) then
    call error_mesg('initialize_corrections','do_water_correction must be .false. in a dry model&
                     & (default is .true.)', FATAL)
  else
    mean_water_previous = &
      mass_weighted_global_integral(grid_tracers(:,:,:,previous,nhum) + delta_t*dt_tracers(:,:,:,nhum), psg(:,:,previous))
  endif
endif

return
end subroutine initialize_corrections 

!================================================================================

subroutine get_surf_geopotential(surf_geopotential_out)
real, intent(out), dimension(:,:) :: surf_geopotential_out
character(len=64) :: chtmp='shape(surf_geopotential)=              should be                '

if(.not.module_is_initialized) then
  call error_mesg('get_surf_geopotential', 'spectral_dynamics_init has not been called.', FATAL)
endif

if(any(shape(surf_geopotential_out) /= shape(surf_geopotential))) then
  write(chtmp(26:37),'(3i4)') shape(surf_geopotential_out)
  write(chtmp(50:61),'(3i4)') shape(surf_geopotential)
  call error_mesg('get_surf_geopotential', 'surf_geopotential has wrong shape. '//chtmp, FATAL)
endif

surf_geopotential_out = surf_geopotential

return
end subroutine get_surf_geopotential
!================================================================================
subroutine get_reference_sea_level_press(reference_sea_level_press_out)
real, intent(out) :: reference_sea_level_press_out

if(.not.module_is_initialized) then
  call error_mesg('get_reference_sea_level_press', 'spectral_dynamics_init has not been called.', FATAL)
endif

reference_sea_level_press_out = reference_sea_level_press

return
end subroutine get_reference_sea_level_press
!================================================================================
subroutine get_use_virtual_temperature(use_virtual_temperature_out)
logical, intent(out) :: use_virtual_temperature_out

if(.not.module_is_initialized) then
  call error_mesg('get_use_virtual_temperature', 'spectral_dynamics_init has not been called.', FATAL)
endif

use_virtual_temperature_out = use_virtual_temperature

return
end subroutine get_use_virtual_temperature
!================================================================================
subroutine get_num_levels(num_levels_out)
integer, intent(out) :: num_levels_out

if(.not.module_is_initialized) then
  call error_mesg('get_num_levels', 'spectral_dynamics_init has not been called.', FATAL)
endif

num_levels_out = num_levels

return
end subroutine get_num_levels
!================================================================================
subroutine get_pk_bk(pk_out, bk_out)
real, intent(out), dimension(:) :: pk_out, bk_out
character(len=32) :: chtmp='size(pk)=      size(bk)=        '

if(.not.module_is_initialized) then
  call error_mesg('get_pk_bk', 'spectral_dynamics_init has not been called.', FATAL)
endif
if(size(pk_out,1) /= size(bk_out,1)) then
  write(chtmp(10:13),'(i4)') size(pk_out,1)
  write(chtmp(25:28),'(i4)') size(bk_out,1)
  call error_mesg('get_pk_bk', 'size(pk) is not equal to size(bk). '//chtmp, FATAL)
endif

pk_out = pk
bk_out = bk

return
end subroutine get_pk_bk
!================================================================================
subroutine complete_update_of_future(psg_in, ug_in, vg_in, tg_in, tracer_attributes, grid_tracers_in)

real,              intent(in), dimension(:,:    ) :: psg_in
real,              intent(in), dimension(:,:,:  ) :: ug_in, vg_in, tg_in
type(tracer_type), intent(in), dimension(:      ) :: tracer_attributes
real,              intent(in), dimension(:,:,:,:) :: grid_tracers_in
real, dimension(size(psg,1), size(psg,2)) :: ln_psg
integer :: ntr

! The time level pointers may be confusing here.
! The future level of the fields are passed to this routine in atmosphere_up,
! and they are loaded into the current level here.
! The reason for this is that the time level pointers in atmosphere_mod have
! not yet been updated, but they have been updated for spectral_dyanmics_mod
! (in Subroutine spectral_dynamics.) The result is that future in
! atmosphere_up points to the same time level as current in this routine.

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

psg(:,:,  current) = psg_in
ug (:,:,:,current) = ug_in
vg (:,:,:,current) = vg_in
tg (:,:,:,current) = tg_in
grid_tracers(:,:,:,current,:) = grid_tracers_in

call vor_div_from_uv_grid(ug(:,:,:,current), vg(:,:,:,current), vors(:,:,:,current), divs(:,:,:,current), triang=triang_trunc)
call trans_spherical_to_grid(vors(:,:,:,current), vorg)
call trans_spherical_to_grid(divs(:,:,:,current), divg)
ln_psg = alog(psg(:,:,current))
call trans_grid_to_spherical(ln_psg, ln_ps(:,:,current))
call trans_grid_to_spherical(tg_in,  ts(:,:,:,current))
do ntr=1,num_tracers
  if(uppercase(trim(tracer_attributes(ntr)%numerical_representation)) == 'SPECTRAL') then
    call trans_grid_to_spherical(grid_tracers(:,:,:,current,ntr), spec_tracers(:,:,:,current,ntr))
  endif
enddo

return
end subroutine complete_update_of_future
!================================================================================
subroutine complete_robert_filter(tracer_attributes)
type(tracer_type), intent(inout), dimension(:) :: tracer_attributes
integer :: ntr

if(robert_complete_for_fields) then
  call error_mesg('complete_robert_filter','This routine should not be called when robert_complete_for_fields=.true.',FATAL)
endif
call leapfrog_2level_B(ln_ps, previous, current, robert_coeff)
call leapfrog_2level_B(vors,  previous, current, robert_coeff)
call leapfrog_2level_B(divs,  previous, current, robert_coeff)
call leapfrog_2level_B(ts,    previous, current, robert_coeff)
robert_complete_for_fields=.true.

if(num_tracers > 0 .and. robert_complete_for_tracers) then
  call error_mesg('complete_robert_filter','This routine should not be called when robert_complete_for_tracers=.true.',FATAL)
endif

do ntr = 1, num_tracers
  if(uppercase(trim(tracer_attributes(ntr)%numerical_representation)) == 'SPECTRAL') then
    call leapfrog_2level_B(spec_tracers(:,:,:,:,ntr), previous, current, tracer_attributes(ntr)%robert_coeff)
  else
    call leapfrog_2level_B(grid_tracers(:,:,:,:,ntr), previous, current, tracer_attributes(ntr)%robert_coeff)
  endif
  robert_complete_for_tracers=.true.
enddo

return
end subroutine complete_robert_filter
!================================================================================

subroutine spectral_dynamics_end(tracer_attributes, Time)

type(tracer_type), intent(in), dimension(:) :: tracer_attributes
type(time_type), intent(in), optional :: Time
integer :: ntr, nt
character(len=64) :: file, tr_name

if(.not.module_is_initialized) return

file='RESTART/spectral_dynamics.res'
call write_data(trim(file), 'previous', previous, no_domain=.true.)
call write_data(trim(file), 'current',  current,  no_domain=.true.)
call write_data(trim(file), 'pk', pk, no_domain=.true.)
call write_data(trim(file), 'bk', bk, no_domain=.true.)
do nt=1,num_time_levels
  call write_data(trim(file), 'vors_real',   real(vors (:,:,:,nt)), spectral_domain)
  call write_data(trim(file), 'vors_imag',  aimag(vors (:,:,:,nt)), spectral_domain)
  call write_data(trim(file), 'divs_real',   real(divs (:,:,:,nt)), spectral_domain)
  call write_data(trim(file), 'divs_imag',  aimag(divs (:,:,:,nt)), spectral_domain)
  call write_data(trim(file), 'ts_real',     real(ts   (:,:,:,nt)), spectral_domain)
  call write_data(trim(file), 'ts_imag',    aimag(ts   (:,:,:,nt)), spectral_domain)
  call write_data(trim(file), 'ln_ps_real',  real(ln_ps(:,:,  nt)), spectral_domain)
  call write_data(trim(file), 'ln_ps_imag', aimag(ln_ps(:,:,  nt)), spectral_domain)
  call write_data(trim(file), 'ug',   ug(:,:,:,nt), grid_domain)
  call write_data(trim(file), 'vg',   vg(:,:,:,nt), grid_domain)
  call write_data(trim(file), 'tg',   tg(:,:,:,nt), grid_domain)
  call write_data(trim(file), 'psg', psg(:,:,  nt), grid_domain)
  do ntr = 1,num_tracers
    tr_name = trim(tracer_attributes(ntr)%name)
    call write_data(trim(file), trim(tr_name), grid_tracers(:,:,:,nt,ntr), grid_domain)
    if(uppercase(trim(tracer_attributes(ntr)%numerical_representation)) == 'SPECTRAL') then
      call write_data(trim(file), trim(tr_name)//'_real',  real(spec_tracers(:,:,:,nt,ntr)), spectral_domain)
      call write_data(trim(file), trim(tr_name)//'_imag', aimag(spec_tracers(:,:,:,nt,ntr)), spectral_domain)
    endif
  enddo
enddo
call write_data(trim(file), 'vorg', vorg, grid_domain)
call write_data(trim(file), 'divg', divg, grid_domain)
call write_data(trim(file), 'surf_geopotential', surf_geopotential, grid_domain)

deallocate(ug, vg, tg, psg)
deallocate(sin_lat, coriolis)
deallocate(pk, bk, dpk, dbk)
deallocate(ln_ps, vors, divs, ts)
deallocate(vorg, divg)
deallocate(surf_geopotential)
deallocate(spec_tracers, grid_tracers)

if(use_implicit) call implicit_end
call spectral_damping_end
call fv_advection_end
call every_step_diagnostics_end(Time)
call spectral_diagnostics_end
call press_and_geopot_end
call transforms_end
call set_domain(grid_domain)
module_is_initialized = .false.

return
end subroutine spectral_dynamics_end
!===================================================================================
subroutine spectral_diagnostics_init(Time)

type(time_type), intent(in) :: Time
real, dimension(lon_max  ) :: lon
real, dimension(lon_max+1) :: lonb
real, dimension(lat_max  ) :: lat
real, dimension(lat_max+1) :: latb
real, dimension(num_levels)   :: p_full, ln_p_full
real, dimension(num_levels+1) :: p_half, ln_p_half
integer, dimension(3) :: axes_3d_half, axes_3d_full
integer :: id_lonb, id_latb, id_phalf, id_lon, id_lat, id_pfull
integer :: id_pk, id_bk, id_zsurf, ntr
real :: rad_to_deg
logical :: used
real,dimension(2) :: vrange
character(len=128) :: tname, longname, units

vrange = (/ -400., 400. /)

rad_to_deg = 180./pi
call get_grid_boundaries(lonb,latb,global=.true.)
call get_deg_lon(lon)
call get_deg_lat(lat)

id_lonb=diag_axis_init('lonb', rad_to_deg*lonb, 'degrees_E', 'x', 'longitude edges', set_name=mod_name, Domain2=grid_domain)
id_latb=diag_axis_init('latb', rad_to_deg*latb, 'degrees_N', 'y', 'latitude edges',  set_name=mod_name, Domain2=grid_domain)
id_lon =diag_axis_init('lon', lon, 'degrees_E', 'x', 'longitude', set_name=mod_name, Domain2=grid_domain, edges=id_lonb)
id_lat =diag_axis_init('lat', lat, 'degrees_N', 'y', 'latitude',  set_name=mod_name, Domain2=grid_domain, edges=id_latb)

call pressure_variables(p_half, ln_p_half, p_full, ln_p_full, reference_sea_level_press)
p_half = .01*p_half
p_full = .01*p_full
id_phalf = diag_axis_init('phalf',p_half,'hPa','z','approx half pressure level',direction=-1,set_name=mod_name)
id_pfull = diag_axis_init('pfull',p_full,'hPa','z','approx full pressure level',direction=-1,set_name=mod_name,edges=id_phalf)

axes_3d_half = (/ id_lon, id_lat, id_phalf /)
axes_3d_full = (/ id_lon, id_lat, id_pfull /)
axis_id(1) = id_lon
axis_id(2) = id_lat
axis_id(3) = id_pfull
axis_id(4) = id_phalf

id_pk = register_static_field(mod_name, 'pk', (/id_phalf/), 'vertical coordinate pressure values', 'pascals')
id_bk = register_static_field(mod_name, 'bk', (/id_phalf/), 'vertical coordinate sigma values', 'none')
id_zsurf = register_static_field(mod_name, 'zsurf', (/id_lon,id_lat/), 'geopotential height at the surface', 'm')

if(id_pk    > 0) used = send_data(id_pk, pk, Time)
if(id_bk    > 0) used = send_data(id_bk, bk, Time)
if(id_zsurf > 0) used = send_data(id_zsurf, surf_geopotential/grav, Time)

id_ps  = register_diag_field(mod_name, &
      'ps', (/id_lon,id_lat/),       Time, 'surface pressure',             'pascals')

id_u   = register_diag_field(mod_name, &
      'ucomp',   axes_3d_full,       Time, 'zonal wind component',         'm/sec',      range=vrange)

id_v   = register_diag_field(mod_name, &
      'vcomp',   axes_3d_full,       Time, 'meridional wind component',    'm/sec',      range=vrange)

id_uu  = register_diag_field(mod_name, &
      'ucomp_sq',axes_3d_full,       Time, 'zonal wind squared',           '(m/sec)**2', range=(/0.,vrange(2)**2/))

id_vv  = register_diag_field(mod_name, &
      'vcomp_sq',axes_3d_full,       Time, 'meridional wind squared',      '(m/sec)**2', range=(/0.,vrange(2)**2/))

id_uv  = register_diag_field(mod_name, &
      'ucomp_vcomp', axes_3d_full, Time, 'zonal wind * meridional wind', '(m/sec)**2', range=(/-vrange(2)**2,vrange(2)**2/))

id_omega_t = register_diag_field(mod_name, &
      'omega_temp',axes_3d_full,     Time, 'dp/dt * temperature',          'Pa*K/sec')

id_wspd= register_diag_field(mod_name, &
      'wspd',    axes_3d_full,       Time, 'wind speed',                   'm/sec',      range=(/0.,vrange(2)/))

id_t   = register_diag_field(mod_name, &
      'temp',    axes_3d_full,       Time, 'temperature',                  'deg_k',      range=valid_range_t)

id_tt  = register_diag_field(mod_name, &
      'temp_sq', axes_3d_full,       Time, 'temperature squared',          'deg_k**2',   range=(/0.,valid_range_t(2)**2/))

id_vor = register_diag_field(mod_name, &
      'vor',     axes_3d_full,       Time, 'vorticity',                    'sec**-1')

id_div = register_diag_field(mod_name, &
      'div',     axes_3d_full,       Time, 'divergence',                   'sec**-1')

id_omega  = register_diag_field(mod_name, &
      'omega',   axes_3d_full,       Time, 'dp/dt vertical velocity',      'Pa/sec')

id_omega_omega = register_diag_field(mod_name, &
      'omega_sq',axes_3d_full,       Time, 'omega squared',                '(Pa/sec)**2')

id_pres_full = register_diag_field(mod_name, &
      'pres_full',    axes_3d_full,       Time, 'pressure at full model levels', 'pascals')

id_pres_half = register_diag_field(mod_name, &
      'pres_half',    axes_3d_half,       Time, 'pressure at half model levels', 'pascals')

id_zfull   = register_diag_field(mod_name, &
      'height',  axes_3d_full,       Time, 'geopotential height at full model levels','m')

id_zhalf   = register_diag_field(mod_name, &
      'height_half',  axes_3d_half,  Time, 'geopotential height at half model levels','m')

id_slp = register_diag_field(mod_name, &
      'slp',(/id_lon,id_lat/),       Time, 'sea level pressure',           'pascals')

if(id_slp > 0) then
  gamma = 0.006
  expf = rdgas*gamma/grav
  expf_inverse = 1./expf
endif

allocate(id_tr(num_tracers))
do ntr=1,num_tracers
  call get_tracer_names(MODEL_ATMOS, ntr, tname, longname, units)
  id_tr(ntr) = register_diag_field(mod_name, tname, axes_3d_full, Time, longname, units)
enddo

id_vort_norm = register_diag_field(mod_name, 'vort_norm', Time, 'vorticity norm', '1/(m*sec)')
id_EKE       = register_diag_field(mod_name, 'EKE', Time, 'eddy kinetic energy', 'J/m^2')

return
end subroutine spectral_diagnostics_init
!===================================================================================
subroutine spectral_diagnostics(Time, p_surf, u_grid, v_grid, t_grid, wg_full, tr_grid, time_level)

type(time_type), intent(in) :: Time
real, intent(in), dimension(is:, js:)          :: p_surf
real, intent(in), dimension(is:, js:, :)       :: u_grid, v_grid, t_grid, wg_full
real, intent(in), dimension(is:, js:, :, :, :) :: tr_grid
integer, intent(in) :: time_level

real, dimension(is:ie, js:je, num_levels)    :: ln_p_full, p_full, z_full, worka3d, workb3d
real, dimension(is:ie, js:je, num_levels+1)  :: ln_p_half, p_half, z_half
real, dimension(is:ie, js:je)                :: t_low, slp, worka2d, workb2d
complex, dimension(ms:me, ns:ne, num_levels) :: vor_spec, div_spec
complex, dimension(ms:me, ns:ne)             :: vorx, vory
logical :: used
integer :: ntr, i, j, k
character(len=8) :: err_msg_1, err_msg_2
real, dimension(lon_max, lat_max) :: work_global
real :: vort_norm, EKE

if(id_ps  > 0)    used = send_data(id_ps,  p_surf, Time)
if(id_u   > 0)    used = send_data(id_u,   u_grid, Time)
if(id_v   > 0)    used = send_data(id_v,   v_grid, Time)
if(id_t   > 0)    used = send_data(id_t,   t_grid, Time)
if(id_vor > 0)    used = send_data(id_vor, vorg, Time)
if(id_div > 0)    used = send_data(id_div, divg, Time)
if(id_omega > 0)  used = send_data(id_omega, wg_full, Time)

if(id_zfull > 0 .or. id_zhalf > 0) then
  call compute_pressures_and_heights(t_grid, p_surf, surf_geopotential, z_full, z_half, p_full, p_half)
else if(id_pres_half > 0 .or. id_pres_full > 0 .or. id_slp > 0) then
  call pressure_variables(p_half, ln_p_half, p_full, ln_p_full, p_surf)
endif

if(id_zfull > 0)   used = send_data(id_zfull,      z_full, Time)
if(id_zhalf > 0)   used = send_data(id_zhalf,      z_half, Time)
if(id_pres_full>0) used = send_data(id_pres_full,  p_full, Time)
if(id_pres_half>0) used = send_data(id_pres_half,  p_half, Time)

if(id_wspd > 0) then
  worka3d = sqrt(u_grid**2 + v_grid**2)
  used = send_data(id_wspd, worka3d, Time)
endif
if(id_uu > 0) then
  worka3d = u_grid**2
  used = send_data(id_uu, worka3d, Time)
endif
if(id_vv > 0) then
  worka3d = v_grid**2
  used = send_data(id_vv, worka3d, Time)
endif
if(id_uv > 0) then
  worka3d = u_grid*v_grid
  used = send_data(id_uv, worka3d, Time)
endif
if(id_tt > 0) then
  worka3d = t_grid**2
  used = send_data(id_tt, worka3d, Time)
endif
if(id_omega_omega > 0) then
  worka3d = wg_full*wg_full
  used = send_data(id_omega_omega, worka3d, Time)
endif
if(id_omega_t > 0) then
  worka3d = wg_full*t_grid
  used = send_data(id_omega_t, worka3d, Time)
endif

if(size(tr_grid,5) /= num_tracers) then
  write(err_msg_1,'(i8)') size(tr_grid,5)
  write(err_msg_2,'(i8)') num_tracers
  call error_mesg('spectral_diagnostics','size(tracers)='//err_msg_1//' Should be='//err_msg_2, FATAL)
endif
do ntr=1,num_tracers
  if(id_tr(ntr) > 0) used = send_data(id_tr(ntr), tr_grid(:,:,:,time_level,ntr), Time)
enddo

if(id_slp > 0) then
  do j=js,je
    do i=is,ie
      do k=1,num_levels
        if(p_full(i,j,k)/p_surf(i,j) > .8) goto 20
      enddo
      call error_mesg('spectral_diagnostics','No sigma values .gt. 0.8  Cannot compute slp',FATAL)
20    continue
      t_low(i,j) = t_grid(i,j,k)*(p_full(i,j,k)/p_surf(i,j))**(-expf)
      slp(i,j) = p_surf(i,j)*((t_low(i,j) + gamma*surf_geopotential(i,j)/grav)/t_low(i,j))**expf_inverse
    enddo
  enddo
  used = send_data(id_slp, slp, Time)
endif

if(interval_alarm(Time, Time_step, Alarm_time, Alarm_interval)) then
  call global_integrals(Time, p_surf, u_grid, v_grid, t_grid, wg_full, tr_grid(:,:,:,time_level,:))
endif

if(id_vort_norm > 0) then
  call vor_div_from_uv_grid(u_grid(:,:,num_levels), v_grid(:,:,num_levels), &
                            vor_spec(:,:,num_levels), div_spec(:,:,num_levels), triang=triang_trunc)
  call compute_gradient_cos(vor_spec(:,:,num_levels), vorx, vory)
  call trans_spherical_to_grid(vorx, worka2d)
  call trans_spherical_to_grid(vory, workb2d)
  call divide_by_cos(worka2d)
  call divide_by_cos(workb2d)
  worka2d = sqrt(worka2d**2 + workb2d**2)
  call mpp_global_field(grid_domain, worka2d, work_global)
  vort_norm = maxval(work_global)
  used = send_data(id_vort_norm, vort_norm, Time)
endif
if(id_EKE > 0) then
  call vor_div_from_uv_grid(u_grid, v_grid, vor_spec, div_spec, triang=triang_trunc)
  if(ms == 0) then
    vor_spec(0,:,:) = cmplx(0.0,0.0)
    div_spec(0,:,:) = cmplx(0.0,0.0)
  endif
  call uv_grid_from_vor_div(vor_spec, div_spec, worka3d, workb3d)
  EKE = mass_weighted_global_integral(.5*(worka3d**2 + workb3d**2), p_surf)
  used = send_data(id_EKE, EKE, Time)
endif

return
end subroutine spectral_diagnostics
!===================================================================================
subroutine global_integrals(Time, p_surf, u_grid, v_grid, t_grid, wg_full, tr_grid)
type(time_type), intent(in) :: Time
real, intent(in), dimension(is:ie, js:je)                          :: p_surf
real, intent(in), dimension(is:ie, js:je, num_levels)              :: u_grid, v_grid, t_grid, wg_full
real, intent(in), dimension(is:ie, js:je, num_levels, num_tracers) :: tr_grid
integer :: year, month, days, hours, minutes, seconds
character(len=4), dimension(12) :: month_name

month_name=(/' Jan',' Feb',' Mar',' Apr',' May',' Jun',' Jul',' Aug',' Sep',' Oct',' Nov',' Dec'/)

if(mpp_pe() == mpp_root_pe()) then
  if(get_calendar_type() == NO_CALENDAR) then
    call get_time(Time, seconds, days)
    write(*,100) days, seconds
  else
    call get_date(Time, year, month, days, hours, minutes, seconds)
    write(*,200) year, month_name(month), days, hours, minutes, seconds
  endif
endif
100 format(' Integration completed through',i6,' days',i6,' seconds')
200 format(' Integration completed through',i5,a4,i3,2x,i2,':',i2,':',i2)

end subroutine global_integrals
!===================================================================================
function get_axis_id()
integer, dimension(4) :: get_axis_id

if(.not.module_is_initialized) then
  call error_mesg('get_axis_id','spectral_diagnostics_init has not been called.', FATAL)
endif
get_axis_id = axis_id
return
end function get_axis_id
!===================================================================================
subroutine spectral_diagnostics_end

if(.not.module_is_initialized) return

deallocate(id_tr)

return
end subroutine spectral_diagnostics_end
!===================================================================================

end module spectral_dynamics_mod
