Remove usage of Fortran copy of field data, part 2

src/soca/AnalyticInit/soca_analytic_mod.F90

@@ -6,6 +6,7 @@
 !> Analytic state/geoval initialization module used for testing only
 module soca_analytic_mod
 
+use atlas_module, only : atlas_field
 use kinds, only : kind_real
 use ufo_geovals_mod, only : ufo_geovals
 use ufo_sampled_locations_mod, only : ufo_sampled_locations
@@ -57,23 +58,30 @@ subroutine soca_analytic_geovals(geovals, locs)
 !!
 !! \see soca_analytic_val
 subroutine soca_analytic_state(state)
-    class(soca_state), intent(inout) :: state
+  class(soca_state), intent(inout) :: state
 
-    integer :: i, j, f, k
-    real(kind=kind_real) :: val
+  integer :: i, j, f, k, idx
 
-    do f = 1, size(state%fields)
-        do i = state%geom%isd, state%geom%ied
-            do j = state%geom%jsd, state%geom%jed
-                do k = 1, state%fields(f)%nz
-                    val = soca_analytic_val(&
-                        state%fields(f)%name, state%fields(f)%lat(i,j), &
-                        state%fields(f)%lon(i,j), k*1.0_kind_real)
-                    state%fields(f)%val(i,j,k) = val
-                end do
-            end do
+  type(atlas_field) :: field
+  real(kind=kind_real), pointer :: fdata(:,:)
+
+  do f=1, size(state%fields)
+    field = state%afieldset%field(state%fields(f)%name)
+    call field%data(fdata)
+
+    do j = state%geom%jsc, state%geom%jec
+      do i = state%geom%isc, state%geom%iec
+        idx = state%geom%atlas_ij2idx(i,j)
+        do k = 1, field%shape(1)
+          fdata(k,idx) = soca_analytic_val(&
+            state%fields(f)%name, state%fields(f)%lat(i,j), &
+            state%fields(f)%lon(i,j), k*1.0_kind_real)
         end do
+      end do
     end do
+    call field%set_dirty()
+    call field%final()
+  end do
 end subroutine
 
 

src/soca/Fields/soca_fields_mod.F90

@@ -95,12 +95,6 @@ module soca_fields_mod
   !>\copybrief soca_field_check_congruent \see soca_field_check_congruent
   procedure :: check_congruent => soca_field_check_congruent
 
-  !>\copybrief soca_field_update_halo \see soca_field_update_halo
-  procedure :: update_halo     => soca_field_update_halo
-
-  !>\copybrief soca_field_stencil_interp \see soca_field_stencil_interp
-  procedure :: stencil_interp  => soca_field_stencil_interp
-
   !>\copybrief soca_field_fill_masked \see soca_field_fill_masked
   procedure :: fill_masked     => soca_field_fill_masked
 
@@ -130,9 +124,6 @@ module soca_fields_mod
   !> \copybrief soca_fields_create \see soca_fields_create
   procedure :: create => soca_fields_create
 
-  !> \copybrief soca_fields_copy \see soca_fields_copy
-  procedure :: copy   => soca_fields_copy
-
   !> \copybrief soca_fields_delete \see soca_fields_delete
   procedure :: delete => soca_fields_delete
 
@@ -184,9 +175,6 @@ module soca_fields_mod
   !> \name misc
   !! \{
 
-  !> \copybrief soca_fields_update_halos \see soca_fields_update_halos
-  procedure :: update_halos => soca_fields_update_halos
-
   !> \copybrief soca_fields_tohpoints \see soca_fields_tohpoints
   procedure :: tohpoints  => soca_fields_tohpoints
   !> \}
@@ -236,18 +224,6 @@ subroutine soca_field_copy(self, rhs)
 end subroutine soca_field_copy
 
 
-! ------------------------------------------------------------------------------
-!> Update the data in the halo region of the field.
-!!
-!! \relates soca_fields_mod::soca_field
-!! \todo have field keep a pointer to its relevant sections of soca_geom?
-subroutine soca_field_update_halo(self, geom)
-  class(soca_field),     intent(inout) :: self
-  type(soca_geom), pointer, intent(in) :: geom !< soca_geom from soca_fields
-
-  call mpp_update_domains(self%val, geom%Domain%mpp_domain)
-end subroutine soca_field_update_halo
-
 ! ------------------------------------------------------------------------------
 !> Make sure the two fields are the same in terms of name, size, shape.
 !!
@@ -273,8 +249,8 @@ end subroutine soca_field_check_congruent
 !!
 !! Interpolation used is inverse distance weidghted, taking into
 !! consideration the mask and using at most 6 neighbors.
