
module vert_diff_driver_mod

!-----------------------------------------------------------------------
!   module performs vertical diffusion of atmospheric variables
!-----------------------------------------------------------------------

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


use diag_manager_mod, only:  register_diag_field, send_data

use time_manager_mod, only:  time_type

use          mpp_mod, only:  input_nml_file
use          fms_mod, only:  file_exist, open_namelist_file, error_mesg,  &
                             check_nml_error, FATAL, mpp_pe, mpp_root_pe, &
                             close_file, write_version_number, stdlog

use    constants_mod, only:  CP_AIR, GRAV

use   field_manager_mod, only: MODEL_ATMOS
use  tracer_manager_mod, only: get_number_tracers, get_tracer_names, &
                               NO_TRACER

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

implicit none
private

public :: vert_diff_driver_down, vert_diff_driver_up,  &
          vert_diff_driver_init, vert_diff_driver_end
public :: surf_diff_type


!-----------------------------------------------------------------------
!---- namelist ----

logical :: do_conserve_energy         = .false.
logical :: do_mcm_no_neg_q            = .false.
logical :: use_virtual_temp_vert_diff = .true.
logical :: do_mcm_plev                = .false.
logical :: do_mcm_vert_diff_tq        = .false.

namelist /vert_diff_driver_nml/ do_conserve_energy,         &
                                do_mcm_no_neg_q,            &
                                use_virtual_temp_vert_diff, &
                                do_mcm_plev, do_mcm_vert_diff_tq

!-----------------------------------------------------------------------
! tracer storage is used 
type :: tracer_storage_type
   integer :: id_tr_dt     = 0 ! diag id of the tracer tendency due 
                               ! to vert diff
   integer :: id_tr_dt_int = 0 ! diag id of the vertically-integrated 
                               ! tracer tendency
   real, pointer :: &
        buffer(:,:,:) => NULL() ! buffer for tendency calculations
end type
type(tracer_storage_type), allocatable :: tr_store(:)

real, allocatable, dimension(:,:,:) :: dt_t_save, dt_q_save

!-------------------- diagnostics fields -------------------------------

integer :: id_tdt_vdif, id_qdt_vdif, id_udt_vdif, id_vdt_vdif,  &
           id_sens_vdif, id_evap_vdif,                          &
           id_tdt_diss_vdif, id_diss_heat_vdif

real :: missing_value = -999.

character(len=9), parameter :: mod_name = 'vert_diff'

!-----------------------------------------------------------------------
!---- version number ----

character(len=128) :: version = '$Id: vert_diff_driver.F90,v 20.0 2013/12/13 23:22:29 fms Exp $'
character(len=128) :: tagname = '$Name: tikal $'

logical :: module_is_initialized = .false.


contains

!#######################################################################

 subroutine vert_diff_driver_down (is, js, Time, delt, p_half, p_full, &
                                   z_full, diff_mom, diff_heat,        &
                                   u, v, t, q, trs,                    &
                                   dtau_du, dtau_dv, tau_x, tau_y,     &
                                   dt_u, dt_v, dt_t, dt_q, dt_trs,     &
                                   Surf_diff, diff_t_clubb, mask, kbot ) !cjg

integer, intent(in)                     :: is, js
type(time_type),   intent(in)           :: Time
real, intent(in)                        :: delt
real, intent(in)   , dimension(:,:,:)   :: p_half, p_full, z_full,  &
                                           diff_mom, diff_heat
real, intent(in),    dimension(:,:,:)   :: u, v, t, q
real, intent(in),    dimension(:,:,:,:) :: trs
real, intent(in),    dimension(:,:)     :: dtau_du, dtau_dv

real, intent(inout), dimension(:,:)     :: tau_x, tau_y
real, intent(inout), dimension(:,:,:)   :: dt_u, dt_v, dt_t, dt_q
real, intent(inout), dimension(:,:,:,:) :: dt_trs

type(surf_diff_type), intent(inout)     :: Surf_diff

