swapsmooth.F90

Go to the documentation of this file.

module swapsmooth_mod
  use monc_component_mod, only : component_descriptor_type
  use state_mod, only : model_state_type, forward_stepping
  use grids_mod, only : x_index, y_index, z_index
  use datadefn_mod, only : default_precision
  implicit none

#ifndef TEST_MODE
  private
#endif

  logical :: mean_profiles_active=.false. 

  public swapsmooth_get_descriptor

contains

  type(component_descriptor_type) function swapsmooth_get_descriptor()
    swapsmooth_get_descriptor%name="swap_smooth"
    swapsmooth_get_descriptor%version=0.1
    swapsmooth_get_descriptor%initialisation=>initialisation_callback
    swapsmooth_get_descriptor%timestep=>timestep_callback
  end function swapsmooth_get_descriptor

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

#ifdef U_ACTIVE
    mean_profiles_active=allocated(current_state%global_grid%configuration%vertical%olubar)
    return
#endif
#ifdef V_ACTIVE
    mean_profiles_active=allocated(current_state%global_grid%configuration%vertical%olvbar)
    return
#endif
    if (current_state%th%active) then
      mean_profiles_active=allocated(current_state%global_grid%configuration%vertical%olthbar)   
      return
    end if
    if (current_state%number_q_fields .gt. 0) then
       mean_profiles_active=allocated(current_state%global_grid%configuration%vertical%olqbar)
       return
    end if
  end subroutine initialisation_callback  

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

    if (.not. current_state%halo_column) then
      if (current_state%field_stepping == forward_stepping) then
        call swap_and_smooth_classic(current_state, .false.)
      else
        ! Centred stepping
        call swap_and_smooth_robert_filter(current_state)
      end if
    end if

    if (mean_profiles_active .and. current_state%last_timestep_column) then
      if (current_state%field_stepping == forward_stepping) then
        call classic_for_average_profiles(current_state, .false.)
      else
        call robert_filter_for_average_profiles(current_state)
      end if
    end if    
  end subroutine timestep_callback

  subroutine swap_and_smooth_classic(current_state, old_smoother)
    type(model_state_type), intent(inout) :: current_state
    logical, intent(in) :: old_smoother

    integer :: y_index, x_index, k, n
    real(kind=DEFAULT_PRECISION) :: c1, c2, existing_value

    if (old_smoother) then
      c1 = 1.0_default_precision - current_state%tsmth                                               
      c2 = 2.0_default_precision * current_state%tsmth - 1.0_default_precision
    else
      c1 = 1.0_default_precision
      c2 = -1.0_default_precision
    end if

    x_index=current_state%column_local_x
    y_index=current_state%column_local_y
    
    do k=1,current_state%global_grid%size(z_index)
#ifdef U_ACTIVE
      existing_value = current_state%u%data(k,y_index,x_index) + current_state%zu%data(k,y_index,x_index)
      current_state%u%data(k,y_index,x_index)=existing_value * c1 + current_state%u%data(k,y_index,x_index) * c2
      current_state%zu%data(k,y_index,x_index)=existing_value - current_state%u%data(k,y_index,x_index)
#endif
#ifdef V_ACTIVE          
      existing_value = current_state%v%data(k,y_index,x_index) + current_state%zv%data(k,y_index,x_index)
      current_state%v%data(k,y_index,x_index)=existing_value * c1 + current_state%v%data(k,y_index,x_index) * c2
      current_state%zv%data(k,y_index,x_index)=existing_value - current_state%v%data(k,y_index,x_index)
#endif        
#ifdef W_ACTIVE
      existing_value = current_state%w%data(k,y_index,x_index) + current_state%zw%data(k,y_index,x_index)
      current_state%w%data(k,y_index,x_index)=existing_value * c1 + current_state%w%data(k,y_index,x_index) * c2
      current_state%zw%data(k,y_index,x_index)=existing_value - current_state%w%data(k,y_index,x_index)
