pdf_analysis.F90

Go to the documentation of this file.

 
module pdf_analysis_mod
  use monc_component_mod, only : component_array_field_type, &
      component_double_data_type, component_descriptor_type, &
      component_field_value_type, component_field_information_type
  use state_mod, only : model_state_type
  use grids_mod, only : x_index, y_index, z_index
  use datadefn_mod, only : default_precision, precision_type
  use optionsdatabase_mod, only : options_has_key, options_get_logical, options_get_integer, options_get_string, options_get_real
  use mpi, only : mpi_sum, mpi_in_place, mpi_int, mpi_real, mpi_double, mpi_comm
  use logging_mod, only : log_info, log_debug, log_error, log_master_log, log_is_master
  use conversions_mod, only : conv_to_string
  implicit none
 
#ifndef TEST_MODE
  private
#endif
 
  integer :: start_x, end_x, start_y, end_y, xsize, ysize
 
  real(kind=default_precision) :: uppercrit, dwnpercrit
  real(kind=default_precision), dimension(:), allocatable :: tmp_all
  integer, dimension(:), allocatable :: gpts_on_proc, & ! number of horizontal grid points on each process
                                        displacements   ! displacement for mpi_gatherv
                                                        ! these are available on all processes
 
  integer :: tpts  ! total number of horizontal grid points on full domain
  integer :: lpts  ! local number of horizontal grid points on 
 
  logical :: show_critical_w  ! stdout diagnostic logical
 
  integer :: diagnostic_generation_frequency
 
  public pdf_analysis_get_descriptor
 
contains
 
 
  type(component_descriptor_type) function pdf_analysis_get_descriptor()
    pdf_analysis_get_descriptor%name="pdf_analysis"
    pdf_analysis_get_descriptor%version=0.1
    pdf_analysis_get_descriptor%initialisation=>init_callback
    pdf_analysis_get_descriptor%timestep=>timestep_callback
    pdf_analysis_get_descriptor%finalisation=>finalisation_callback
 
    pdf_analysis_get_descriptor%field_value_retrieval=>field_value_retrieval_callback
    pdf_analysis_get_descriptor%field_information_retrieval=>field_information_retrieval_callback
    allocate(pdf_analysis_get_descriptor%published_fields(2))
    pdf_analysis_get_descriptor%published_fields(1)="critical_updraft_local"
    pdf_analysis_get_descriptor%published_fields(2)="critical_downdraft_local"
 
  end function pdf_analysis_get_descriptor
 
 
  subroutine init_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state
    integer :: ierr, inc
 
    tpts = current_state%global_grid%size(x_index)*current_state%global_grid%size(y_index)
 
    start_x = current_state%local_grid%local_domain_start_index(x_index)
    end_x   = current_state%local_grid%local_domain_end_index(x_index)
    start_y = current_state%local_grid%local_domain_start_index(y_index)
    end_y   = current_state%local_grid%local_domain_end_index(y_index)
 
    xsize = end_x - start_x + 1
    ysize = end_y - start_y + 1
 
    lpts = xsize*ysize
 
    uppercrit  = options_get_real(current_state%options_database, "uppercrit")
    dwnpercrit = options_get_real(current_state%options_database, "dwnpercrit")
 
    show_critical_w = options_get_logical(current_state%options_database, "show_critical_w")
 
!    if (current_state%parallel%my_rank == 0) 
     allocate(tmp_all(tpts))
