Fixed a small bug & updated with latest develop branch

src/soca/Geometry/soca_geom_mod.F90

@@ -21,7 +21,7 @@ module soca_geom_mod
                        restore_state, query_initialized, &
                        free_restart_type, save_restart
 use mpp_domains_mod, only : mpp_get_compute_domain, mpp_get_data_domain, &
-                            mpp_update_domains
+                            mpp_get_global_domain, mpp_update_domains
 use fms_mod,         only : write_data, read_data
 use fms_io_mod,      only : fms_io_init, fms_io_exit
 use fckit_geometry_module, only: sphere_distance
@@ -39,8 +39,11 @@ module soca_geom_mod
     integer :: nzo, nzi, nzs, ncat
     integer :: isc, iec, jsc, jec  !< indices of compute domain
     integer :: isd, ied, jsd, jed  !< indices of data domain
+    integer :: isg, ieg, jsg, jeg  !< indices of global domain
     integer :: iscl, iecl, jscl, jecl  !< indices of local compute domain
     integer :: isdl, iedl, jsdl, jedl  !< indices of local data domain
+    real(kind=kind_real), allocatable, dimension(:)   :: lonh, lath
+    real(kind=kind_real), allocatable, dimension(:)   :: lonq, latq   
     real(kind=kind_real), allocatable, dimension(:,:) :: lon, lat !< Tracer point grid
     real(kind=kind_real), allocatable, dimension(:,:) :: lonu, latu !< Zonal velocity grid
     real(kind=kind_real), allocatable, dimension(:,:) :: lonv, latv !< Meridional velocity grid
@@ -51,6 +54,7 @@ module soca_geom_mod
     real(kind=kind_real), allocatable, dimension(:,:) :: mask2dv   !< v        "   . 0 = land 1 = ocean
     real(kind=kind_real), allocatable, dimension(:,:) :: cell_area
     real(kind=kind_real), allocatable, dimension(:,:) :: rossby_radius
+    real(kind=kind_real), allocatable, dimension(:,:,:) :: h
     logical :: save_local_domain = .false. ! If true, save the local geometry for each pe.
     character(len=:), allocatable :: geom_grid_file
     type(fckit_mpi_comm) :: f_comm
@@ -139,6 +143,10 @@ end subroutine geom_init
 subroutine geom_end(self)
   class(soca_geom), intent(out)  :: self
 
+  if (allocated(self%lonh))          deallocate(self%lonh)
+  if (allocated(self%lath))          deallocate(self%lath)
+  if (allocated(self%lonq))          deallocate(self%lonq)
+  if (allocated(self%latq))          deallocate(self%latq)          
   if (allocated(self%lon))           deallocate(self%lon)
   if (allocated(self%lat))           deallocate(self%lat)
   if (allocated(self%lonu))          deallocate(self%lonu)
@@ -152,6 +160,7 @@ subroutine geom_end(self)
   if (allocated(self%mask2dv))       deallocate(self%mask2dv)
   if (allocated(self%cell_area))     deallocate(self%cell_area)
   if (allocated(self%rossby_radius)) deallocate(self%rossby_radius)
+  if (allocated(self%h))             deallocate(self%h)
   nullify(self%Domain)
 
 end subroutine geom_end
@@ -177,6 +186,10 @@ subroutine geom_clone(self, other)
 
   ! Allocate and clone geometry
   call geom_allocate(other)
+  other%lonh = self%lonh
+  other%lath = self%lath
+  other%lonq = self%lonq
+  other%latq = self%latq  
   other%lon = self%lon
   other%lat = self%lat
   other%lonu = self%lonu
@@ -190,6 +203,7 @@ subroutine geom_clone(self, other)
   other%mask2dv = self%mask2dv
   other%cell_area = self%cell_area
   other%rossby_radius = self%rossby_radius
+  other%h = self%h
 
 end subroutine geom_clone
 
@@ -202,6 +216,10 @@ subroutine geom_gridgen(self)
 
   ! Generate grid
   call soca_mom6_init(mom6_config, partial_init=.true.)
+  self%lonh = mom6_config%grid%gridlont
+  self%lath = mom6_config%grid%gridlatt
+  self%lonq = mom6_config%grid%gridlonb
+  self%latq = mom6_config%grid%gridlatb  
   self%lon = mom6_config%grid%GeoLonT
   self%lat = mom6_config%grid%GeoLatT
   self%lonu = mom6_config%grid%geoLonCu
@@ -216,6 +234,7 @@ subroutine geom_gridgen(self)
   self%mask2du = mom6_config%grid%mask2dCu
   self%mask2dv = mom6_config%grid%mask2dCv
   self%cell_area  = mom6_config%grid%areaT
