Land-Sea mask in LFRic2UM

applications/lfricinputs/source/common/lfricinp_masked_field_adjust_type_mod.f90

@@ -39,27 +39,35 @@ module lfricinp_masked_field_adjust_type_mod
 
   ! Indices of adjusted points on the destination grid/mesh
   integer(kind=int32), allocatable :: adjusted_dst_indices_1D(:)
+  integer(kind=int32), allocatable :: adjusted_dst_indices_2D(:,:)
   ! Map that links the which single source point has been selected to replace
   ! the adjusted destination data point
   integer(kind=int32), allocatable :: adjusted_dst_to_src_map_2D(:,:)
+  integer(kind=int32), allocatable :: adjusted_dst_to_src_map_1D(:)
 
   ! Flag to check whether the adjustment type has been initialised
   logical :: initialised = .false.
 
   ! Destination mask - logical is true when point is VALID, logical false
   ! point should be ignored/masked out
   logical, allocatable :: dst_mask_1D(:)
+  logical, allocatable :: dst_mask_2D(:,:)
 
 contains
   procedure :: find_adjusted_points_src_2d_dst_1d
   procedure :: apply_masked_adjustment_src_2d_dst_1d
+  procedure :: find_adjusted_points_src_1d_dst_2d
+  procedure :: apply_masked_adjustment_src_1d_dst_2d
 
 end type lfricinp_masked_field_adjust_type
 
 contains
-!---------------------------------------------------------
+!----------------------------------------------------------------
 ! Start of type bound procedures
-!---------------------------------------------------------
+!----------------------------------------------------------------
+!----------------------------------------------------------------
+! Start of 2-dimensional source to 1-dimensional field procedures
+!----------------------------------------------------------------
 !> @brief   Finds all masked field destination points that will require
 !!          post regridding adjustment
 !> @param[in] src_mask  Mask that applies to the source grid
@@ -152,7 +160,7 @@ subroutine find_adjusted_points_src_2d_dst_1d(self, src_mask, dst_mask, weights)
 
 end subroutine find_adjusted_points_src_2d_dst_1d
 
-!---------------------------------------------------------
+!----------------------------------------------------------------
 !> @brief Applies post regridding adjustments to masked field destination points
 !> @param[in]     src   2d source data array
 !> @param[inout]  dst   1d destination array
@@ -199,6 +207,150 @@ subroutine apply_masked_adjustment_src_2d_dst_1d(self, src, dst)
 
 end subroutine apply_masked_adjustment_src_2d_dst_1d
 
