module idealized_moist_phys_mod

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

use fms_mod, only: write_version_number, file_exist, close_file, stdlog, error_mesg, FATAL

use           constants_mod, only: grav

use        time_manager_mod, only: time_type, get_time, operator( + )

use    vert_turb_driver_mod, only: vert_turb_driver_init, vert_turb_driver, vert_turb_driver_end

use           vert_diff_mod, only: vert_diff_init, gcm_vert_diff_down, gcm_vert_diff_up, vert_diff_end, surf_diff_type

use two_stream_gray_rad_mod, only: two_stream_gray_rad_init, two_stream_gray_rad_down, two_stream_gray_rad_up, two_stream_gray_rad_end

use         mixed_layer_mod, only: mixed_layer_init, mixed_layer, mixed_layer_end

use         lscale_cond_mod, only: lscale_cond_init, lscale_cond, lscale_cond_end

use qe_moist_convection_mod, only: qe_moist_convection_init, qe_moist_convection, qe_moist_convection_end

use        diag_manager_mod, only: register_diag_field, send_data

use          transforms_mod, only: get_grid_domain

use   spectral_dynamics_mod, only: get_axis_id, get_num_levels, get_surf_geopotential

use        surface_flux_mod, only: surface_flux

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

character(len=128) :: version= &
'$Id: idealized_moist_phys.F90,v 20.0 2013/12/13 23:25:39 fms Exp $'

character(len=128) :: tagname= &
'$Name: tikal $'
character(len=10), parameter :: mod_name='atmosphere'

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

public :: idealized_moist_phys_init , idealized_moist_phys , idealized_moist_phys_end

logical :: module_is_initialized =.false.
logical :: turb = .false.
logical :: do_virtual = .false. ! whether virtual temp used in gcm_vert_diff
logical :: lwet_convection = .false.
logical :: two_stream = .true.
logical :: mixed_layer_bc = .false.
real :: roughness_heat = 0.05
real :: roughness_moist = 0.05
real :: roughness_mom = 0.05

namelist / idealized_moist_phys_nml / turb, lwet_convection, roughness_heat, two_stream, mixed_layer_bc, &
                                      roughness_moist, roughness_mom, do_virtual

real, allocatable, dimension(:,:)   ::                                        &
     z_surf,               &   ! surface height
     t_surf,               &   ! surface temperature
     q_surf,               &   ! surface moisture
     u_surf,               &   ! surface U wind
     v_surf,               &   ! surface V wind
     rough_mom,            &   ! momentum roughness length for surface_flux
     rough_heat,           &   ! heat roughness length for surface_flux
     rough_moist,          &   ! moisture roughness length for surface_flux
     gust,                 &   ! gustiness constant
     z_pbl,                &   ! gustiness constant
     flux_t,               &   ! surface sensible heat flux
     flux_q,               &   ! surface moisture flux
     flux_r,               &   ! surface radiation flux
     flux_u,               &   ! surface flux of zonal mom.
     flux_v,               &   ! surface flux of meridional mom.
     drag_m,               &   ! momentum drag coefficient
     drag_t,               &   ! heat drag coefficient
     drag_q,               &   ! moisture drag coefficient
     w_atm,                &   ! wind speed
     ustar,                &   ! friction velocity
     bstar,                &   ! buoyancy scale
     qstar,                &   ! moisture scale
     dhdt_surf,            &   ! d(sensible heat flux)/d(surface temp)
     dedt_surf,            &   ! d(latent heat flux)/d(surface temp)???
     dedq_surf,            &   ! d(latent heat flux)/d(surface moisture)???
     drdt_surf,            &   ! d(upward longwave)/d(surface temp)
     dhdt_atm,             &   ! d(sensible heat flux)/d(atmos.temp)
     dedq_atm,             &   ! d(latent heat flux)/d(atmospheric mixing rat.)
     dtaudv_atm,           &   ! d(stress component)/d(atmos wind)
     dtaudu_atm,           &   ! d(stress component)/d(atmos wind)
     fracland,             &   ! fraction of land in gridbox
     rough                     ! roughness for vert_turb_driver


