its working..

Author: mlee03

horiz_interp/horiz_interp_conserve.F90

@@ -214,11 +214,10 @@ subroutine horiz_interp_read_weights_conserve(Interp, weight_filename, weight_fi
     integer, intent(in) :: isw, iew, jsw, jew
     integer, intent(in), optional :: src_tile
 
-    integer :: i, j, ncells, mpi_ncells
-    integer :: tile_id(1) !< tile id for saving xgrid for the correct tile
-    integer :: istart(1), iend(1), i_dst, j_dst, index
-    real(8), allocatable :: dst_area2(:,:), dst_area1(:), read1(:), read1_element
-    integer, allocatable :: tile1(:), indices(:), read2(:,:)
+    integer :: i, j, ncells, domain_ncells
+    integer :: istart, iend, i_dst, j_dst, index
+    real(8), allocatable :: dst_area2(:,:), dst_area1(:), read1(:), xarea(:)
+    integer, allocatable :: tile1(:), read2(:,:)
     logical, allocatable :: mask(:)
 
     type(FmsNetcdfFile_t) :: weight_fileobj !< FMS2io fileob for the weight file
@@ -234,91 +233,91 @@ subroutine horiz_interp_read_weights_conserve(Interp, weight_filename, weight_fi
       ! get ncells
       call get_dimension_size(weight_fileobj, "ncells", ncells)
 
+      istart = 1
+      iend = ncells
+
       !get section of xgrid on src_tile
       if(present(src_tile)) then
         allocate(tile1(ncells))
         call read_data(weight_fileobj, "tile1", tile1)
-        istart = FINDLOC(tile1, src_tile)
-        iend = FINDLOC(tile1, src_tile, back=.true.)
-        ncells = iend(1) - istart(1) + 1
+        !find istart
+        do i=1, ncells
+          if(tile1(i) == src_tile) then
+            istart = i
+            exit
+          end if
+        end do
+        !find iend
+        do i=istart, ncells
+          if(tile1(i) /= src_tile) then
+            iend = i
+            exit
+          end if
+        end do
+        ncells = iend - istart + 1
         deallocate(tile1)
-      else
-        istart(1) = 1
-        iend(1) = ncells
       end if
 
       ! allocate arrays for reading data
-      allocate(read2(2, ncells), indices(ncells))
+      allocate(read2(2, ncells))
 
       ! get section of xgrid for the specified window (compute domain) on the tgt grid
-      call read_data(weight_fileobj, "tile2_cell", read2, corner=[1,istart(1)], edge_lengths=[2,iend(1)])
+      call read_data(weight_fileobj, "tile2_cell", read2, corner=[1,istart], edge_lengths=[2,ncells])
+
+      ! get xgrid indices, used copilot for this section of code
       allocate(mask(ncells))
-      mpi_ncells = 0
-      mask = (read2(1,:) >= isw) .and. (read2(1,:) <= iew) .and. (read2(2,:) >= jsw) .and. (read2(2,:) <= jew)
-      do i = 1, ncells
-        if (mask(i)) then
-          mpi_ncells = mpi_ncells + 1
-          indices(mpi_ncells) = i
-        end if
-      end do
-      deallocate(mask)
+      mask = (read2(1,:) >= isw .and. read2(1,:) <= iew .and. read2(2,:) >= jsw .and. read2(2,:) <= jew)
+
+      domain_ncells = count(mask)
 
       ! allocate data to store xgrid
-      allocate(Interp%i_src(mpi_ncells))
-      allocate(Interp%j_src(mpi_ncells))
-      allocate(Interp%i_dst(mpi_ncells))
-      allocate(Interp%j_dst(mpi_ncells))
-      allocate(Interp%horizInterpReals8_type%area_frac_dst(mpi_ncells))
-
-      !save dst parent cell indices on pe domain
-      do i=1, mpi_ncells
-        index = indices(i)
-        Interp%i_dst(i) = read2(1,index) - isw + 1
-        Interp%j_dst(i) = read2(2,index) - jsw + 1
-      end do
+      allocate(Interp%i_src(domain_ncells))
+      allocate(Interp%j_src(domain_ncells))
+      allocate(Interp%i_dst(domain_ncells))
+      allocate(Interp%j_dst(domain_ncells))
+      allocate(Interp%horizInterpReals8_type%area_frac_dst(domain_ncells))
+
+      Interp%i_dst = pack(read2(1,:), mask) - isw + 1
+      Interp%j_dst = pack(read2(2,:), mask) - jsw + 1
 
       !save src parent cell indices
-      call read_data(weight_fileobj, "tile1_cell", read2, corner=[1,istart(1)], edge_lengths=[2,iend(1)])
-      do i=1, mpi_ncells
-        index = indices(i)
-        Interp%i_src(i) = read2(1, index)
-        Interp%j_src(i) = read2(2, index)
-      end do
+      call read_data(weight_fileobj, "tile1_cell", read2, corner=[1,istart], edge_lengths=[2,ncells])
+      Interp%i_src = pack(read2(1,:), mask)
+      Interp%j_src = pack(read2(2,:), mask)
 
       deallocate(read2)
 
       ! allocate arrays to compute weights
-      allocate(read1(ncells), dst_area1(mpi_ncells), dst_area2(nlon_dst, nlat_dst))
+      allocate(read1(ncells), dst_area1(domain_ncells), dst_area2(nlon_dst, nlat_dst), xarea(domain_ncells))
 
       ! read xgrid area
-      call read_data(weight_fileobj, "xgrid_area", read1, corner=[istart(1)], edge_lengths=[iend(1)])
+      call read_data(weight_fileobj, "xgrid_area", read1, corner=[istart], edge_lengths=[ncells])
+
+      xarea = pack(read1, mask)
 
       !sum over xgrid area to get destination grid area
       dst_area2 = 0.0
-      do i = 1, mpi_ncells
-        index = indices(i)
+      do i = 1, domain_ncells
         i_dst = Interp%i_dst(i)
         j_dst = Interp%j_dst(i)
-        read1_element = read1(index)
-        dst_area2(i_dst, j_dst) =+ read1_element
+        dst_area2(i_dst, j_dst) = dst_area2(i_dst, j_dst) + xarea(i)
         dst_area1(i) = dst_area2(i_dst, j_dst)
-        Interp%horizInterpReals8_type%area_frac_dst(i) = read1_element
       end do
 
-      Interp%horizInterpReals8_type%area_frac_dst = dst_area1/Interp%horizInterpReals8_type%area_frac_dst
+      Interp%horizInterpReals8_type%area_frac_dst = xarea/dst_area1
 
       deallocate(read1)
       deallocate(dst_area1)
       deallocate(dst_area2)
-      deallocate(indices)
+      deallocate(xarea)
 
       call close_file(weight_fileobj)
 
     else
       call mpp_error(FATAL, "cannot open weight file")
     end if
 
-    Interp%nxgrid = mpi_ncells
+    Interp%nxgrid = domain_ncells
     Interp%nlon_src = nlon_src
     Interp%nlat_src = nlat_src
     Interp%nlon_dst = nlon_dst

libFMS.F90

@@ -402,7 +402,8 @@ module fms
   !> horiz_interp
   use horiz_interp_mod, only: fms_horiz_interp => horiz_interp, fms_horiz_interp_new => horiz_interp_new, &
                               fms_horiz_interp_del => horiz_interp_del, fms_horiz_interp_init => horiz_interp_init, &
-                              fms_horiz_interp_end => horiz_interp_end
+                              fms_horiz_interp_end => horiz_interp_end, &
+                              fms_horiz_interp_read_weights_conserve => horiz_interp_read_weights_conserve
   use horiz_interp_type_mod, only: FmsHorizInterp_type => horiz_interp_type, &
                               assignment(=), fms_horiz_interp_type_stats => stats
   !! used via horiz_interp

test_fms/horiz_interp/test_horiz_interp_read.F90

@@ -16,10 +16,10 @@ program main
 
   implicit none
 
-  integer, parameter :: nlon_dst = 32 !< number of lon cells in destination grid
-  integer, parameter :: nlat_dst = 16 !< number of lat cells in destionation grid
-  integer, parameter :: nlon_src = nlon_dst !< number of lon cells in src grid
-  integer, parameter :: nlat_src = nlat_dst !< number of lat cells in src_grid
+  integer, parameter :: nlon_src = 32 !< number of lon cells in src grid
+  integer, parameter :: nlat_src = 16 !< number of lat cells in src_grid
+  integer, parameter :: nlon_dst = nlon_src*2 !< number of lon cells in destination grid
+  integer, parameter :: nlat_dst = nlat_src*2 !< number of lat cells in destionation grid
   integer, parameter :: ncells = nlon_dst*nlat_dst !< number of exchange grid cells
 
   type(horiz_interp_type) :: Interp !< Interp to be populated
@@ -28,7 +28,7 @@ program main
   type(domain2d) :: domain !< domain
   integer :: global_indices(4) = [1, nlon_dst, 1, nlat_dst] !< global_indices
   integer :: layout(2) !< layout
-  integer :: isc, iec, jsc, jec !< compute domain indices on dst grid
+  integer :: isc, iec, jsc, jec, pe !< compute domain indices on dst grid
   integer :: nlon_dst_d, nlat_dst_d !< dst grid size on pe domain
 
   real(8), allocatable :: data_src(:,:), data_dst(:,:) !< data for remapping
@@ -41,34 +41,36 @@ program main
 
   ! get domain
   call mpp_define_layout(global_indices, mpp_npes(), layout)
-  call mpp_define_domains(global_indices, layout, domain, tile_id=1)
+  call mpp_define_domains(global_indices, layout, domain)
   call mpp_get_compute_domain(domain, isc, iec, jsc, jec, xsize=nlon_dst_d, ysize=nlat_dst_d)
 
+  pe = mpp_pe()
+
   ! write remap file to test
-  if(mpp_pe() == mpp_root_pe()) call write_remap(nlon_dst, nlat_dst, ncells)
+  if(pe == mpp_root_pe()) call write_remap(nlon_dst, nlat_dst, ncells)
   call mpp_sync()
 
   ! allocate src data
   allocate(data_src(nlon_src, nlat_src))
   do j=1, nlat_src
     do i=1, nlon_src
-      data_src(i,j) = j*nlon_src + i
+      data_src(i,j) = pe*10000 + j*100 + i
     end do
   end do
 
   ! allocate dst data on domain
-  allocate(data_dst(isc:iec, jsc:jec))
+  allocate(data_dst(nlon_dst_d, nlat_dst_d))
   data_dst = -99.0
 
   ! read remap file
   call horiz_interp_read_weights_conserve(Interp, remap, "fregrid", nlon_src, nlat_src, &
-    nlon_dst_d, nlat_dst_d, isc, iec, jsc, jec, 1)
+    nlon_dst_d, nlat_dst_d, isc, iec, jsc, jec, mpp_pe()+1)
 
   ! remap
   call horiz_interp(Interp, data_src, data_dst)
 
   ! check answers
-  if(.not.all(data_src(isc:iec, jsc:jec) == data_dst)) then
+  if(.not.all(data_src == data_dst)) then
     call mpp_error(FATAL, "error remapping with read-in exchange grid")
   end if
 
@@ -85,7 +87,7 @@ subroutine write_remap(nlon_dst, nlat_dst, ncells)
     integer, allocatable :: tile1(:), tile1_cell(:, :), tile2_cell(:, :)
     real(r8_kind), allocatable :: lon_dst(:,:), lat_dst(:,:), xgrid_area(:)
 
-    integer :: i, j, ij
+    integer :: i, j, ij, itile_x, itile_y, count
 
     allocate(lon_dst(nlon_dst+1, nlat_dst+1), lat_dst(nlon_dst+1, nlat_dst+1))
     ! create grid in radians
@@ -106,15 +108,24 @@ subroutine write_remap(nlon_dst, nlat_dst, ncells)
 
     ! xgrid
     ij = 1
-    do j=1, nlat_src
-      do i=1, nlon_src
-        tile1_cell(1, ij) = i
-        tile1_cell(2, ij) = j
-        ij = ij + 1
+    do itile_y=0, 1
+      do itile_x=0, 1
+        do j=1, nlat_src
+          do i=1, nlon_src
+            tile1_cell(1, ij) = i
+            tile1_cell(2, ij) = j
+            tile2_cell(1, ij) = itile_x*nlon_src + i
+            tile2_cell(2, ij) = itile_y*nlat_src + j
+            ij = ij + 1
+          end do
+        end do
       end do
     end do
-    tile2_cell = tile1_cell
-    tile1 = 1
+
+    tile1(1:ncells/4) = 1
+    tile1(ncells/4+1:ncells/2) = 2
+    tile1(ncells/2+1:ncells*3/4) = 3
+    tile1(ncells*3/4+1:ncells) = 4
 
    !write remap file
     if(open_file(fileobj, "remap.nc", "overwrite")) then