!-->cjg
real, intent(in)   , dimension(:,:,:), optional :: diff_t_clubb
!<--cjg
real   , intent(in), dimension(:,:,:), optional :: mask
integer, intent(in), dimension(:,:),   optional :: kbot

real, dimension(size(t,1),size(t,2),size(t,3)) :: tt, dpg, q_2
real, dimension(size(t,1),size(t,2),size(t,3)) :: dissipative_heat
integer :: k, ntp, tr
logical :: used
real, dimension(size(t,1),size(t,2)) :: diag2
integer :: ie, je

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

  if (.not. module_is_initialized) call error_mesg       &
                  ('vert_diff_driver_mod',  &
                   'vert_diff_driver_init must be called first', FATAL)

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

  ntp = size(dt_trs,4) ! number of prognostic tracers
  if (size(trs,4) < ntp) call error_mesg ('vert_diff_driver', &
             'Number of tracers .lt. number of tracer tendencies',FATAL)

  ie = is + size(t,1) -1
  je = js + size(t,2) -1


    if(do_mcm_vert_diff_tq) then
      dt_t_save(is:ie,js:je,:) = dt_t
      dt_q_save(is:ie,js:je,:) = dt_q
      dt_t = 0.0
      dt_q = 0.0
    endif

!-----------------------------------------------------------------------
!---- to do diagnostics on dt_t, dt_q, dt_u, and dt_v at this point add 
!-----in the negative value of the field.  Note that the multiplication
!---- by 2 is necessary to get the correct value because sum_diag_phys
!---- is called twice for this single field
!-----------------------------------------------------------------------




!------- diagnostics for uwnd_diff -------
    if ( id_udt_vdif > 0 ) then
       used = send_data ( id_udt_vdif, -2.*dt_u, Time, is, js, 1, &
                          rmask=mask )
    endif

!------- diagnostics for vwnd_diff -------
    if ( id_vdt_vdif > 0 ) then
       used = send_data ( id_vdt_vdif, -2.*dt_v, Time, is, js, 1, &
                          rmask=mask )
    endif


!------- diagnostics for dt/dt_diff -------
    if ( id_tdt_vdif > 0 ) then
       used = send_data ( id_tdt_vdif, -3.*dt_t, Time, is, js, 1, &
                           rmask=mask )
    endif
    
!------- diagnostics for dq/dt_diff -------
    if ( id_qdt_vdif > 0 ) then
       used = send_data ( id_qdt_vdif, -2.*dt_q, Time, is, js, 1, &
                          rmask=mask )
    endif

    ! store values of tracer tendencies before the vertical diffusion, if 
    ! requested -- availability of storage serves as an indicatior that 
    ! storing is necessary
    do tr = 1,size(tr_store(:))
       if( associated(tr_store(tr)%buffer) ) &
            tr_store(tr)%buffer(is:ie,js:je,:) = dt_trs(:,:,:,tr)
    enddo

!-----------------------------------------------------------------------
!---- local temperature ----
!     (avoids divid by zero when computing virtual temp)

   tt = t
   if (present(mask)) where (mask < 0.5) tt = 200.

!-----------------------------------------------------------------------
!---- momentum diffusion ----
!---- heat/moisture diffusion (down only) ----
!---- tracer diffusion (no surface flux) ----

 q_2 = q
 if (do_mcm_no_neg_q) then
   where (q_2 < 0.0)  q_2 = 0.0
 endif

!--> cjg
 call gcm_vert_diff_down (is, js, delt, u, v, tt, q_2, trs(:,:,:,1:ntp), &
                          diff_mom, diff_heat,                           &
                          p_half, p_full, z_full,                        &
                          tau_x, tau_y, dtau_du, dtau_dv,                &
                          dt_u, dt_v, dt_t, dt_q, dt_trs(:,:,:,1:ntp),   &
                          dissipative_heat, Surf_diff,  diff_t_clubb, kbot ) !cjg