real, allocatable, dimension(:,:,:) ::                                        &
     diff_m,               &   ! momentum diffusion coeff.
     diff_t,               &   ! temperature diffusion coeff.
     tdtlw,                &   ! place holder. appears in calling arguments of vert_turb_driver but not used unless do_edt=.true. -- pjp
     diss_heat,            &   ! heat dissipated by vertical diffusion
     non_diff_dt_ug,       &   ! zonal wind tendency except from vertical diffusion
     non_diff_dt_vg,       &   ! merid. wind tendency except from vertical diffusion
     non_diff_dt_tg,       &   ! temperature tendency except from vertical diffusion
     non_diff_dt_qg,       &   ! moisture tendency except from vertical diffusion
     conv_dt_tg,           &   ! temperature tendency from convection
     conv_dt_qg,           &   ! moisture tendency from convection
     cond_dt_tg,           &   ! temperature tendency from condensation
     cond_dt_qg                ! moisture tendency from condensation


logical, allocatable, dimension(:,:) ::                                       &
     avail,                &   ! generate surf. flux (all true)
     land,                 &   ! land points (all false)
     coldT,                &   ! should precipitation be snow at this point
     convect                   ! place holder. appears in calling arguments of vert_turb_driver but not used unless do_entrain=.true. -- pjp

real, allocatable, dimension(:,:) ::                                          &
     klzbs,                &   ! stored level of zero buoyancy values
     cape,                 &   ! convectively available potential energy
     cin,                  &   ! convective inhibition (this and the above are before the adjustment)
     invtau_q_relaxation,  &   ! temperature relaxation time scale
     invtau_t_relaxation,  &   ! humidity relaxation time scale
     rain,                 &   !
     snow

real, allocatable, dimension(:,:,:) :: &
     t_ref,          &   ! relaxation temperature for bettsmiller scheme
     q_ref               ! relaxation moisture for bettsmiller scheme

real, allocatable, dimension(:,:) :: &
     net_surf_sw_down,  &   ! net sw flux at surface
     surf_lw_down           ! downward lw flux at surface

integer ::           &
     id_diff_dt_ug,  &   ! zonal wind tendency from vertical diffusion
     id_diff_dt_vg,  &   ! merid. wind tendency from vertical diffusion
     id_diff_dt_tg,  &   ! temperature tendency from vertical diffusion
     id_diff_dt_qg,  &   ! moisture tendency from vertical diffusion
     id_conv_rain,   &   ! rain from convection
     id_cond_rain,   &   ! rain from condensation
     id_conv_dt_tg,  &   ! temperature tendency from convection
     id_conv_dt_qg,  &   ! temperature tendency from convection
     id_cond_dt_tg,  &   ! temperature tendency from convection
     id_cond_dt_qg       ! temperature tendency from convection

integer, allocatable, dimension(:,:) :: convflag ! indicates which qe convection subroutines are used
real,    allocatable, dimension(:,:) :: rad_lat   

type(surf_diff_type) :: Tri_surf ! used by gcm_vert_diff

logical :: used, doing_edt, doing_entrain
integer, dimension(4) :: axes
integer :: is, ie, js, je, num_levels, nsphum, dt_integer
real :: dt_real
type(time_type) :: Time_step

!=================================================================================================================================
contains
!=================================================================================================================================

subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lat_2d, rad_lonb_2d, rad_latb_2d, t_surf_init)
type(time_type), intent(in) :: Time, Time_step_in
integer, intent(in) :: nhum
real, intent(in), dimension(:,:) :: rad_lat_2d, rad_lonb_2d, rad_latb_2d, t_surf_init

integer :: io, nml_unit, stdlog_unit, seconds, days

if(module_is_initialized) return

call write_version_number(version, tagname)