-subroutine soca_field_stencil_interp(self, geom, fromto)
-  class(soca_field), intent(inout) :: self
+subroutine soca_field_stencil_interp(field, geom, fromto)
+  real(kind=kind_real), allocatable, intent(inout) :: field(:,:,:)
   class(soca_geom),    intent(in) :: geom   !< geometry
   character(len=4),     intent(in) :: fromto !< "u2h", "v2h"
 
@@ -289,7 +265,7 @@ subroutine soca_field_stencil_interp(self, geom, fromto)
   real(kind=kind_real), allocatable :: masksrc_local(:,:), maskdst_local(:,:)
 
   ! Initialize temporary arrays
-  allocate(val_tmp, mold=self%val)
+  allocate(val_tmp, mold=field)
   val_tmp = 0_kind_real
 
   ! Identify source and destination grids
@@ -319,7 +295,7 @@ subroutine soca_field_stencil_interp(self, geom, fromto)
   end select
 
   ! Interpolate
-  allocate(val(6,self%nz))
+  allocate(val(6,size(field, 3)))
   do j = geom%jsc, geom%jec
      do i = geom%isc, geom%iec
         ! destination on land, skip
@@ -334,12 +310,12 @@ subroutine soca_field_stencil_interp(self, geom, fromto)
            if (masksrc_local(ij(1,sti), ij(2,sti)) == 0_kind_real) cycle
 
            ! outcroping of layers, skip
-           if (abs(self%val(ij(1,sti), ij(2,sti),1)) > val_max) cycle
+           if (abs(field(ij(1,sti), ij(2,sti),1)) > val_max) cycle
 
            ! store the valid neighbors
            lon_src(nn) = lonsrc_local(ij(1,sti), ij(2,sti))
            lat_src(nn) = latsrc_local(ij(1,sti), ij(2,sti))
-           val(nn,:) = self%val(ij(1,sti), ij(2,sti),:)
+           val(nn,:) = field(ij(1,sti), ij(2,sti),:)
            nn = nn + 1
         end do
         nn = nn - 1
@@ -352,7 +328,7 @@ subroutine soca_field_stencil_interp(self, geom, fromto)
         end if
      end do
   end do
-  self%val = val_tmp
+  field = val_tmp
 
 end subroutine soca_field_stencil_interp
 
@@ -512,42 +488,6 @@ subroutine soca_fields_delete(self)
 end subroutine
 
 
-! ------------------------------------------------------------------------------
-!> Copy the contents of \p rhs to \p self.
-!!
-!! \p self will be initialized with the variable names in \p rhs if
-!! not already initialized.
-!!
-!! \see soca_fields_init_vars
-!! \relates soca_fields_mod::soca_fields
-subroutine soca_fields_copy(self, rhs)
-  class(soca_fields), intent(inout) :: self
-  class(soca_fields),  intent(in)    :: rhs !< fields to copy from
-
-  character(len=:), allocatable :: vars_str(:)
-  integer :: i
-  type(soca_field), pointer :: rhs_fld
-
-  ! initialize the variables based on the names in rhs
-  if (.not. allocated(self%fields)) then
-    self%geom => rhs%geom
-    allocate(character(len=1024) :: vars_str(size(rhs%fields)))
-    do i=1, size(vars_str)
-      vars_str(i) = rhs%fields(i)%name
-    end do
-    call soca_fields_init_vars(self, vars_str)
-  end if
-
-  ! copy values from rhs to self, only if the variable exists
-  !  in self
-  do i=1,size(self%fields)
-    call rhs%get(self%fields(i)%name, rhs_fld)
-    call self%fields(i)%copy(rhs_fld)
-  end do
-
-end subroutine
-
-
 ! ------------------------------------------------------------------------------
 !> Get a pointer to the soca_field with the given name.
 !!
@@ -595,32 +535,28 @@ function soca_fields_has(self, name) result(res)
 end function
 
 
-! ------------------------------------------------------------------------------
-!> Update the halo region of all fields.
-!!
-!! \relates soca_fields_mod::soca_fields
-subroutine soca_fields_update_halos(self)
-  class(soca_fields), intent(inout) :: self
-  integer :: i
-
-  do i=1,size(self%fields)
-    call self%fields(i)%update_halo(self%geom)
-  end do
-end subroutine soca_fields_update_halos
-
-
 ! ------------------------------------------------------------------------------
 !> Set the value of all fields to one.
 !!
 !! \relates soca_fields_mod::soca_fields
 subroutine soca_fields_ones(self)
   class(soca_fields), intent(inout) :: self
+  type(atlas_field) :: field
+  real(kind=kind_real), pointer :: fdata(:,:)
   integer :: i
 