+  self%h = mom6_config%MOM_CSp%h
 
   ! Get Rossby Radius
   call geom_rossby_radius(self)
@@ -237,11 +256,16 @@ subroutine geom_allocate(self)
   call geom_get_domain_indices(self, "compute", self%isc, self%iec, self%jsc, self%jec)
   call geom_get_domain_indices(self, "data", isd, ied, jsd, jed)
   self%isd = isd ;  self%ied = ied ; self%jsd = jsd; self%jed = jed
+  call geom_get_domain_indices(self, "global", self%isg, self%ieg, self%jsg, self%jeg)
   call geom_get_domain_indices(self, "compute", self%iscl, self%iecl, self%jscl, self%jecl, local=.true.)
   call geom_get_domain_indices(self, "data", self%isdl, self%iedl, self%jsdl, self%jedl, local=.true.)
   nzo = self%nzo
 
   ! Allocate arrays on compute domain
+  allocate(self%lonh(self%isg:self%ieg));        self%lonh = 0.0_kind_real
+  allocate(self%lath(self%jsg:self%jeg));        self%lath = 0.0_kind_real
+  allocate(self%lonq(self%isg:self%ieg));        self%lonq = 0.0_kind_real
+  allocate(self%latq(self%jsg:self%jeg));        self%latq = 0.0_kind_real  
   allocate(self%lon(isd:ied,jsd:jed));           self%lon = 0.0_kind_real
   allocate(self%lat(isd:ied,jsd:jed));           self%lat = 0.0_kind_real
   allocate(self%lonu(isd:ied,jsd:jed));          self%lonu = 0.0_kind_real
@@ -258,6 +282,7 @@ subroutine geom_allocate(self)
 
   allocate(self%cell_area(isd:ied,jsd:jed));     self%cell_area = 0.0_kind_real
   allocate(self%rossby_radius(isd:ied,jsd:jed)); self%rossby_radius = 0.0_kind_real
+  allocate(self%h(isd:ied,jsd:jed,1:nzo));       self%h = 0.0_kind_real
 
 end subroutine geom_allocate
 
@@ -314,6 +339,26 @@ subroutine geom_write(self)
 
   ! Save global domain
   call fms_io_init()
+  idr_geom = register_restart_field(geom_restart, &
+                                   &self%geom_grid_file, &
+                                   &'lonh', &
+                                   &self%lonh(:), &
+                                   domain=self%Domain%mpp_domain)
+  idr_geom = register_restart_field(geom_restart, &
+                                   &self%geom_grid_file, &
+                                   &'lath', &
+                                   &self%lath(:), &
+                                   domain=self%Domain%mpp_domain)
+  idr_geom = register_restart_field(geom_restart, &
+                                   &self%geom_grid_file, &
+                                   &'lonq', &
+                                   &self%lonq(:), &
+                                   domain=self%Domain%mpp_domain)
+  idr_geom = register_restart_field(geom_restart, &
+                                   &self%geom_grid_file, &
+                                   &'latq', &
+                                   &self%latq(:), &
+                                   domain=self%Domain%mpp_domain)  
   idr_geom = register_restart_field(geom_restart, &
                                    &self%geom_grid_file, &
                                    &'lon', &
@@ -379,6 +424,11 @@ subroutine geom_write(self)
                                    &'mask2dv', &
                                    &self%mask2dv(:,:), &
                                    domain=self%Domain%mpp_domain)
+  idr_geom = register_restart_field(geom_restart, &
+                                   &self%geom_grid_file, &
+                                   &'h', &
+                                   &self%h(:,:,:), &
+                                   domain=self%Domain%mpp_domain)                           
   call save_restart(geom_restart, directory='')
   call free_restart_type(geom_restart)
   call fms_io_exit()
@@ -407,6 +457,26 @@ subroutine geom_read(self)
   type(restart_file_type) :: geom_restart
 
   call fms_io_init()
+  idr_geom = register_restart_field(geom_restart, &
+                                   &self%geom_grid_file, &
+                                   &'lonh', &
+                                   &self%lonh(:), &
+                                   domain=self%Domain%mpp_domain)
+  idr_geom = register_restart_field(geom_restart, &
+                                   &self%geom_grid_file, &
+                                   &'lath', &
+                                   &self%lath(:), &
+                                   domain=self%Domain%mpp_domain)
+  idr_geom = register_restart_field(geom_restart, &
+                                   &self%geom_grid_file, &
+                                   &'lonq', &
+                                   &self%lonq(:), &
+                                   domain=self%Domain%mpp_domain)
+  idr_geom = register_restart_field(geom_restart, &
+                                   &self%geom_grid_file, &
+                                   &'latq', &
+                                   &self%latq(:), &
+                                   domain=self%Domain%mpp_domain)  
   idr_geom = register_restart_field(geom_restart, &
                                    &self%geom_grid_file, &
                                    &'lon', &
@@ -472,6 +542,11 @@ subroutine geom_read(self)
                                    &'mask2dv', &
                                    &self%mask2dv(:,:), &
                                    domain=self%Domain%mpp_domain)