#ifdef INTERNAL_FILE_NML
   read (input_nml_file, nml=idealized_moist_phys_nml, iostat=io)
#else  
   if ( file_exist('input.nml') ) then
      nml_unit = open_namelist_file()
      read (nml_unit, idealized_moist_phys_nml, iostat=io)
      call close_file(nml_unit)
   endif
#endif
stdlog_unit = stdlog()
write(stdlog_unit, idealized_moist_phys_nml)

nsphum = nhum
Time_step = Time_step_in
call get_time(Time_step, seconds, days)
dt_integer   = 86400*days + seconds
dt_real      = float(dt_integer)

call get_grid_domain(is, ie, js, je)
call get_num_levels(num_levels)

allocate(rad_lat     (is:ie, js:je)); rad_lat = rad_lat_2d
allocate(z_surf      (is:ie, js:je))
allocate(t_surf      (is:ie, js:je))
allocate(q_surf      (is:ie, js:je)); q_surf = 0.0
allocate(u_surf      (is:ie, js:je)); u_surf = 0.0
allocate(v_surf      (is:ie, js:je)); v_surf = 0.0
allocate(rough_mom   (is:ie, js:je)); rough_mom = roughness_mom
allocate(rough_heat  (is:ie, js:je)); rough_heat = roughness_heat
allocate(rough_moist (is:ie, js:je)); rough_moist = roughness_moist
allocate(gust        (is:ie, js:je)); gust = 1.0
allocate(z_pbl       (is:ie, js:je))
allocate(flux_t      (is:ie, js:je))
allocate(flux_q      (is:ie, js:je))
allocate(flux_r      (is:ie, js:je))
allocate(flux_u      (is:ie, js:je))
allocate(flux_v      (is:ie, js:je))
allocate(drag_m      (is:ie, js:je))
allocate(drag_t      (is:ie, js:je))
allocate(drag_q      (is:ie, js:je))
allocate(w_atm       (is:ie, js:je))
allocate(ustar       (is:ie, js:je))
allocate(bstar       (is:ie, js:je))
allocate(qstar       (is:ie, js:je))
allocate(dhdt_surf   (is:ie, js:je))
allocate(dedt_surf   (is:ie, js:je))
allocate(dedq_surf   (is:ie, js:je))
allocate(drdt_surf   (is:ie, js:je))
allocate(dhdt_atm    (is:ie, js:je))
allocate(dedq_atm    (is:ie, js:je))
allocate(dtaudv_atm  (is:ie, js:je))
allocate(dtaudu_atm  (is:ie, js:je))
allocate(land        (is:ie, js:je)); land = .false.
allocate(avail       (is:ie, js:je)); avail = .true.
allocate(fracland    (is:ie, js:je)); fracland = 0.0
allocate(rough       (is:ie, js:je))
allocate(diff_t      (is:ie, js:je, num_levels))
allocate(diff_m      (is:ie, js:je, num_levels))
allocate(diss_heat   (is:ie, js:je, num_levels))
allocate(tdtlw       (is:ie, js:je, num_levels)); tdtlw = 0.0

allocate(non_diff_dt_ug  (is:ie, js:je, num_levels))
allocate(non_diff_dt_vg  (is:ie, js:je, num_levels))
allocate(non_diff_dt_tg  (is:ie, js:je, num_levels))
allocate(non_diff_dt_qg  (is:ie, js:je, num_levels))

allocate(net_surf_sw_down        (is:ie, js:je))
allocate(surf_lw_down            (is:ie, js:je))
allocate(conv_dt_tg  (is:ie, js:je, num_levels))
allocate(conv_dt_qg  (is:ie, js:je, num_levels))
allocate(cond_dt_tg  (is:ie, js:je, num_levels))
allocate(cond_dt_qg  (is:ie, js:je, num_levels))