!    else
!      allocate(tmp_all(1))
!    end if
 
    allocate(gpts_on_proc(current_state%parallel%processes))
    call mpi_allgather(lpts, 1, mpi_int, gpts_on_proc, 1, mpi_int, current_state%parallel%monc_communicator, ierr)
 
    allocate(displacements(current_state%parallel%processes)) 
    displacements(1) = 0
    do inc = 2, current_state%parallel%processes
      displacements(inc) = displacements(inc-1) + gpts_on_proc(inc-1)
    end do ! loop over processes
 
    if (.not. current_state%continuation_run) then
      allocate(current_state%global_grid%configuration%vertical%w_dwn(current_state%local_grid%size(z_index)),&
               current_state%global_grid%configuration%vertical%w_up(current_state%local_grid%size(z_index)))
    end if
 
    ! Save the sampling_frequency to force diagnostic calculation on select time steps
    diagnostic_generation_frequency=options_get_integer(current_state%options_database, "sampling_frequency")
 
 end subroutine init_callback
 
 
  subroutine timestep_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state
 
 
    if (mod(current_state%timestep, diagnostic_generation_frequency) == 0) call calculate_w_percentiles(current_state)
 
  end subroutine timestep_callback
 
 
  subroutine finalisation_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state
 
    if (allocated(tmp_all)) deallocate(tmp_all)    
  end subroutine finalisation_callback
 
 
  subroutine calculate_w_percentiles(current_state)
    type(model_state_type), target, intent(inout) :: current_state
    real(kind=default_precision), dimension(lpts) :: tmp_var
 
    integer :: i, j, k, num_neg, num_pos, dd_thresh_pos, ud_thresh_pos
    integer :: max_up_k, min_dwn_k 
    real(kind=default_precision), dimension((lpts+1)/2) :: t
    real(kind=default_precision), dimension((tpts+1)/2) :: tall
    real(kind=default_precision)                        :: max_up, min_dwn, &
                                                           max_up_th, min_dwn_th
    integer :: ierr
    real(kind=default_precision), dimension(ysize,xsize) :: l2d ! local 2d data
 
    max_up_th  = 0.0_default_precision
    min_dwn_th = 0.0_default_precision
    max_up     = max_up_th
    min_dwn    = min_dwn_th
    max_up_k   = 0
    min_dwn_k  = 0
 
    current_state%global_grid%configuration%vertical%w_dwn(:) = 0.0_default_precision 
    current_state%global_grid%configuration%vertical%w_up(:)  = 0.0_default_precision
 
    do k = 2, current_state%local_grid%size(z_index)
 
       l2d(:,:)=0.5_default_precision*(   current_state%w%data(k  , start_y:end_y, start_x:end_x)   &
                                        + current_state%w%data(k-1, start_y:end_y, start_x:end_x)    )
 
       tmp_var=pack(l2d,.true.)
 
       call mergesort(tmp_var,lpts,t)
 
       call mpi_gatherv(tmp_var, lpts, precision_type, tmp_all, gpts_on_proc, displacements, precision_type, &
                        0, current_state%parallel%monc_communicator, ierr )
 
       if (current_state%parallel%my_rank == 0) then
 
         call mergesort(tmp_all,tpts,tall)
 
         num_neg = count(tmp_all < 0.0_default_precision)
         num_pos = count(tmp_all > 0.0_default_precision)
 
         dd_thresh_pos = int(num_neg * dwnpercrit) 
         ud_thresh_pos = tpts - int(num_pos * uppercrit) + 1
 
         if ( dd_thresh_pos == 0 ) dd_thresh_pos = 1
         if ( ud_thresh_pos == 0 .or. num_pos == 0 ) ud_thresh_pos = tpts
 
         current_state%global_grid%configuration%vertical%w_dwn(k) = tmp_all(dd_thresh_pos)
         current_state%global_grid%configuration%vertical%w_up(k)  = tmp_all(ud_thresh_pos)
 
         if (show_critical_w) then
           if ( tmp_all(dd_thresh_pos) < min_dwn_th ) then
             min_dwn_th = tmp_all(dd_thresh_pos)
             min_dwn = tmp_all(1)  ! sorted array
             min_dwn_k = k
           end if 
           if ( tmp_all(ud_thresh_pos) > max_up_th ) then
             max_up_th = tmp_all(ud_thresh_pos)
             max_up = tmp_all(tpts)  ! sorted array
             max_up_k = k
           end if
         end if ! show_critical_w
 
       end if ! global operations section
 
    end do ! loop over k
 
    call mpi_bcast(current_state%global_grid%configuration%vertical%w_dwn(:), current_state%local_grid%size(z_index), &
                   precision_type, 0, current_state%parallel%monc_communicator, ierr)
    call mpi_bcast(current_state%global_grid%configuration%vertical%w_up(:),  current_state%local_grid%size(z_index), &
                   precision_type, 0, current_state%parallel%monc_communicator, ierr)
 
 
    if (show_critical_w) then
      call log_master_log(log_info, 'Time:  '//trim(conv_to_string(current_state%time))//' s')
      call log_master_log(log_info, 'Maximum updraft threshold:   '&
                          //trim(conv_to_string(max_up_th))//' found at level '//trim(conv_to_string(max_up_k)) )
      call log_master_log(log_info, 'Maximum updraft:   '&
                          //trim(conv_to_string(max_up))//' at level '//trim(conv_to_string(max_up_k)) )
      call log_master_log(log_info, 'Minimum downdraft threshold:   '&
                          //trim(conv_to_string(min_dwn_th))//' found at level '//trim(conv_to_string(min_dwn_k)) )   
      call log_master_log(log_info, 'Minimum downdraft:   '&
                          //trim(conv_to_string(min_dwn))//' at level '//trim(conv_to_string(min_dwn_k)) )
    end if ! show_critical_w
 
  end subroutine calculate_w_percentiles   
 
 
  !  https://rosettacode.org/wiki/Sorting_algorithms/Merge_sort#Fortran
  !  and modified to match local type and renamed to avoid confusion with intrinsic merge
  !  All parameters based on MergeSort.  No need to modify anything.
  subroutine mergesortmerge(A,NA,B,NB,C,NC)
 
    integer, intent(in) :: NA,NB,NC                              ! Normal usage: NA+NB = NC
    real(kind=default_precision), intent(in out) :: a(na)        ! B overlays C(NA+1:NC)
    real(kind=default_precision), intent(in)     :: b(nb)
    real(kind=default_precision), intent(in out) :: c(nc)
 
    integer :: I,J,K
 
    i = 1; j = 1; k = 1;
    do while(i <= na .and. j <= nb)
      if (a(i) <= b(j)) then
         c(k) = a(i)
         i = i+1
      else
         c(k) = b(j)
         j = j+1
      endif
      k = k + 1
    enddo
    do while (i <= na)
      c(k) = a(i)
      i = i + 1
      k = k + 1
    enddo
    return
 
  end subroutine mergesortmerge
 
  !  https://rosettacode.org/wiki/Sorting_algorithms/Merge_sort#Fortran
  !  and modified to match local type
  !! @A array of values to be sorted, returned sorted
  !! @N size of A
  !! @T I don't really understand T
  recursive subroutine mergesort(A,N,T)
 
    integer, intent(in) :: n
    real(kind=default_precision), dimension(N), intent(in out) :: a
    real(kind=default_precision), dimension((N+1)/2), intent (out) :: t
 
    integer :: na,nb
    real(kind=default_precision) :: v
 
    if (n < 2) return
    if (n == 2) then
      if (a(1) > a(2)) then
        v = a(1)
        a(1) = a(2)
        a(2) = v
      endif
      return
    endif      
    na=(n+1)/2
    nb=n-na
 
    call mergesort(a,na,t)
    call mergesort(a(na+1),nb,t)
 
    if (a(na) > a(na+1)) then
      t(1:na)=a(1:na)
      call mergesortmerge(t,na,a(na+1),nb,a,n)
    endif
    return
 
  end subroutine mergesort
 
 
  subroutine field_information_retrieval_callback(current_state, name, field_information)
    type(model_state_type), target, intent(inout) :: current_state
    character(len=*), intent(in) :: name
    type(component_field_information_type), intent(out) :: field_information
 
    field_information%field_type=component_array_field_type
    field_information%number_dimensions=1
    field_information%dimension_sizes(1)=current_state%local_grid%size(z_index)
    field_information%data_type=component_double_data_type
    field_information%enabled=.true.
 
  end subroutine field_information_retrieval_callback
 
 
 
  subroutine field_value_retrieval_callback(current_state, name, field_value)
    type(model_state_type), target, intent(inout) :: current_state
    character(len=*), intent(in) :: name
    type(component_field_value_type), intent(out) :: field_value
 
    if      (name .eq. "critical_updraft_local") then
      allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)), &
               source=current_state%global_grid%configuration%vertical%w_up(:))
    else if (name .eq. "critical_downdraft_local") then
      allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)), &
               source=current_state%global_grid%configuration%vertical%w_dwn(:))
    end if 
  end subroutine field_value_retrieval_callback
 
end module pdf_analysis_mod