#endif        
      if (current_state%th%active) then
        existing_value = current_state%th%data(k,y_index,x_index) + current_state%zth%data(k,y_index,x_index)
        current_state%th%data(k,y_index,x_index)=existing_value * c1 + current_state%th%data(k,y_index,x_index) * c2
        current_state%zth%data(k,y_index,x_index)=existing_value - current_state%th%data(k,y_index,x_index)
      end if
      do n=1,current_state%number_q_fields
        if (current_state%q(n)%active) then
          existing_value = current_state%q(n)%data(k,y_index,x_index) + current_state%zq(n)%data(k,y_index,x_index)
          current_state%q(n)%data(k,y_index,x_index)=existing_value * c1 + current_state%q(n)%data(k,y_index,x_index) * c2
          current_state%zq(n)%data(k,y_index,x_index)=existing_value - current_state%q(n)%data(k,y_index,x_index)
        end if
      end do
    end do
  end subroutine swap_and_smooth_classic

  subroutine swap_and_smooth_robert_filter(current_state)
    type(model_state_type), intent(inout) :: current_state

    integer :: y_index, x_index, k, n
    real(kind=DEFAULT_PRECISION) :: c1, c2, existing_value

    x_index=current_state%column_local_x
    y_index=current_state%column_local_y

    c1 = 1.0_default_precision - 2.0_default_precision*current_state%tsmth                                               
    c2 = current_state%tsmth
   
    do k=1,current_state%global_grid%size(z_index)
#ifdef U_ACTIVE
      existing_value = current_state%u%data(k,y_index,x_index)
      current_state%u%data(k,y_index,x_index)=current_state%zu%data(k,y_index,x_index)
      current_state%zu%data(k,y_index,x_index)=c1*existing_value+c2*(current_state%u%data(k, y_index, x_index)+&
           current_state%savu%data(k,y_index,x_index) -current_state%ugal)
#endif
#ifdef V_ACTIVE   
      existing_value = current_state%v%data(k,y_index,x_index)
      current_state%v%data(k,y_index,x_index)=current_state%zv%data(k,y_index,x_index)
      current_state%zv%data(k,y_index,x_index)=c1*existing_value+c2*(current_state%v%data(k, y_index, x_index)+&
           current_state%savv%data(k,y_index,x_index)-current_state%vgal)
#endif
#ifdef W_ACTIVE
      existing_value = current_state%w%data(k,y_index,x_index)
      current_state%w%data(k,y_index,x_index)=current_state%zw%data(k,y_index,x_index)
      current_state%zw%data(k,y_index,x_index)=c1*existing_value+c2*(current_state%w%data(k, y_index, x_index)+&
           current_state%savw%data(k,y_index,x_index))
#endif
      if (current_state%th%active) then
        ! Uses the partial smooth of theta from stepfields 
        existing_value = current_state%zth%data(k,y_index,x_index)
        current_state%zth%data(k,y_index,x_index)=current_state%th%data(k,y_index,x_index) + current_state%tsmth * existing_value
        current_state%th%data(k,y_index,x_index)=existing_value
      end if
      do n=1, current_state%number_q_fields
        if (current_state%q(n)%active) then
          ! Uses the partial smooth of q from stepfields 
          existing_value = current_state%zq(n)%data(k,y_index,x_index)
          current_state%zq(n)%data(k,y_index,x_index)=current_state%q(n)%data(k,y_index,x_index)+&
               current_state%tsmth * existing_value
          current_state%q(n)%data(k,y_index,x_index)=existing_value
        end if
      end do
    end do
  end subroutine swap_and_smooth_robert_filter

  subroutine classic_for_average_profiles(current_state, old_smoother)
    type(model_state_type), intent(inout) :: current_state
    logical, intent(in) :: old_smoother

    integer :: k, n
    real(kind=DEFAULT_PRECISION) :: c1, c2

    if (old_smoother) then
      c1 = 1.0_default_precision - current_state%tsmth                                               
      c2 = 2.0_default_precision * current_state%tsmth - 1.0_default_precision
    else
      c1 = 1.0_default_precision
      c2 = -1.0_default_precision
    end if

    do k=1,current_state%global_grid%size(z_index)
#ifdef U_ACTIVE
      current_state%global_grid%configuration%vertical%olzubar(k)=current_state%global_grid%configuration%vertical%olubar(k) +&
           current_state%global_grid%configuration%vertical%olzubar(k)
      current_state%global_grid%configuration%vertical%olubar(k)=current_state%global_grid%configuration%vertical%olzubar(k) *&
           c1 + current_state%global_grid%configuration%vertical%olubar(k) * c2
      current_state%global_grid%configuration%vertical%olzubar(k)=current_state%global_grid%configuration%vertical%olzubar(k) -&
           current_state%global_grid%configuration%vertical%olubar(k)