allocate(coldT        (is:ie, js:je)); coldT = .false.
allocate(klzbs        (is:ie, js:je))
allocate(cape         (is:ie, js:je))
allocate(cin          (is:ie, js:je))
allocate(invtau_q_relaxation  (is:ie, js:je))
allocate(invtau_t_relaxation  (is:ie, js:je))
allocate(rain         (is:ie, js:je)); rain = 0.0
allocate(snow         (is:ie, js:je)); snow = 0.0
allocate(convflag     (is:ie, js:je))
allocate(convect      (is:ie, js:je)); convect = .false.

allocate(t_ref (is:ie, js:je, num_levels)); t_ref = 0.0
allocate(q_ref (is:ie, js:je, num_levels)); q_ref = 0.0

call get_surf_geopotential(z_surf)
z_surf = z_surf/grav

if(mixed_layer_bc) then
  ! need an initial condition for the mixed layer temperature
  ! may be overwritten by restart file
  ! choose an unstable initial condition to allow moisture
  ! to quickly enter the atmosphere avoiding problems with the convection scheme
  t_surf = t_surf_init + 1.0
  call mixed_layer_init(is, ie, js, je, num_levels, t_surf, get_axis_id(), Time) ! t_surf is intent(inout)
endif

if(turb) then
! need to call vert_diff_init even if using gcm_vert_diff (rather than
! gcm_vert_diff_down) because the variable sphum is not initialized
! otherwise in the vert_diff module
   call vert_diff_init (Tri_surf, ie-is+1, je-js+1, num_levels, .true., do_virtual)
end if

call lscale_cond_init()

axes = get_axis_id()

id_cond_dt_qg = register_diag_field(mod_name, 'dt_qg_condensation',        &
     axes(1:3), Time, 'Moisture tendency from condensation','kg/kg/s')
id_cond_dt_tg = register_diag_field(mod_name, 'dt_tg_condensation',        &
     axes(1:3), Time, 'Temperature tendency from condensation','K/s')
id_cond_rain = register_diag_field(mod_name, 'condensation_rain',          &
     axes(1:2), Time, 'Rain from condensation','kg/m/m/s')

if(lwet_convection) then
   call qe_moist_convection_init()
   id_conv_dt_qg = register_diag_field(mod_name, 'dt_qg_convection',          &
        axes(1:3), Time, 'Moisture tendency from convection','kg/kg/s')
   id_conv_dt_tg = register_diag_field(mod_name, 'dt_tg_convection',          &
        axes(1:3), Time, 'Temperature tendency from convection','K/s')
   id_conv_rain = register_diag_field(mod_name, 'convection_rain',            &
        axes(1:2), Time, 'Rain from convection','kg/m/m/s')
endif

if(two_stream) call two_stream_gray_rad_init(is, ie, js, je, num_levels, get_axis_id(), Time)

if(turb) then
   call vert_turb_driver_init (rad_lonb_2d, rad_latb_2d, ie-is+1,je-js+1, &
                 num_levels,get_axis_id(),Time, doing_edt, doing_entrain)

   axes = get_axis_id()
   id_diff_dt_ug = register_diag_field(mod_name, 'dt_ug_diffusion',        &
        axes(1:3), Time, 'zonal wind tendency from diffusion','m/s^2')
   id_diff_dt_vg = register_diag_field(mod_name, 'dt_vg_diffusion',        &
        axes(1:3), Time, 'meridional wind tendency from diffusion','m/s^2')
   id_diff_dt_tg = register_diag_field(mod_name, 'dt_tg_diffusion',        &
        axes(1:3), Time, 'temperature diffusion tendency','T/s')
   id_diff_dt_qg = register_diag_field(mod_name, 'dt_qg_diffusion',        &
        axes(1:3), Time, 'moisture diffusion tendency','T/s')
endif

end subroutine idealized_moist_phys_init
!=================================================================================================================================
subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg, grid_tracers, &
                                previous, current, dt_ug, dt_vg, dt_tg, dt_tracers)