+  idr_geom = register_restart_field(geom_restart, &
+                                   &self%geom_grid_file, &
+                                   &'h', &
+                                   &self%h(:,:,:), &
+                                   domain=self%Domain%mpp_domain)                           
   call restore_state(geom_restart, directory='')
   call free_restart_type(geom_restart)
   call fms_io_exit()
@@ -488,9 +563,11 @@ subroutine geom_get_domain_indices(self, domain_type, is, ie, js, je, local)
 
   integer :: isc, iec, jsc, jec
   integer :: isd, ied, jsd, jed
+  integer :: isg, ieg, jsg, jeg 
 
   call mpp_get_compute_domain(self%Domain%mpp_domain,isc,iec,jsc,jec)
   call mpp_get_data_domain(self%Domain%mpp_domain,isd,ied,jsd,jed)
+  call mpp_get_global_domain(self%Domain%mpp_domain, isg, ieg, jsg, jeg)
   if (present(local)) then
      isc = isc - (isd-1) ; iec = iec - (isd-1) ; ied = ied - (isd-1) ; isd = 1
      jsc = jsc - (jsd-1) ; jec = jec - (jsd-1) ; jed = jed - (jsd-1) ; jsd = 1
@@ -501,7 +578,8 @@ subroutine geom_get_domain_indices(self, domain_type, is, ie, js, je, local)
      is = isc; ie = iec; js = jsc; je = jec;
   case ("data")
      is = isd; ie = ied; js = jsd; je = jed;
-
+  case ("global")
+     is = isg; ie = ieg; js = jsg; je = jeg;
   end select
 
 end subroutine geom_get_domain_indices

src/soca/State/soca_state.interface.F90

@@ -267,7 +267,12 @@ subroutine soca_state_change_resol_c(c_key_fld,c_key_rhs) bind(c,name='soca_stat
     call soca_state_registry%get(c_key_rhs,rhs)
 
     ! TODO implement a proper change of resolution, just copying for now
-    call fld%copy(rhs)
+    if (size(fld%geom%lon,1)==size(rhs%geom%lon,1) .and. size(fld%geom%lat,2)==size(rhs%geom%lat,2) .and. &
+      fld%geom%nzo==rhs%geom%nzo ) then
+      call fld%copy(rhs)
+    else
+      call fld%convert(rhs)      
+    endif
 
 end subroutine soca_state_change_resol_c
 

src/soca/State/soca_state_mod.F90

@@ -8,6 +8,7 @@ module soca_state_mod
 use soca_geom_mod
 use soca_fields_mod
 use soca_increment_mod
+use soca_convert_stats_mod
 use oops_variables_mod
 use soca_geostrophy_mod
 use kinds, only: kind_real
@@ -29,7 +30,8 @@ module soca_state_mod
 
   ! misc
   procedure :: rotate => soca_state_rotate
-
+  procedure :: convert => soca_state_convert
+
 end type
 
 !------------------------------------------------------------------------------
@@ -235,5 +237,23 @@ subroutine soca_state_diff_incr(x1, x2, inc)
   end do
 end subroutine soca_state_diff_incr
 
+! ------------------------------------------------------------------------------
+!> example of usage of ConvertStat app
+subroutine soca_state_convert(self, rhs)
+  class(soca_state), intent(inout) :: self ! target
+  class(soca_state), intent(in)   :: rhs   ! source
+  integer :: n
+  type(soca_convertstat_type) :: convert_stat
+  type(soca_field), pointer :: field1, field2, hocn1, hocn2
+
+  call rhs%get("hocn", hocn1) ; call self%get("hocn", hocn2)
+  call convert_stat%setup(rhs%geom, self%geom, hocn1, hocn2)
+  do n = 1, size(rhs%fields)
+    field1 => rhs%fields(n) ; call self%get(trim(field1%name),field2)
+    if (field1%io_file=="ocn" .or. field1%io_file=="sfc" .or. field1%io_file=="ice")  &
+            call convert_stat%change_resol(field1, field2, rhs%geom, self%geom)
+  end do !n
+  call convert_stat%clean()
+end subroutine soca_state_convert
 
 end module

src/soca/Utils/CMakeLists.txt

@@ -1,4 +1,5 @@
 soca_target_sources(
+  soca_convert_stats_mod.F90 	
   soca_geostrophy_mod.F90
   soca_omb_stats_mod.F90
   soca_utils.F90