!<--cjg

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!---- to do diagnostics on dt_u, and dt_v at this point add 
!-----in the value of the field.  Note that the multiplication
!---- by 2 is necessary to get the correct value because sum_diag_phys
!---- is called twice for this single field
!-----------------------------------------------------------------------

!------ preliminary calculations for vert integrals -----
    if ( id_sens_vdif > 0 .or. id_evap_vdif > 0 .or. id_diss_heat_vdif > 0 ) then
            do k = 1, size(p_half,3)-1
               dpg(:,:,k) = (p_half(:,:,k+1)-p_half(:,:,k))/GRAV
            enddo
            if (present(mask)) dpg = dpg*mask
    endif
    

!------- diagnostics for sens_diff -------
    if ( id_sens_vdif > 0 ) then
!         --- compute column changes ---
          diag2 = sum( dt_t*dpg, 3 )
          used = send_data ( id_sens_vdif, -2.*CP_AIR*diag2, Time, is, js )
    endif

!------- diagnostics for evap_diff -------
    if ( id_evap_vdif > 0 ) then
!         --- compute column changes ---
          diag2 = sum( dt_q*dpg, 3 )
          used = send_data ( id_evap_vdif, -2.*diag2, Time, is, js )
    endif
    
    

!------- diagnostics for uwnd_diff -------
    if ( id_udt_vdif > 0 ) then
       used = send_data ( id_udt_vdif, 2.*dt_u, Time, is, js, 1, &
                          rmask=mask )
    endif

!------- diagnostics for vwnd_diff -------
    if ( id_vdt_vdif > 0 ) then
       used = send_data ( id_vdt_vdif, 2.*dt_v, Time, is, js, 1, &
                          rmask=mask )
    endif
    
!------- diagnostics for dissipative heating -------
    if ( id_tdt_diss_vdif > 0 ) then
       used = send_data ( id_tdt_diss_vdif, dissipative_heat, Time, &
                          is, js, 1, &
                          rmask=mask)
    endif

!------- diagnostics for dt/dt_diff -------
    if ( id_tdt_vdif > 0 ) then
       used = send_data ( id_tdt_vdif, -3.*dissipative_heat, Time, is, js, 1, &
                           rmask=mask )
    endif

!------- diagnostics for vertically integrated dissipative heating -------
    if ( id_diss_heat_vdif > 0 ) then
          diag2 = sum( CP_AIR*dissipative_heat*dpg, 3 )
          used = send_data ( id_diss_heat_vdif, diag2, Time, is, js )
    endif

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

 end subroutine vert_diff_driver_down

!#######################################################################

 subroutine vert_diff_driver_up (is, js, Time, delt, p_half, &
                                 Surf_diff, dt_t, dt_q, dt_tr, mask, kbot)

 integer,           intent(in)            :: is, js
 type(time_type),   intent(in)            :: Time
 real,    intent(in)                      :: delt
 real,    intent(in),    dimension(:,:,:) :: p_half
 type(surf_diff_type),   intent(in)       :: Surf_diff
 real,    intent(inout), dimension(:,:,:) :: dt_t, dt_q
 real,    intent(inout), dimension(:,:,:,:) :: dt_tr
 real   , intent(in), dimension(:,:,:), optional :: mask
 integer, intent(in),    dimension(:,:), optional :: kbot

 integer :: k, tr
 logical :: used
 real, dimension(size(p_half,1),size(p_half,2),size(p_half,3)-1) :: dpg
 real, dimension(size(p_half,1),size(p_half,2)) :: diag2
 integer :: ie, je

!-----------------------------------------------------------------------
    ie = is + size(p_half,1) -1
    je = js + size(p_half,2) -1
!-----------------------------------------------------------------------

    call gcm_vert_diff_up (is, js, delt, Surf_diff, dt_t, dt_q, dt_tr, kbot)