type(time_type),            intent(in)    :: Time
real, dimension(:,:,:,:),   intent(in)    :: p_half, p_full, z_half, z_full, ug, vg, tg
real, dimension(:,:,:,:,:), intent(in)    :: grid_tracers
integer,                    intent(in)    :: previous, current
real, dimension(:,:,:),     intent(inout) :: dt_ug, dt_vg, dt_tg
real, dimension(:,:,:,:),   intent(inout) :: dt_tracers

real :: delta_t
real, dimension(size(ug,1), size(ug,2), size(ug,3)) :: tg_tmp, qg_tmp

if(current == previous) then
   delta_t = dt_real
else
   delta_t = 2*dt_real
endif

if (lwet_convection) then
   rain = 0.0; snow = 0.0
   call qe_moist_convection ( delta_t,              tg(:,:,:,previous),      &
    grid_tracers(:,:,:,previous,nsphum),        p_full(:,:,:,previous),      &
                          p_half(:,:,:,previous),                coldT,      &
                                 rain,                            snow,      &
                           conv_dt_tg,                      conv_dt_qg,      &
                                q_ref,                        convflag,      &
                                klzbs,                            cape,      &
                                  cin,             invtau_q_relaxation,      &
                  invtau_t_relaxation,                           t_ref)

   tg_tmp = conv_dt_tg + tg(:,:,:,previous)
   qg_tmp = conv_dt_qg + grid_tracers(:,:,:,previous,nsphum)
!  note the delta's are returned rather than the time derivatives

   conv_dt_tg = conv_dt_tg/delta_t
   conv_dt_qg = conv_dt_qg/delta_t
   rain       = rain/delta_t

   dt_tg = dt_tg + conv_dt_tg
   dt_tracers(:,:,:,nsphum) = dt_tracers(:,:,:,nsphum) + conv_dt_qg

   if(id_conv_dt_qg > 0) used = send_data(id_conv_dt_qg, conv_dt_qg, Time)
   if(id_conv_dt_tg > 0) used = send_data(id_conv_dt_tg, conv_dt_tg, Time)
   if(id_conv_rain  > 0) used = send_data(id_conv_rain, rain, Time)

else

   tg_tmp = tg(:,:,:,previous)
   qg_tmp = grid_tracers(:,:,:,previous,nsphum)

endif

rain = 0.0
call lscale_cond (         tg_tmp,                          qg_tmp,        &
           p_full(:,:,:,previous),          p_half(:,:,:,previous),        &
                            coldT,                            rain,        &
                             snow,                      cond_dt_tg,        &
                       cond_dt_qg )
                                                                          
cond_dt_tg = cond_dt_tg/delta_t
cond_dt_qg = cond_dt_qg/delta_t
rain       = rain/delta_t
                                                                             
dt_tg = dt_tg + cond_dt_tg
dt_tracers(:,:,:,nsphum) = dt_tracers(:,:,:,nsphum) + cond_dt_qg
                                                                               
if(id_cond_dt_qg > 0) used = send_data(id_cond_dt_qg, cond_dt_qg, Time)
if(id_cond_dt_tg > 0) used = send_data(id_cond_dt_tg, cond_dt_tg, Time)
if(id_cond_rain  > 0) used = send_data(id_cond_rain, rain, Time)

! Begin the radiation calculation by computing downward fluxes.
! This part of the calculation does not depend on the surface temperature.

if(two_stream) then
   call two_stream_gray_rad_down(is, js, Time, &
                       rad_lat(:,:),           &
                       p_half(:,:,:,current),  &
                       tg(:,:,:,previous),     &
                       net_surf_sw_down(:,:),  &
                       surf_lw_down(:,:))
end if

if(.not.mixed_layer_bc) then
                                                                               
!!$! infinite heat capacity
!    t_surf = surface_temperature_forced(rad_lat)
!!$! no heat capacity:
!!$   t_surf = tg(:,:,num_levels,previous)
                                                                              
!!$! surface temperature has same potential temp. as lowest layer:
!!$  t_surf = surface_temperature(tg(:,:,:,previous), p_full(:,:,:,current), p_half(:,:,:,current))
end if

