pressuresource.F90

Go to the documentation of this file.

module pressuresource_mod
  use monc_component_mod, only : component_descriptor_type
  use grids_mod, only : z_index, y_index, x_index
  use state_mod, only : model_state_type
  use datadefn_mod, only : default_precision, precision_type
  use registry_mod, only : is_component_enabled
  use logging_mod, only : log_error, log_master_log
  use mpi, only : mpi_request_null
  implicit none

#ifndef TEST_MODE
  private
#endif

  real(kind=DEFAULT_PRECISION), dimension(:,:), allocatable :: send_buffer_x, send_buffer_y

  public pressuresource_get_descriptor
contains

  type(component_descriptor_type) function pressuresource_get_descriptor()
    pressuresource_get_descriptor%name="psrce"
    pressuresource_get_descriptor%version=0.1
    pressuresource_get_descriptor%initialisation=>initialisation_callback
    pressuresource_get_descriptor%timestep=>timestep_callback
    pressuresource_get_descriptor%finalisation=>finalisation_callback
  end function pressuresource_get_descriptor

  subroutine initialisation_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state

    if (.not. is_component_enabled(current_state%options_database, "diverr")) then
      call log_master_log(log_error, "The pressure source component requires the diverr component to be enabled")
    end if  

#ifdef U_ACTIVE
    allocate(send_buffer_x(current_state%local_grid%size(z_index)-1, current_state%local_grid%size(y_index)), &
         current_state%psrce_recv_buffer_x(current_state%local_grid%size(z_index)-1, current_state%local_grid%size(y_index)))
#endif
#ifdef V_ACTIVE
    allocate(send_buffer_y(current_state%local_grid%size(z_index)-1, current_state%local_grid%size(x_index)), &
         current_state%psrce_recv_buffer_y(current_state%local_grid%size(z_index)-1, current_state%local_grid%size(x_index)))
#endif
    current_state%psrce_x_hs_send_request=mpi_request_null
    current_state%psrce_y_hs_send_request=mpi_request_null
    current_state%psrce_x_hs_recv_request=mpi_request_null
    current_state%psrce_y_hs_recv_request=mpi_request_null
  end subroutine initialisation_callback

  subroutine timestep_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state

    if (current_state%first_timestep_column) call register_neighbouring_pressure_data_recv(current_state)
    if (.not. current_state%halo_column) call calculate_psrce(current_state)        
    if (current_state%last_timestep_column) call send_neighbouring_pressure_data(current_state)
  end subroutine timestep_callback

  subroutine finalisation_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state

    if (allocated(send_buffer_x)) deallocate(send_buffer_x)
    if (allocated(current_state%psrce_recv_buffer_x)) deallocate(current_state%psrce_recv_buffer_x)
    if (allocated(send_buffer_y)) deallocate(send_buffer_y)
    if (allocated(current_state%psrce_recv_buffer_y)) deallocate(current_state%psrce_recv_buffer_y)
  end subroutine finalisation_callback      

  subroutine calculate_psrce(current_state)
    type(model_state_type), target, intent(inout) :: current_state

    integer :: k, local_y, local_x, corrected_y, corrected_x
    logical :: last_x, last_y

    local_y=current_state%column_local_y
    local_x=current_state%column_local_x
    last_y = local_y == current_state%local_grid%local_domain_end_index(y_index)
    last_x = local_x == current_state%local_grid%local_domain_end_index(x_index)
    if (last_x .or. last_y) then
      corrected_x=local_x-current_state%local_grid%halo_size(x_index)
      corrected_y=local_y-current_state%local_grid%halo_size(y_index)
    end if
    do k=2,current_state%local_grid%size(z_index)
#ifdef W_ACTIVE
      current_state%p%data(k, local_y, local_x)=current_state%p%data(k, local_y, local_x)+&
           4.0_default_precision*(current_state%global_grid%configuration%vertical%tzc2(k)*&
           current_state%sw%data(k, local_y, local_x)-&
           current_state%global_grid%configuration%vertical%tzc1(k)*current_state%sw%data(k-1, local_y, local_x))      