!-----------------------------------------------------------------------
!---- to do diagnostics on dt_t and dt_q at this point add in the
!---- the postive value of the field.  Note that the multiplication
!---- by 2 is necessary to get the correct value because sum_diag_phys
!---- is called twice for this single field
!-----------------------------------------------------------------------
!------- diagnostics for dt/dt_diff -------

    if ( id_tdt_vdif > 0 ) then
       used = send_data ( id_tdt_vdif, 3.*dt_t, Time, is, js, 1, &
                          rmask=mask )
    endif

!------- diagnostics for dq/dt_diff -------
    if ( id_qdt_vdif > 0 ) then
       used = send_data ( id_qdt_vdif, 2.*dt_q, Time, is, js, 1, &
                          rmask=mask )
    endif

!------ preliminary calculations for vert integrals -----
    if ( id_sens_vdif > 0 .or. id_evap_vdif > 0 ) then
            do k = 1, size(p_half,3)-1
               dpg(:,:,k) = (p_half(:,:,k+1)-p_half(:,:,k))/GRAV
            enddo
            if (present(mask)) dpg = dpg*mask
    endif

!------- diagnostics for sens_diff -------
    if ( id_sens_vdif > 0 ) then
!         --- compute column changes ---
          diag2 = sum( dt_t*dpg, 3 )
          used = send_data ( id_sens_vdif, 2.*CP_AIR*diag2, Time, is, js )
    endif

!------- diagnostics for evap_diff -------
    if ( id_evap_vdif > 0 ) then
!         --- compute column changes ---
          diag2 = sum( dt_q*dpg, 3 )
          used = send_data ( id_evap_vdif, 2.*diag2, Time, is, js )
    endif

!------- diagnostics of tracer tendencies ------- 
    do tr = 1, size(tr_store(:))
       if(tr_store(tr)%id_tr_dt > 0) then
          used = send_data(tr_store(tr)%id_tr_dt, &
               dt_tr(:,:,:,tr)-tr_store(tr)%buffer(is:ie,js:je,:),Time,is,js)
       endif
       
       if(tr_store(tr)%id_tr_dt_int > 0) then
          diag2 = sum((dt_tr(:,:,:,tr)-tr_store(tr)%buffer(is:ie,js:je,:))*dpg,3)
          used = send_data(tr_store(tr)%id_tr_dt_int, diag2, Time, is, js)
       endif
    enddo

    if(do_mcm_vert_diff_tq) then
      dt_t = dt_t + dt_t_save(is:ie,js:je,:)
      dt_q = dt_q + dt_q_save(is:ie,js:je,:)
    endif

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

 end subroutine vert_diff_driver_up

!#######################################################################

 subroutine vert_diff_driver_init ( Surf_diff, idim, jdim, kdim,  &
                                    axes, Time, do_clubb )  !cjg

 type(surf_diff_type), intent(inout) :: Surf_diff
 integer             , intent(in)    :: idim, jdim, kdim, axes(4)
 type(time_type)     , intent(in)    :: Time
!-->cjg
 integer, intent(in)                 :: do_clubb
!<--cjg

 integer :: unit, io, ierr, tr, logunit
 integer :: ntprog ! number of prognostic tracers in the atmosphere
 character(len=32)  :: name, units ! name of the tracer
 character(len=128) :: longname    ! long name of the tracer

!-----------------------------------------------------------------------
!------ read namelist ------

#ifdef INTERNAL_FILE_NML
   read (input_nml_file, nml=vert_diff_driver_nml, iostat=io)
   ierr = check_nml_error(io,'vert_diff_driver_nml')
#else   
   if ( file_exist('input.nml')) then
      unit = open_namelist_file ()
      ierr=1; do while (ierr /= 0)
         read  (unit, nml=vert_diff_driver_nml, iostat=io, end=10)
         ierr = check_nml_error(io,'vert_diff_driver_nml')
      enddo
 10   call close_file (unit)
   endif
#endif

!--------- write version number and namelist ------------------

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