-!---------------------------------------------------------
+!----------------------------------------------------------------
+! Start of 1-dimensional source to 2-dimensional field procedures
+!----------------------------------------------------------------
+!> @brief   Finds all masked field destination points that will require
+!!          post regridding adjustment
+!> @param[in] src_mask  Mask that applies to the source grid
+!> @param[in] dst_mask  Mask that applies to the destination grid
+!> @param[in] weights   Regridding weights for masked data
+subroutine find_adjusted_points_src_1d_dst_2d(self, src_mask, dst_mask, weights)
+!
+! Argument(s)
+!
+class(lfricinp_masked_field_adjust_type)       :: self
+logical,                            intent(in) :: src_mask(:)
+logical,                            intent(in) :: dst_mask(:,:)
+type(lfricinp_regrid_weights_type), intent(in) :: weights
+
+!
+! Local variables
+!
+integer(kind=int32)              :: i_x, i_y, j, l, w
+integer(kind=int32)              :: src_index, dst_index1, dst_index2
+integer(kind=int32), allocatable :: dst_point_contrb_record(:,:)
+real(kind=real64)                :: weight_value
+integer(kind=int32), parameter   :: unchecked = 0, src_mask_contrb_only = 1,   &
+                                    off_src_mask_contrb = 2
+logical                          :: l_on_src_mask, l_on_dst_mask
+
+allocate(self%dst_mask_2D, mold=dst_mask)
+self%dst_mask_2D(:,:) = dst_mask(:,:)
+
+! Initialise arrays that records whether dst points had any or no
+! contribution from on mask src points and the src point data
+! to replace the dst point data
+allocate(dst_point_contrb_record(size(dst_mask,1),size(dst_mask,2)))
+dst_point_contrb_record = unchecked
+
+! Loop over remap matrix, considering only non-zero weight elements, to
+! determine whether a dst point has contribution from any off mask src points
+do w = 1, weights%num_wgts
+  do l = 1, weights%num_links
+
+    dst_index1 = weights%dst_address_2d(l,1)
+    dst_index2 = weights%dst_address_2d(l,2)
+    l_on_dst_mask = dst_mask(dst_index1, dst_index2)
+
+    weight_value = abs(weights%remap_matrix(w,l))
+    if (weight_value > 0.0_real64 .and. l_on_dst_mask) then
+
+      src_index = weights%src_address(l)
+      l_on_src_mask= src_mask(src_index)
+
+      ! Update records on whether the dst point has contributions from only
+      ! masked src points, or some off mask source points.
+      select case (dst_point_contrb_record(dst_index1,dst_index2))
+
+        case (unchecked)
+          if (l_on_src_mask) then
+            dst_point_contrb_record(dst_index1,dst_index2) = src_mask_contrb_only
+          else
+            dst_point_contrb_record(dst_index1,dst_index2) = off_src_mask_contrb
+          end if
+
+        case (src_mask_contrb_only)
+          if ( .not. l_on_src_mask) then
+            dst_point_contrb_record(dst_index1,dst_index2) = off_src_mask_contrb
+          end if
+
+      end select
+
+    end if
+  end do
+end do
+
+! Set number of part resolved dst points.
+self%num_adjusted_points = count((dst_point_contrb_record ==                   &
+                                       off_src_mask_contrb))
+write (log_scratch_space, '(A,I0)') "Number of adjusted points = ",            &
+     self%num_adjusted_points
+call log_event(log_scratch_space, LOG_LEVEL_INFO)
+
+! Generate array of dst indices that requires post regridding masked adjustment
+allocate(self%adjusted_dst_indices_2D(self%num_adjusted_points,2))
+j = 0
+do i_y = 1, size(dst_mask,2)
+  do i_x = 1, size(dst_mask,1)
+  if (dst_point_contrb_record(i_x,i_y) == off_src_mask_contrb) then
+    j = j + 1
+    self%adjusted_dst_indices_2D(j,1) = i_x
+    self%adjusted_dst_indices_2D(j,2) = i_y
+  end if
+end do
+
+deallocate(dst_point_contrb_record)
+
+end subroutine find_adjusted_points_src_1d_dst_2d
+
+!----------------------------------------------------------------
+!> @brief Applies post regridding adjustments to masked field destination points
+!> @param[in]     src   1d source data array
+!> @param[inout]  dst   2d destination array
+subroutine apply_masked_adjustment_src_1d_dst_2d(self, src, dst)
+!
+! Uses real arrays, could be overloaded for different types,
+! precision and shapes
+!
+use lfricinp_um_parameters_mod, only: um_rmdi
+
+!
+! Argument(s)
+!
+real(kind=real64), intent(in) :: src(:)
+real(kind=real64), intent(in out) :: dst(:,:)
+class(lfricinp_masked_field_adjust_type) :: self
+!
+! Local variables
+!
+integer(kind=int32) :: i, i_x, i_y
+
+! Check if masked field adjust type has been initialised. If not
+! report a warning
+if (self%initialised) then
+
+  do i = 1, self%num_adjusted_points
+    dst(self%adjusted_dst_indices_2D(i,1),self%adjusted_dst_indices_2D(i,2)) = &
+                                         src(self%adjusted_dst_to_src_map_1D(i))
+  end do
+
+  do i_y = 1, size(dst,2)
+    do i_x = 1, size(dst,1)
+      if (.not. self%dst_mask_2D(i_x,i_y)) then
+        dst(i_x,i_y) = um_rmdi
+      end if
+  end do
+
+else
+
+  log_scratch_space = 'Masked field adjustment type not initialised.'
+  call log_event(log_scratch_space, LOG_LEVEL_ERROR)
+
+end if
+
+end subroutine apply_masked_adjustment_src_1d_dst_2d
+
+!----------------------------------------------------------------
 
 end module lfricinp_masked_field_adjust_type_mod

