lateral_bcs.F90

Go to the documentation of this file.

module lateral_bcs_mod
  use monc_component_mod, only : component_descriptor_type
  use state_mod, only : model_state_type
  use grids_mod, only : x_index, y_index, z_index
  use logging_mod, only : log_debug, log_get_logging_level, log_log
  use datadefn_mod, only : default_precision
  implicit none

#ifndef TEST_MODE
  private
#endif  

  INTEGER :: lbc_type_x, lbc_type_y ! probably need to add these to state

  INTEGER :: idepth, jdepth ! depth of lateral bcs in i and j directions

  INTEGER, PARAMETER :: lbc_periodic = 1
  INTEGER, PARAMETER :: lbc_rigid    = 2
  INTEGER, PARAMETER :: lbc_blend   = 3

  public lateral_bcs_get_descriptor

contains

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

  subroutine initialisation_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state
    
    ! Hard wiring these for now
    idepth=3
    jdepth=2

    lbc_type_x=lbc_rigid
    lbc_type_y=lbc_periodic

  end subroutine initialisation_callback  

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

    INTEGER :: istart, iend, jstart, jend, kstart, kend

    INTEGER :: i,j,k,iq

    LOGICAL :: is_east, is_west, is_north, is_south

    if (current_state%timestep==1) return

    is_west=current_state%parallel%wrapped_around(x_index,1) 
    is_east=current_state%parallel%wrapped_around(x_index,2) 
    is_south=current_state%parallel%wrapped_around(y_index,1) 
    is_north=current_state%parallel%wrapped_around(y_index,2) 

    ! k points
    if (is_west .or. is_east .or. is_north .or. is_south)then
      kstart=current_state%local_grid%local_domain_start_index(z_index)
      kend=current_state%local_grid%local_domain_end_index(z_index)
    end if

    ! Eastern boundary
    if (is_east)then

      iend=current_state%local_grid%local_domain_end_index(x_index) 
      istart = iend - idepth + 1
      jstart=current_state%local_grid%local_domain_start_index(y_index)
      jend=current_state%local_grid%local_domain_end_index(y_index)

      if (lbc_type_x==lbc_periodic)then
        ! how should this be done best Luis?
      else if (lbc_type_x==lbc_rigid)then
        ! No communications needed assuming idepth and jdepth fit 
        ! on a single pe (what if they don't?)
        call apply_rigid(current_state,istart,iend,jstart,jend,kstart,kend)
        ! propagate winds to halos
        do i=iend+1,iend+current_state%local_grid%halo_size(x_index)
          current_state%u%data(:,:,i) = current_state%u%data(:,:,iend)
          current_state%zu%data(:,:,i) = current_state%zu%data(:,:,iend)
          current_state%w%data(:,:,i) = current_state%w%data(:,:,iend)
          current_state%zw%data(:,:,i) = current_state%zw%data(:,:,iend)
        end do
        
      else if (lbc_type_x==lbc_blend)then
        ! To be coded....
      end if

    end if

    ! Western boundary
    if (is_west)then

      istart=current_state%local_grid%local_domain_start_index(x_index)
      iend=istart + idepth - 1
      jstart=current_state%local_grid%local_domain_start_index(y_index)
      jend=current_state%local_grid%local_domain_end_index(y_index)


      if (lbc_type_x==lbc_periodic)then
        ! how should this be done best Luis?
      else if (lbc_type_x==lbc_rigid)then
        ! No communications needed assuming idepth and jdepth fit 
        ! on a single pe (what if they don't?)
        call apply_rigid(current_state,istart,iend,jstart,jend,kstart,kend)
        ! propagate winds to halos
        do i=istart-current_state%local_grid%halo_size(x_index),istart-1
          current_state%u%data(:,:,i) = current_state%u%data(:,:,istart)
          current_state%zu%data(:,:,i) = current_state%zu%data(:,:,istart)
          current_state%w%data(:,:,i) = current_state%w%data(:,:,istart)
          current_state%zw%data(:,:,i) = current_state%zw%data(:,:,istart)
        end do
      else if (lbc_type_x==lbc_blend)then
        ! To be coded....
      end if
    end if
      
    
    ! Northern boundary
    if (is_north)then

      istart=current_state%local_grid%local_domain_start_index(x_index)
      iend=current_state%local_grid%local_domain_end_index(x_index)
      jend=current_state%local_grid%local_domain_end_index(y_index)
      jstart=jend - jdepth + 1

      if (lbc_type_y==lbc_periodic)then
        ! how should this be done best Luis?
      else if (lbc_type_y==lbc_rigid)then
        ! No communications needed assuming idepth and jdepth fit 
        ! on a single pe (what if they don't?)
        call apply_rigid(current_state,istart,iend,jstart,jend,kstart,kend)
        do j=jend+1,jend+current_state%local_grid%halo_size(y_index)
          current_state%v%data(:,:,j) = current_state%v%data(:,:,jend)
          current_state%zv%data(:,:,j) = current_state%zv%data(:,:,jend)
          current_state%w%data(:,:,j) = current_state%w%data(:,:,jend)
          current_state%zw%data(:,:,j) = current_state%zw%data(:,:,jend)
        end do
      else if (lbc_type_y==lbc_blend)then
        ! To be coded....
      end if

    end if
      

    ! Southern boundary
    if (is_south)then

      istart=current_state%local_grid%local_domain_start_index(x_index)
      iend=current_state%local_grid%local_domain_end_index(x_index)
      jstart=current_state%local_grid%local_domain_start_index(y_index)
      jend=jstart + jdepth - 1

      if (lbc_type_y==lbc_periodic)then
        ! how should this be done best Luis?
      else if (lbc_type_y==lbc_rigid)then
        ! No communications needed assuming idepth and jdepth fit 
        ! on a single pe (what if they don't?)
        call apply_rigid(current_state,istart,iend,jstart,jend,kstart,kend)
        do j=jstart-current_state%local_grid%halo_size(y_index),jstart-1
          current_state%v%data(:,:,j) = current_state%v%data(:,:,jstart)
          current_state%zv%data(:,:,j) = current_state%zv%data(:,:,jstart)
          current_state%w%data(:,:,j) = current_state%w%data(:,:,jstart)
          current_state%zw%data(:,:,j) = current_state%zw%data(:,:,jstart)
        end do
      else if (lbc_type_y==lbc_blend)then
        ! To be coded....
      end if

    end if
  end subroutine timestep_callback

  subroutine apply_rigid(current_state,istart,iend,jstart,jend,kstart,kend)
    type(model_state_type), target, intent(inout) :: current_state

    INTEGER :: istart, iend, jstart, jend, kstart, kend

    INTEGER :: i,j,k,iq

    do i=istart, iend
      do j=jstart, jend
        do k=kstart,kend
          ! Update to old values, i.e. fixed in time
#ifdef U_ACTIVE
!          current_state%u%data(k,j,i) = current_state%zu%data(k,j,i)
#endif
#ifdef V_ACTIVE
!          current_state%v%data(k,j,i) = current_state%zv%data(k,j,i)
#endif
#ifdef W_ACTIVE
!          current_state%w%data(k,j,i) = current_state%zw%data(k,j,i)
#endif
          current_state%th%data(k,j,i) = current_state%zth%data(k,j,i)
        end do
      end do
    end do
    do iq = 1, current_state%number_q_fields
      do i=istart, iend
        do j=jstart, jend
          do k=kstart,kend
            current_state%q(iq)%data(k,j,i) = current_state%zq(iq)%data(k,j,i) 
          end do
        end do
      end do
    end do

  end subroutine apply_rigid

end module lateral_bcs_mod