!-------- initialize gcm vertical diffusion ------

   call vert_diff_init (Surf_diff, idim, jdim, kdim, do_conserve_energy, &
                        use_virtual_temp_vert_diff, do_mcm_plev, do_clubb) !cjg

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

   if(do_mcm_vert_diff_tq) then
     allocate(dt_t_save(idim,jdim,kdim)) ; dt_t_save = 0.0
     allocate(dt_q_save(idim,jdim,kdim)) ; dt_q_save = 0.0
   endif

!--------------- initialize diagnostic fields --------------------

   id_tdt_vdif = &
   register_diag_field ( mod_name, 'tdt_vdif', axes(1:3), Time, &
                        'Temperature tendency from vert diff',  &
                        'deg_K/s', missing_value=missing_value  )

   id_qdt_vdif = &
   register_diag_field ( mod_name, 'qdt_vdif', axes(1:3), Time, &
                        'Spec humidity tendency from vert diff',&
                        'kg/kg/s', missing_value=missing_value  )

   id_udt_vdif = &
   register_diag_field ( mod_name, 'udt_vdif', axes(1:3), Time, &
                        'Zonal wind tendency from vert diff',   &
                        'm/s2', missing_value=missing_value     )

   id_vdt_vdif = &
   register_diag_field ( mod_name, 'vdt_vdif', axes(1:3), Time,    &
                        'Meridional wind tendency from vert diff', &
                        'm/s2', missing_value=missing_value        )

   id_sens_vdif = &
   register_diag_field ( mod_name, 'sens_vdif', axes(1:2), Time,  &
                        'Integrated heat flux from vert diff',    &
                        'W/m2' )

   id_evap_vdif = &
   register_diag_field ( mod_name, 'evap_vdif', axes(1:2), Time,    &
                        'Integrated moisture flux from vert diff',  &
                        'kg/m2/s' )

   id_tdt_diss_vdif = &
   register_diag_field ( mod_name, 'tdt_diss_vdif', axes(1:3), Time,  &
                        'Dissipative heating from vert_diff', 'deg_K/s', &
                         missing_value=missing_value  ) 

   id_diss_heat_vdif = &
   register_diag_field ( mod_name, 'diss_heat_vdif', axes(1:2), Time,  &
                        'Integrated dissipative heating from vert diff',  &
                        'W/m2' )

   ! initialize diagnostics tracers
   call get_number_tracers(MODEL_ATMOS, num_prog=ntprog)
   allocate(tr_store(ntprog))
   do tr = 1,ntprog
      call get_tracer_names( MODEL_ATMOS, tr, name, longname, units )
      tr_store(tr)%id_tr_dt = &
        register_diag_field ( mod_name, trim(name)//'dt_vdif', axes(1:3), Time, &
           'Tendency of '//trim(longname)//' from vert diff', trim(units)//'/s')
      tr_store(tr)%id_tr_dt_int = &
        register_diag_field ( mod_name, trim(name)//'dt_vint_vdif', axes(1:2), Time, &
           'Integrated tendency of '//trim(longname)//' from vert diff',&
           trim(units)//' kg/(m2 s)')

      if(tr_store(tr)%id_tr_dt>0 .or.tr_store(tr)%id_tr_dt_int>0 ) &
           allocate(tr_store(tr)%buffer(idim,jdim,kdim))
   enddo

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

   module_is_initialized = .true.

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

 end subroutine vert_diff_driver_init

!#######################################################################

 subroutine vert_diff_driver_end

   integer :: tr ! tracer index

   call vert_diff_end
   if(do_mcm_vert_diff_tq) deallocate(dt_t_save, dt_q_save)
   ! deallocate tracer diagnostics storage
   do tr = 1,size(tr_store(:))
      if(associated(tr_store(tr)%buffer)) &
           deallocate(tr_store(tr)%buffer)
   enddo
   deallocate(tr_store)

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

   module_is_initialized = .false.

!-----------------------------------------------------------------------
 end subroutine vert_diff_driver_end

!#######################################################################

end module vert_diff_driver_mod