#endif
#ifdef U_ACTIVE
      current_state%p%data(k, local_y, local_x)=current_state%p%data(k, local_y, local_x)+&
           current_state%global_grid%configuration%horizontal%cx * current_state%su%data(k, local_y, local_x)
#endif
#ifdef V_ACTIVE
      current_state%p%data(k, local_y, local_x)=current_state%p%data(k, local_y, local_x)+&
           current_state%global_grid%configuration%horizontal%cy * current_state%sv%data(k, local_y, local_x)
#endif
#ifdef U_ACTIVE
      if (local_x .gt. 3) then
        current_state%p%data(k, local_y, local_x)=current_state%p%data(k, local_y, local_x)-&
             current_state%global_grid%configuration%horizontal%cx * current_state%su%data(k, local_y, local_x-1)
      end if
      
      if (last_x) then
        send_buffer_x(k-1, corrected_y)=&
             current_state%global_grid%configuration%horizontal%cx * current_state%su%data(k, local_y, local_x)
      end if
#endif
#ifdef V_ACTIVE
      if (local_y .gt. 3 .and. local_x .gt. 3) then
        current_state%p%data(k, local_y, local_x)=current_state%p%data(k, local_y, local_x)-&
             current_state%global_grid%configuration%horizontal%cy * current_state%sv%data(k, local_y-1, local_x)      
      end if
      if (last_y) then
        send_buffer_y(k-1, corrected_x)=&
             current_state%global_grid%configuration%horizontal%cy * current_state%sv%data(k, local_y, local_x)
      end if
#endif
    end do
  end subroutine calculate_psrce

  subroutine send_neighbouring_pressure_data(current_state)
    type(model_state_type), target, intent(inout) :: current_state

    integer :: ierr

#ifdef U_ACTIVE
    if (current_state%local_grid%neighbours(x_index,3) .eq. current_state%parallel%my_rank) then
      current_state%psrce_recv_buffer_x=send_buffer_x
    else
      call mpi_isend(send_buffer_x, size(send_buffer_x), precision_type, current_state%local_grid%neighbours(x_index,3), &
           10, current_state%parallel%neighbour_comm, current_state%psrce_x_hs_send_request, ierr)
    end if
#endif
#ifdef V_ACTIVE
    if (current_state%local_grid%neighbours(y_index,3) .eq. current_state%parallel%my_rank) then
      current_state%psrce_recv_buffer_y=send_buffer_y
    else
      call mpi_isend(send_buffer_y, size(send_buffer_y), precision_type, current_state%local_grid%neighbours(y_index,3), &
           10, current_state%parallel%neighbour_comm, current_state%psrce_y_hs_send_request, ierr)
    end if
#endif
  end subroutine send_neighbouring_pressure_data

  subroutine register_neighbouring_pressure_data_recv(current_state)
    type(model_state_type), target, intent(inout) :: current_state

    integer :: ierr

#ifdef U_ACTIVE
    if (current_state%local_grid%neighbours(x_index,2) .ne. current_state%parallel%my_rank) then
      call mpi_irecv(current_state%psrce_recv_buffer_x, size(current_state%psrce_recv_buffer_x), precision_type, &
           current_state%local_grid%neighbours(x_index,2), 10, current_state%parallel%neighbour_comm, &
           current_state%psrce_x_hs_recv_request, ierr)
    end if    
#endif
#ifdef V_ACTIVE
    if (current_state%local_grid%neighbours(y_index,2) .ne. current_state%parallel%my_rank) then
      call mpi_irecv(current_state%psrce_recv_buffer_y, size(current_state%psrce_recv_buffer_y), precision_type, &
           current_state%local_grid%neighbours(y_index,2), 10, current_state%parallel%neighbour_comm, &
           current_state%psrce_y_hs_recv_request, ierr)
    end if    
#endif
  end subroutine register_neighbouring_pressure_data_recv
end module pressuresource_mod