#endif
#ifdef V_ACTIVE
      current_state%global_grid%configuration%vertical%olzvbar(k)=current_state%global_grid%configuration%vertical%olvbar(k) +&
           current_state%global_grid%configuration%vertical%olzvbar(k)
      current_state%global_grid%configuration%vertical%olvbar(k)=current_state%global_grid%configuration%vertical%olzvbar(k) *&
           c1 + current_state%global_grid%configuration%vertical%olvbar(k) * c2
      current_state%global_grid%configuration%vertical%olzvbar(k)=current_state%global_grid%configuration%vertical%olzvbar(k) -&
           current_state%global_grid%configuration%vertical%olvbar(k)
#endif
      if (current_state%th%active) then
        current_state%global_grid%configuration%vertical%olzthbar(k)=current_state%global_grid%configuration%vertical%olthbar(k)+&
             current_state%global_grid%configuration%vertical%olzthbar(k)
        current_state%global_grid%configuration%vertical%olthbar(k)=current_state%global_grid%configuration%vertical%olzthbar(k)*&
             c1 + current_state%global_grid%configuration%vertical%olthbar(k) * c2
        current_state%global_grid%configuration%vertical%olzthbar(k)=&
             current_state%global_grid%configuration%vertical%olzthbar(k)-&
             current_state%global_grid%configuration%vertical%olthbar(k)
      end if
      if (current_state%number_q_fields .gt. 0) then
        do n=1, current_state%number_q_fields
          current_state%global_grid%configuration%vertical%olzqbar(k,n)=&
               current_state%global_grid%configuration%vertical%olqbar(k,n)+&
               current_state%global_grid%configuration%vertical%olzqbar(k,n)
          current_state%global_grid%configuration%vertical%olqbar(k,n)=&
               current_state%global_grid%configuration%vertical%olzqbar(k,n)*&
               c1 + current_state%global_grid%configuration%vertical%olqbar(k,n) * c2
          current_state%global_grid%configuration%vertical%olzqbar(k,n)=&
               current_state%global_grid%configuration%vertical%olzqbar(k,n)-&
               current_state%global_grid%configuration%vertical%olqbar(k,n)
        end do
      end if      
    end do    
  end subroutine classic_for_average_profiles  

  subroutine robert_filter_for_average_profiles(current_state)
    type(model_state_type), intent(inout) :: current_state

    integer :: k, n
    real(kind=DEFAULT_PRECISION) :: c1, c2, existing_value

    c1 = 1.0_default_precision - 2.0_default_precision*current_state%tsmth                                               
    c2 = current_state%tsmth

    do k=1,current_state%global_grid%size(z_index)
#ifdef U_ACTIVE
      existing_value=current_state%global_grid%configuration%vertical%olubar(k)
      current_state%global_grid%configuration%vertical%olubar(k)=current_state%global_grid%configuration%vertical%olzubar(k)
      current_state%global_grid%configuration%vertical%olzubar(k)=c1*existing_value+c2*&
           (current_state%global_grid%configuration%vertical%olubar(k) + &
           current_state%global_grid%configuration%vertical%savolubar(k))
#endif
#ifdef V_ACTIVE
      existing_value=current_state%global_grid%configuration%vertical%olvbar(k)
      current_state%global_grid%configuration%vertical%olvbar(k)=current_state%global_grid%configuration%vertical%olzvbar(k)
      current_state%global_grid%configuration%vertical%olzvbar(k)=c1*existing_value+c2*&
           (current_state%global_grid%configuration%vertical%olvbar(k) + &
           current_state%global_grid%configuration%vertical%savolvbar(k))
#endif
      if (current_state%th%active) then
        existing_value=current_state%global_grid%configuration%vertical%olzthbar(k)
        current_state%global_grid%configuration%vertical%olzthbar(k)=&
             current_state%global_grid%configuration%vertical%olthbar(k) + current_state%tsmth * existing_value
        current_state%global_grid%configuration%vertical%olthbar(k)=existing_value
      end if
      if (current_state%number_q_fields .gt. 0) then
        do n=1, current_state%number_q_fields
          existing_value=current_state%global_grid%configuration%vertical%olzqbar(k,n)
          current_state%global_grid%configuration%vertical%olzqbar(k,n)=&
               current_state%global_grid%configuration%vertical%olqbar(k,n) + current_state%tsmth * existing_value
          current_state%global_grid%configuration%vertical%olqbar(k,n)=existing_value
        end do
      end if
    end do
  end subroutine robert_filter_for_average_profiles  
end module swapsmooth_mod