+  ! TODO delete
   do i = 1, size(self%fields)
     self%fields(i)%val = 1.0_kind_real
   end do
 
+  do i = 1, self%afieldset%size()
+    field = self%afieldset%field(i)
+    call field%data(fdata)
+    fdata(:,:) = 1.0_kind_real
+    call field%final()
+  end do
+
 end subroutine soca_fields_ones
 
 
@@ -630,12 +566,22 @@ end subroutine soca_fields_ones
 !! \relates soca_fields_mod::soca_fields
 subroutine soca_fields_zeros(self)
   class(soca_fields), intent(inout) :: self
+  type(atlas_field) :: field
+  real(kind=kind_real), pointer :: fdata(:,:)
   integer :: i
 
+  ! TODO delete
   do i = 1, size(self%fields)
     self%fields(i)%val = 0.0_kind_real
   end do
 
+  do i = 1, self%afieldset%size()
+    field = self%afieldset%field(i)
+    call field%data(fdata)
+    fdata(:,:) = 0.0_kind_real
+    call field%final()
+  end do
+
 end subroutine soca_fields_zeros
 
 
@@ -1259,30 +1205,57 @@ end subroutine soca_fields_write_rst
 subroutine soca_fields_tohpoints(self)
   class(soca_fields), intent(inout) :: self !< self
 
-  integer :: i
-  real(kind=kind_real), allocatable :: val(:,:,:)
-  real(kind=kind_real), pointer :: lon_out(:,:) => null()
-  real(kind=kind_real), pointer :: lat_out(:,:) => null()
+  integer :: i,j,k,n,idx
   character(len=4) :: fromto
 
-  ! Associate lon_out and lat_out with the h-grid
-  lon_out => self%geom%lon
-  lat_out => self%geom%lat
+  type(atlas_field) :: afield
+  real(kind=kind_real), pointer :: adata(:,:)
+  real(kind=kind_real), allocatable :: fdata(:,:,:)
 
   ! Apply interpolation to all fields, when necessary
-  do i=1,size(self%fields)
+  do n=1,size(self%fields)
     ! Check if already on h-points
-    if (self%fields(i)%metadata%grid == 'h') cycle
+    if (self%fields(n)%metadata%grid == 'h') cycle
 
     ! Interpolate to different location of the stencil
-    fromto = self%fields(i)%metadata%grid//'toh'
-    call self%fields(i)%stencil_interp(self%geom, fromto)
-    call self%fields(i)%update_halo(self%geom)
+    fromto = self%fields(n)%metadata%grid//'toh'
+
+    ! convert from atlas to 3d fortran field...
+    ! because I don't want to fully refactor stencil interpolation
+    allocate(fdata(self%geom%isd:self%geom%ied, self%geom%jsd:self%geom%jed, self%fields(n)%nz))
+    afield = self%aFieldset%field(self%fields(n)%name)
+    call afield%data(adata)
+    do j=self%geom%jsc, self%geom%jec
+      do i=self%geom%isc, self%geom%iec
+        idx = self%geom%atlas_ij2idx(i,j)
+        do k=1,afield%shape(1)
+          fdata(i,j,k) = adata(k, idx)
+        end do
+      end do
+    end do
+    call mpp_update_domains(fdata, self%geom%Domain%mpp_domain)
+
+    ! interp
+    call soca_field_stencil_interp(fdata, self%geom, fromto)
+
+    !copy back to atlas
+    do j=self%geom%jsc, self%geom%jec
+      do i=self%geom%isc, self%geom%iec
+        idx = self%geom%atlas_ij2idx(i,j)
+        do k=1,afield%shape(1)
+          adata(k, idx) = fdata(i,j,k)
+        end do
+      end do
+    end do
+    deallocate(fdata)
+
+    call afield%set_dirty()
+    call afield%final()
 
     ! Update grid location to h-points
-    self%fields(i)%metadata%grid = 'h'
-    self%fields(i)%lon => self%geom%lon
-    self%fields(i)%lat => self%geom%lat
+    self%fields(n)%metadata%grid = 'h'
+    self%fields(n)%lon => self%geom%lon
+    self%fields(n)%lat => self%geom%lat
  end do
 
 end subroutine soca_fields_tohpoints
@@ -1471,11 +1444,11 @@ subroutine soca_fields_sync_from_atlas(self)
         field%val(i,j,:) = real_ptr(:, self%geom%atlas_ij2idx(i,j))
       end do
     end do
-    call field%update_halo(self%geom)
+    call mpp_update_domains(field%val, self%geom%Domain%mpp_domain)
+
   call afield%final()
   end do
 
-
   call self%update_metadata()
 end subroutine