call surface_flux(                                                          &
                  tg(:,:,num_levels,previous),                              &
 grid_tracers(:,:,num_levels,previous,nsphum),                              &
                  ug(:,:,num_levels,previous),                              &
                  vg(:,:,num_levels,previous),                              &
               p_full(:,:,num_levels,current),                              &
   z_full(:,:,num_levels,current)-z_surf(:,:),                              &
             p_half(:,:,num_levels+1,current),                              &
                                  t_surf(:,:),                              &
                                  t_surf(:,:),                              &
                                  q_surf(:,:),                              & ! is intent(inout)
                                  u_surf(:,:),                              &
                                  v_surf(:,:),                              &
                               rough_mom(:,:),                              &
                              rough_heat(:,:),                              &
                             rough_moist(:,:),                              &
                               rough_mom(:,:),                              & ! using rough_mom in place of rough_scale -- pjp
                                    gust(:,:),                              &
                                  flux_t(:,:),                              & ! is intent(out)
                                  flux_q(:,:),                              & ! is intent(out)
                                  flux_r(:,:),                              & ! is intent(out)
                                  flux_u(:,:),                              & ! is intent(out)
                                  flux_v(:,:),                              & ! is intent(out)
                                  drag_m(:,:),                              & ! is intent(out)
                                  drag_t(:,:),                              & ! is intent(out)
                                  drag_q(:,:),                              & ! is intent(out)
                                   w_atm(:,:),                              & ! is intent(out)
                                   ustar(:,:),                              & ! is intent(out)
                                   bstar(:,:),                              & ! is intent(out)
                                   qstar(:,:),                              & ! is intent(out)
                               dhdt_surf(:,:),                              & ! is intent(out)
                               dedt_surf(:,:),                              & ! is intent(out)
                               dedq_surf(:,:),                              & ! is intent(out)
                               drdt_surf(:,:),                              & ! is intent(out)
                                dhdt_atm(:,:),                              & ! is intent(out)
                                dedq_atm(:,:),                              & ! is intent(out)
                              dtaudu_atm(:,:),                              & ! is intent(out)
                              dtaudv_atm(:,:),                              & ! is intent(out)
                                      delta_t,                              &
                                    land(:,:),                              &
                               .not.land(:,:),                              &
                                   avail(:,:)  )

! Now complete the radiation calculation by computing the upward and net fluxes.

if(two_stream) then
   call two_stream_gray_rad_up(is, js, Time, &
                     rad_lat(:,:),           &
                     p_half(:,:,:,current),  &
                     t_surf(:,:),            &
                     tg(:,:,:,previous),     &
                     dt_tg(:,:,:))
end if

if(turb) then

   call vert_turb_driver(            1,                              1, &
                                  Time,                 Time+Time_step, &
                               delta_t, tdtlw(:,:,:),    fracland(:,:), &
                 p_half(:,:,:,current),          p_full(:,:,:,current), &
                 z_half(:,:,:,current),          z_full(:,:,:,current), &
                            ustar(:,:),                     bstar(:,:), &
                            qstar(:,:),                     rough(:,:), &
                          rad_lat(:,:),                   convect(:,:), &
                    ug(:,:,:,current ),             vg(:,:,:,current ), &
                    tg(:,:,:,current ),                                 &
    grid_tracers(:,:,:,current,nsphum),  grid_tracers(:,:,:,current,:), &
                    ug(:,:,:,previous),                                 &
                    vg(:,:,:,previous),             tg(:,:,:,previous), &
   grid_tracers(:,:,:,previous,nsphum), grid_tracers(:,:,:,previous,:), &
                          dt_ug(:,:,:),                   dt_vg(:,:,:), &
                          dt_tg(:,:,:),       dt_tracers(:,:,:,nsphum), &
                   dt_tracers(:,:,:,:),                  diff_t(:,:,:), &
                         diff_m(:,:,:),                      gust(:,:), &
                            z_pbl(:,:) )