applications/lfricinputs/source/lfric2um.f90

@@ -19,6 +19,9 @@ program lfric2um
 use lfric2um_namelists_mod,               only: lfric2um_nl_fname,             &
                                                 lfric2um_config,               &
                                                 required_lfric_namelists
+
+use lfric2um_init_masked_field_adjustments_mod,                                &
+                                   only: lfric2um_init_masked_field_adjustments
 use lfric2um_initialise_um_mod,           only: lfric2um_initialise_um,        &
                                                 um_output_file
 use lfric2um_initialise_lfric2um_mod,     only: lfric2um_initialise_lfric2um
@@ -63,6 +66,9 @@ program lfric2um
 !==========================================================================
 ! lfric2um main loop
 !==========================================================================
+! Now initialise masked points that requires post regridding adjustments
+call lfric2um_init_masked_field_adjustments()
+
 ! Main loop over fields to be read, regridded and written to output dump
 call lfric2um_main_loop()
 

applications/lfricinputs/source/lfric2um/lfric2um_init_masked_field_adjustments_mod.f90

@@ -0,0 +1,122 @@
+! *****************************COPYRIGHT*******************************
+! (C) Crown copyright Met Office. All rights reserved.
+! For further details please refer to the file LICENCE
+! which you should have received as part of this distribution.
+! *****************************COPYRIGHT*******************************
+module lfric2um_init_masked_field_adjustments_mod
+
+! Intrinsic modules
+use, intrinsic :: iso_fortran_env, only: int32, int64, real64
+
+implicit none
+
+private
+
+public :: lfric2um_init_masked_field_adjustments
+
+contains
+
+!> @brief     Initialise the masked field adjustments ready for post-processing
+!> @details   Read land/sea masks from the UM file and LFRic ancils and use this
+!!            information to determine which points are not completely on one
+!!            or the other and so need post-processing after regridding.
+!!            Nearest neighbour algorithms are used to determine which data
+!!            should be used instead.
+subroutine lfric2um_init_masked_field_adjustments()
+
+use lfricinp_get_latlon_mod,           only: get_um_grid_coords
+use lfricinp_nearest_neighbour_mod,    only: find_nn_on_lfric_mesh
+use lfricinp_masks_mod,                only: lfricinp_init_masks,              &
+                                             lfricinp_finalise_masks,          &
+                                             um_land_mask,                     &
+                                             um_maritime_mask,                 &
+                                             lfric_land_mask,                  &
+                                             lfric_maritime_mask
+use lfricinp_stashmaster_mod,          only: stashcode_land_frac
+
+use lfric2um_masked_field_adjustments_mod, only: land_field_adjustments,       &
+                                                 maritime_field_adjustments
+use lfric2um_regrid_weights_mod,           only: get_weights
+
+implicit none
+
+! Local variables
+character(len=1),    parameter :: lfric_land_mask_grid_type = 'p'
+integer(kind=int32)            :: cell_lid_nn, idx, idy, i
+real(kind=real64)              :: lon, lat
+
+! Initialise masks
+call lfricinp_init_masks(stashcode_land_frac)
+
+!
+! Initialise land field adjustments
+!
+if (allocated(um_land_mask) .and. allocated(lfric_land_mask)) then
+
+  ! Find indices of UM NN land point values that will require adjustment to
+  ! land points on LFRic mesh
+  call land_field_adjustments%find_adjusted_points_src_1d_dst_2d(              &
+                                      src_mask=lfric_land_mask,                &
+                                      dst_mask=um_land_mask,                   &
+                                      weights=get_weights(stashcode_land_frac))
+  !
+  ! Set map from UM NN land point indices to LFRic adjusted land point indices
+  allocate(land_field_adjustments%adjusted_dst_to_src_map_1D(                  &
+                                    land_field_adjustments%num_adjusted_points))
+  do i = 1, land_field_adjustments%num_adjusted_points
+    idx = land_field_adjustments%adjusted_dst_indices_2D(i,1)
+    idy = land_field_adjustments%adjusted_dst_indices_2D(i,2)
+    call get_um_grid_coords(idx, idy, lon, lat)
+    call find_nn_on_lfric_mesh(lfric_mask=lfric_land_mask,                     &
+                               lfric_mask_grid_type=lric_land_mask_grid_type,  &
+                               lon_ref=lon,                                    &
+                               lat_ref=lat,                                    &
+                               cell_lid_nn=cell_lid_nn)
+    land_field_adjustments%adjusted_dst_to_src_map_1D(i) = cell_lid_nn
+  end do
+  !
+  ! Set initialisation flag of land field adjustment
+  land_field_adjustments%initialised = .true.
+
+end if
+
+!
+! Initialise maritime field adjustments
+!
+if (allocated(um_maritime_mask) .and. allocated(lfric_maritime_mask)) then
+
+  ! Find indices of maritime points on LFRic mesh that will require adjustment
+  ! to UM NN maritime point values
+  call maritime_field_adjustments%find_adjusted_points_src_1d_dst_2d(          &
+                                      src_mask=lfric_maritime_mask,            &
+                                      dst_mask=um_maritime_mask,               &
+                                      weights=get_weights(stashcode_land_frac))
+  !
+  ! Set map from LFRic adjusted maritime point indices to UM NN maritime point
+  ! indices
+  allocate(maritime_field_adjustments%adjusted_dst_to_src_map_2D(              &
+                              maritime_field_adjustments%num_adjusted_points,2 &
+                                                            ))
+  do i = 1, maritime_field_adjustments%num_adjusted_points
+    idx = maritime_field_adjustments%adjusted_dst_indices_2D(i,1)
+    idy = maritime_field_adjustments%adjusted_dst_indices_2D(i,2)
+    call get_um_grid_coords(idx, idy, lon, lat)
+    call find_nn_on_lfric_mesh(lfric_mask=lfric_maritime_mask,                 &
+                               lfric_mask_grid_type=lfric_land_mask_grid_type, &
+                               lon_ref=lon,                                    &
+                               lat_ref=lat,                                    &
+                               cell_lid_nn=cell_lid_nn)
+    maritime_field_adjustments%adjusted_dst_to_src_map_1D(i) = cell_lid_nn
+  end do
+  !
+  ! Set initialisation flag of maritime field adjustment
+  maritime_field_adjustments%initialised = .true.
+
+end if
+
+! Finalise masks
+call lfricinp_finalise_masks
+
+end subroutine lfric2um_init_masked_field_adjustments
+
+end module lfric2um_init_masked_field_adjustments_mod