!
!! Don't zero these derivatives as the surface flux depends implicitly
!! on the lowest level values
!! However it should be noted that these derivatives do not take into
!! account the change in the Monin-Obukhov coefficients, and so are not
!! very accurate.
!
!!$   dtaudv_atm = 0.0
!!$   dhdt_atm   = 0.0
!!$   dedq_atm   = 0.0

   if(.not.mixed_layer_bc) then
     call error_mesg('atmosphere','no diffusion implentation for non-mixed layer b.c.',FATAL)
   endif


! We must use gcm_vert_diff_down and _up rather than gcm_vert_diff as the surface flux
! depends implicitly on the surface values

!
! Don't want to do time splitting for the implicit diffusion step in case
! of compensation of the tendencies
!
   non_diff_dt_ug  = dt_ug
   non_diff_dt_vg  = dt_vg
   non_diff_dt_tg  = dt_tg
   non_diff_dt_qg  = dt_tracers(:,:,:,nsphum)

   call gcm_vert_diff_down (1, 1,                                          &
                            delta_t,             ug(:,:,:,previous),       &
                            vg(:,:,:,previous),  tg(:,:,:,previous),       &
                            grid_tracers(:,:,:,previous,nsphum),           &
                            grid_tracers(:,:,:,previous,:), diff_m(:,:,:), &
                            diff_t(:,:,:),          p_half(:,:,:,current), &
                            p_full(:,:,:,current),  z_full(:,:,:,current), &
                            flux_u(:,:),                      flux_v(:,:), &
                            dtaudu_atm(:,:),              dtaudv_atm(:,:), &
                            dt_ug(:,:,:),                    dt_vg(:,:,:), &
                            dt_tg(:,:,:),        dt_tracers(:,:,:,nsphum), &
                            dt_tracers(:,:,:,:),         diss_heat(:,:,:), &
                            Tri_surf)
!
! update surface temperature
!
   call mixed_layer(                                                       &
                              Time,                                        &
                              t_surf(:,:),                                 & ! t_surf is intent(inout)
                              flux_t(:,:),                                 &
                              flux_q(:,:),                                 &
                              flux_r(:,:),                                 &
                                  dt_real,                                 &
                    net_surf_sw_down(:,:),                                 &
                        surf_lw_down(:,:),                                 &
                            Tri_surf,                                      & ! Tri_surf is intent(inout)
                           dhdt_surf(:,:),                                 &
                           dedt_surf(:,:),                                 &
                           dedq_surf(:,:),                                 &
                           drdt_surf(:,:),                                 &
                            dhdt_atm(:,:),                                 &
                            dedq_atm(:,:))


   call gcm_vert_diff_up (1, 1, delta_t, Tri_surf, dt_tg(:,:,:), dt_tracers(:,:,:,nsphum), dt_tracers(:,:,:,:))

   if(id_diff_dt_ug > 0) used = send_data(id_diff_dt_ug, dt_ug - non_diff_dt_ug, Time)
   if(id_diff_dt_vg > 0) used = send_data(id_diff_dt_vg, dt_vg - non_diff_dt_vg, Time)
   if(id_diff_dt_tg > 0) used = send_data(id_diff_dt_tg, dt_tg - non_diff_dt_tg, Time)
   if(id_diff_dt_qg > 0) used = send_data(id_diff_dt_qg, dt_tracers(:,:,:,nsphum) - non_diff_dt_qg, Time)

endif ! if(turb) then

end subroutine idealized_moist_phys
!=================================================================================================================================
subroutine idealized_moist_phys_end

if(two_stream)      call two_stream_gray_rad_end
if(lwet_convection) call qe_moist_convection_end
if(turb) then
   call vert_diff_end
   call vert_turb_driver_end
endif
call lscale_cond_end
if(mixed_layer_bc)  call mixed_layer_end(t_surf)

end subroutine idealized_moist_phys_end
!=================================================================================================================================

end module idealized_moist_phys_mod