applications/lfricinputs/source/lfric2um/lfric2um_initialise_um_mod.f90

@@ -193,7 +193,7 @@ subroutine lfric2um_set_integer_constants(um_output_file, um_grid, lfric_mesh)
 ! hardcode for now - doesn't exist in lfric yet
 int_constants(ih_1_c_rho_level) = 30
 !  hardcode for now - aquaplanet
-int_constants(ih_land_points) = 0
+! int_constants(ih_land_points) = 0
 int_constants(ih_ozone_levels) = int_constants(ih_model_levels)
 int_constants(ih_convect_levels) = 0
 

applications/lfricinputs/source/lfric2um/lfric2um_masked_field_adjustments_mod.f90

@@ -0,0 +1,18 @@
+! *****************************COPYRIGHT*******************************
+! (C) Crown copyright Met Office. All rights reserved.
+! For further details please refer to the file LICENCE
+! which you should have received as part of this distribution.
+! *****************************COPYRIGHT*******************************
+module lfric2um_masked_field_adjustments_mod
+
+use lfricinp_masked_field_adjust_type_mod, &
+                                        only: lfricinp_masked_field_adjust_type
+
+implicit none
+
+private
+
+type(lfricinp_masked_field_adjust_type), public  :: land_field_adjustments,    &
+                                                    maritime_field_adjustments
+
+end module lfric2um_masked_field_adjustments_mod