MOM_oda_driver.F90

!> Interfaces for MOM6 ensembles and data assimilation.
module mom_oda_driver_mod
 
! This file is part of MOM6. see LICENSE.md for the license.
use fms_mod, only : open_namelist_file, close_file, check_nml_error
use fms_mod, only : error_mesg, fatal
use mpp_mod, only : stdout, stdlog, mpp_error, npes=>mpp_npes,pe=>mpp_pe
use mpp_mod, only : set_current_pelist => mpp_set_current_pelist
use mpp_mod, only : set_root_pe => mpp_set_root_pe
use mpp_mod, only : mpp_sync_self, mpp_sum, get_pelist=>mpp_get_current_pelist, mpp_root_pe
use mpp_mod, only : set_stack_size=>mpp_set_stack_size, broadcast=>mpp_broadcast
use mpp_io_mod, only : io_set_stack_size=>mpp_io_set_stack_size
use mpp_io_mod, only : mpp_single,mpp_multi
use mpp_domains_mod, only : domain2d, mpp_global_field
use mpp_domains_mod, only : mpp_get_compute_domain, mpp_get_data_domain
use mpp_domains_mod, only : mpp_redistribute, mpp_broadcast_domain
use mpp_domains_mod, only : set_domains_stack_size=>mpp_domains_set_stack_size
use diag_manager_mod, only : register_diag_field, diag_axis_init, send_data
use ensemble_manager_mod, only : get_ensemble_id, get_ensemble_size
use ensemble_manager_mod, only : get_ensemble_pelist, get_ensemble_filter_pelist
use time_manager_mod, only : time_type, decrement_time, increment_time
use time_manager_mod, only : get_date, operator(>=),operator(/=),operator(==),operator(<)
use constants_mod, only : radius, epsln
! ODA Modules
use ocean_da_types_mod, only : grid_type, ocean_profile_type, ocean_control_struct
use ocean_da_core_mod, only : ocean_da_core_init, get_profiles
!use eakf_oda_mod, only : ensemble_filter
use write_ocean_obs_mod, only : open_profile_file
use write_ocean_obs_mod, only : write_profile,close_profile_file
use kdtree, only : kd_root !# JEDI
! MOM Modules
use mom_io, only : slasher, mom_read_data
use mom_diag_mediator, only : diag_ctrl, set_axes_info
use mom_error_handler, only : fatal, warning, mom_error, mom_mesg, is_root_pe
use mom_get_input, only : get_mom_input, directories
use mom_grid, only : ocean_grid_type, mom_grid_init
use mom_grid_initialize, only : set_grid_metrics
use mom_hor_index, only : hor_index_type, hor_index_init
use mom_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
use mom_transcribe_grid, only : copy_dyngrid_to_mom_grid, copy_mom_grid_to_dyngrid
use mom_fixed_initialization, only : mom_initialize_fixed, mom_initialize_topography
use mom_coord_initialization, only : mom_initialize_coord
use mom_file_parser, only : read_param, get_param, param_file_type
use mom_string_functions, only : lowercase
use mom_ale, only : ale_cs, ale_initthicknesstocoord, ale_init, ale_updateverticalgridtype
use mom_domains, only : mom_domains_init, mom_domain_type, clone_mom_domain
use mom_remapping, only : remapping_cs, initialize_remapping, remapping_core_h
use mom_regridding, only : regridding_cs, initialize_regridding
use mom_regridding, only : regridding_main, set_regrid_params
use mom_unit_scaling, only : unit_scale_type, unit_scaling_init
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type, verticalgridinit
 
implicit none ; private
 
public :: init_oda, oda_end, set_prior_tracer, get_posterior_tracer
public :: set_analysis_time, oda, save_obs_diff, apply_oda_tracer_increments
 
#include <MOM_memory.h>
 
!> Control structure that contains a transpose of the ocean state across ensemble members.
type, public :: oda_cs ; private
  type(ocean_control_struct), pointer :: ocean_prior=> null() !< ensemble ocean prior states in DA space
  type(ocean_control_struct), pointer :: ocean_posterior=> null() !< ensemble ocean posterior states
                                                                  !! or increments to prior in DA space
  integer :: nk !< number of vertical layers used for DA
  type(ocean_grid_type), pointer :: grid => null() !< MOM6 grid type and decomposition for the DA
  type(ptr_mpp_domain), pointer, dimension(:) :: domains => null() !< Pointer to mpp_domain objects
                                                                       !! for ensemble members
  type(verticalgrid_type), pointer :: gv => null() !< vertical grid for DA
  type(unit_scale_type), pointer :: &
    us => null()    !< structure containing various unit conversion factors for DA
 
  type(domain2d), pointer :: mpp_domain => null() !< Pointer to a mpp domain object for DA
  type(grid_type), pointer :: oda_grid !< local tracer grid
  real, pointer, dimension(:,:,:) :: h => null() !<layer thicknesses [H ~> m or kg m-2] for DA
  type(thermo_var_ptrs), pointer :: tv => null() !< pointer to thermodynamic variables
  integer :: ni          !< global i-direction grid size
  integer :: nj          !< global j-direction grid size
  logical :: reentrant_x !< grid is reentrant in the x direction
  logical :: reentrant_y !< grid is reentrant in the y direction
  logical :: tripolar_n !< grid is folded at its north edge
  logical :: symmetric !< Values at C-grid locations are symmetric
  integer :: assim_method !< Method: NO_ASSIM,EAKF_ASSIM or OI_ASSIM
  integer :: ensemble_size !< Size of the ensemble
  integer :: ensemble_id = 0 !< id of the current ensemble member
  integer, pointer, dimension(:,:) :: ensemble_pelist !< PE list for ensemble members
  integer, pointer, dimension(:) :: filter_pelist !< PE list for ensemble members
  integer :: assim_frequency !< analysis interval in hours
  ! Profiles local to the analysis domain
  type(ocean_profile_type), pointer :: profiles => null() !< pointer to linked list of all available profiles
  type(ocean_profile_type), pointer :: cprofiles => null()!< pointer to linked list of current profiles
  type(kd_root), pointer :: kdroot => null() !< A structure for storing nearest neighbors
  type(ale_cs), pointer :: ale_cs=>null() !< ALE control structure for DA
  logical :: use_ale_algorithm !< true is using ALE remapping
  type(regridding_cs) :: regridcs !< ALE control structure for regridding
  type(remapping_cs) :: remapcs !< ALE control structure for remapping
  type(time_type) :: time !< Current Analysis time
  type(diag_ctrl) :: diag_cs !<Diagnostics control structure
end type oda_cs
 
!> A structure with a pointer to a domain2d, to allow for the creation of arrays of pointers.
type :: ptr_mpp_domain
  type(domain2d), pointer :: mpp_domain => null() !< pointer to an mpp domain2d
end type ptr_mpp_domain
 
!>@{  DA parameters
integer, parameter :: no_assim = 0, oi_assim=1, eakf_assim=2
!!@}
 
contains
 
!> initialize First_guess (prior) and Analysis grid
!! information for all ensemble members
subroutine init_oda(Time, G, GV, CS)
 
  type(time_type), intent(in) :: time !< The current model time.
  type(ocean_grid_type), pointer :: g !< domain and grid information for ocean model
  type(verticalgrid_type), intent(in) :: gv   !< The ocean's vertical grid structure
  type(oda_cs), pointer, intent(inout) :: cs  !< The DA control structure
 
! Local variables
  type(thermo_var_ptrs) :: tv_dummy
  type(dyn_horgrid_type), pointer :: dg=> null()
  type(hor_index_type), pointer :: hi=> null()
  type(directories) :: dirs
 
  type(grid_type), pointer :: t_grid !< global tracer grid
  real, dimension(:,:), allocatable :: global2d, global2d_old
  real, dimension(:), allocatable :: lon1d, lat1d, glon1d, glat1d
  type(param_file_type) :: pf
  integer :: n, m, k, i, j, nk
  integer :: is,ie,js,je,isd,ied,jsd,jed
  integer :: stdout_unit
  character(len=32) :: assim_method
  integer :: npes_pm, ens_info(6), ni, nj
  character(len=128) :: mesg
  character(len=32) :: fldnam
  character(len=30) :: coord_mode
  character(len=200) :: inputdir, basin_file
  logical :: reentrant_x, reentrant_y, tripolar_n, symmetric
 
  if (associated(cs)) call mpp_error(fatal,'Calling oda_init with associated control structure')
  allocate(cs)
! Use ens1 parameters , this could be changed at a later time
! if it were desirable to have alternate parameters, e.g. for the grid
! for the analysis
  call get_mom_input(pf,dirs,ensemble_num=0)
 
  call unit_scaling_init(pf, cs%US)
 
  call get_param(pf, "MOM", "ASSIM_METHOD", assim_method,  &
       "String which determines the data assimilation method "//&
       "Valid methods are: \'EAKF\',\'OI\', and \'NO_ASSIM\'", default='NO_ASSIM')
  call get_param(pf, "MOM", "ASSIM_FREQUENCY", cs%assim_frequency,  &
       "data assimilation frequency in hours")
  call get_param(pf, "MOM", "USE_REGRIDDING", cs%use_ALE_algorithm , &
                "If True, use the ALE algorithm (regridding/remapping).\n"//&
                "If False, use the layered isopycnal algorithm.", default=.false. )
  call get_param(pf, "MOM", "REENTRANT_X", cs%reentrant_x, &
       "If true, the domain is zonally reentrant.", default=.true.)
  call get_param(pf, "MOM", "REENTRANT_Y", cs%reentrant_y, &
       "If true, the domain is meridionally reentrant.", &
       default=.false.)
  call get_param(pf,"MOM", "TRIPOLAR_N", cs%tripolar_N, &
       "Use tripolar connectivity at the northern edge of the "//&
       "domain.  With TRIPOLAR_N, NIGLOBAL must be even.", &
       default=.false.)
  call get_param(pf,"MOM", "NIGLOBAL", cs%ni, &
       "The total number of thickness grid points in the "//&
       "x-direction in the physical domain.")
  call get_param(pf,"MOM", "NJGLOBAL", cs%nj, &
       "The total number of thickness grid points in the "//&
       "y-direction in the physical domain.")
  call get_param(pf, 'MOM', "INPUTDIR", inputdir)
  inputdir = slasher(inputdir)
 
  select case(lowercase(trim(assim_method)))
  case('eakf')
      cs%assim_method = eakf_assim
  case('oi')
     cs%assim_method = oi_assim
  case('no_assim')
      cs%assim_method = no_assim
  case default
      call mpp_error(fatal,'Invalid assimilation method provided')
  end select
 
  ens_info = get_ensemble_size()
  cs%ensemble_size = ens_info(1)
  npes_pm=ens_info(3)
  cs%ensemble_id = get_ensemble_id()
  !! Switch to global pelist
  allocate(cs%ensemble_pelist(cs%ensemble_size,npes_pm))
  allocate(cs%filter_pelist(cs%ensemble_size*npes_pm))
  call get_ensemble_pelist(cs%ensemble_pelist,'ocean')
  call get_ensemble_filter_pelist(cs%filter_pelist,'ocean')
 
  call set_current_pelist(cs%filter_pelist)
 
  allocate(cs%domains(cs%ensemble_size))
  cs%domains(cs%ensemble_id)%mpp_domain => g%Domain%mpp_domain
  do n=1,cs%ensemble_size
    if (.not. associated(cs%domains(n)%mpp_domain)) allocate(cs%domains(n)%mpp_domain)
    call set_root_pe(cs%ensemble_pelist(n,1))
    call mpp_broadcast_domain(cs%domains(n)%mpp_domain)
  enddo
  call set_root_pe(cs%filter_pelist(1))
  allocate(cs%Grid)
  ! params NIHALO_ODA, NJHALO_ODA set the DA halo size
  call mom_domains_init(cs%Grid%Domain,pf,param_suffix='_ODA')
  allocate(hi)
  call hor_index_init(cs%Grid%Domain, hi, pf, &
                      local_indexing=.false.)  ! Use global indexing for DA
  call verticalgridinit( pf, cs%GV, cs%US )
  allocate(dg)
  call create_dyn_horgrid(dg, hi)
  call clone_mom_domain(cs%Grid%Domain, dg%Domain,symmetric=.false.)
  call set_grid_metrics(dg,pf)
  call mom_initialize_topography(dg%bathyT,dg%max_depth,dg,pf)
  call mom_initialize_coord(cs%GV, cs%US, pf, .false., &
           dirs%output_directory, tv_dummy, dg%max_depth)
  call ale_init(pf, cs%GV, cs%US, dg%max_depth, cs%ALE_CS)
  call mom_grid_init(cs%Grid, pf, global_indexing=.true.)
  call ale_updateverticalgridtype(cs%ALE_CS,cs%GV)
  call copy_dyngrid_to_mom_grid(dg, cs%Grid)
  cs%mpp_domain => cs%Grid%Domain%mpp_domain
  cs%Grid%ke = cs%GV%ke
  cs%nk = cs%GV%ke
  ! initialize storage for prior and posterior
  allocate(cs%Ocean_prior)
  call init_ocean_ensemble(cs%Ocean_prior,cs%Grid,cs%GV,cs%ensemble_size)
  allocate(cs%Ocean_posterior)
  call init_ocean_ensemble(cs%Ocean_posterior,cs%Grid,cs%GV,cs%ensemble_size)
  allocate(cs%tv)
 
  call get_param(pf, 'oda_driver', "REGRIDDING_COORDINATE_MODE", coord_mode, &
       "Coordinate mode for vertical regridding.", &
       default="ZSTAR", fail_if_missing=.false.)
  call initialize_regridding(cs%regridCS, cs%GV, cs%US, dg%max_depth,pf,'oda_driver',coord_mode,'','')
  call initialize_remapping(cs%remapCS,'PLM')
  call set_regrid_params(cs%regridCS, min_thickness=0.)
  call mpp_get_data_domain(g%Domain%mpp_domain,isd,ied,jsd,jed)
  if (.not. associated(cs%h)) then
    allocate(cs%h(isd:ied,jsd:jed,cs%GV%ke)); cs%h(:,:,:)=0.0
    ! assign thicknesses
    call ale_initthicknesstocoord(cs%ALE_CS,g,cs%GV,cs%h)
  endif
  allocate(cs%tv%T(isd:ied,jsd:jed,cs%GV%ke)); cs%tv%T(:,:,:)=0.0
  allocate(cs%tv%S(isd:ied,jsd:jed,cs%GV%ke)); cs%tv%S(:,:,:)=0.0
 
  call set_axes_info(cs%Grid, cs%GV, cs%US, pf, cs%diag_cs, set_vertical=.true.)
  do n=1,cs%ensemble_size
    write(fldnam,'(a,i2.2)') 'temp_prior_',n
    cs%Ocean_prior%id_t(n)=register_diag_field('ODA',trim(fldnam),cs%diag_cs%axesTL%handles,time, &
                                               'ocean potential temperature','degC')
    write(fldnam,'(a,i2.2)') 'salt_prior_',n
    cs%Ocean_prior%id_s(n)=register_diag_field('ODA',trim(fldnam),cs%diag_cs%axesTL%handles,time, &
                                               'ocean salinity','psu')
    write(fldnam,'(a,i2.2)') 'temp_posterior_',n
    cs%Ocean_posterior%id_t(n)=register_diag_field('ODA',trim(fldnam),cs%diag_cs%axesTL%handles,time, &
                                               'ocean potential temperature','degC')
    write(fldnam,'(a,i2.2)') 'salt_posterior_',n
    cs%Ocean_posterior%id_s(n)=register_diag_field('ODA',trim(fldnam),cs%diag_cs%axesTL%handles,time, &
                                               'ocean salinity','psu')
  enddo
 
  call mpp_get_data_domain(cs%mpp_domain,isd,ied,jsd,jed)
  allocate(cs%oda_grid)
  cs%oda_grid%x => cs%Grid%geolonT
  cs%oda_grid%y => cs%Grid%geolatT
 
  call get_param(pf, 'oda_driver', "BASIN_FILE", basin_file, &
          "A file in which to find the basin masks, in variable 'basin'.", &
          default="basin.nc")
  basin_file = trim(inputdir) // trim(basin_file)
  allocate(cs%oda_grid%basin_mask(isd:ied,jsd:jed))
  cs%oda_grid%basin_mask(:,:) = 0.0
  call mom_read_data(basin_file,'basin',cs%oda_grid%basin_mask,cs%Grid%domain, timelevel=1)
 
!    get global grid information from ocean_model
  allocate(t_grid)
  allocate(t_grid%x(cs%ni,cs%nj))
  allocate(t_grid%y(cs%ni,cs%nj))
  allocate(t_grid%basin_mask(cs%ni,cs%nj))
  call mpp_global_field(cs%mpp_domain, cs%Grid%geolonT, t_grid%x)
  call mpp_global_field(cs%mpp_domain, cs%Grid%geolatT, t_grid%y)
  call mpp_global_field(cs%mpp_domain, cs%oda_grid%basin_mask, t_grid%basin_mask)
  t_grid%ni = cs%ni
  t_grid%nj = cs%nj
  t_grid%nk = cs%nk
  allocate(t_grid%mask(cs%ni,cs%nj,cs%nk))
  allocate(t_grid%z(cs%ni,cs%nj,cs%nk))
  allocate(global2d(cs%ni,cs%nj))
  allocate(global2d_old(cs%ni,cs%nj))
  t_grid%mask(:,:,:) = 0.0
  t_grid%z(:,:,:) = 0.0
 
  do k = 1, cs%nk
    call mpp_global_field(g%Domain%mpp_domain, cs%h(:,:,k), global2d)
    do i=1, cs%ni; do j=1, cs%nj
      if ( global2d(i,j) > 1 ) then
        t_grid%mask(i,j,k) = 1.0
      endif
    enddo ; enddo
    if (k == 1) then
      t_grid%z(:,:,k) = global2d/2
    else
      t_grid%z(:,:,k) = t_grid%z(:,:,k-1) + (global2d + global2d_old)/2
    endif
    global2d_old = global2d
  enddo
 
  call ocean_da_core_init(cs%mpp_domain, t_grid, cs%Profiles, time)
 
  cs%Time=time
  !! switch back to ensemble member pelist
  call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
end subroutine init_oda
 
!> Copy ensemble member tracers to ensemble vector.
subroutine set_prior_tracer(Time, G, GV, h, tv, CS)
  type(time_type), intent(in)    :: time !< The current model time
  type(ocean_grid_type), pointer :: g !< domain and grid information for ocean model
  type(verticalgrid_type),               intent(in)    :: gv   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h    !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),                 intent(in) :: tv   !< A structure pointing to various thermodynamic variables
 
  type(oda_cs), pointer :: cs !< ocean DA control structure
  real, dimension(:,:,:), allocatable :: t, s
  type(ocean_grid_type), pointer :: grid=>null()
  integer :: i,j, m, n, ss
  integer :: is, ie, js, je
  integer :: isc, iec, jsc, jec
  integer :: isd, ied, jsd, jed
  integer :: id
  logical :: used
 
  ! return if not time for analysis
  if (time < cs%Time) return
 
  if (.not. associated(cs%Grid)) call mom_error(fatal,'ODA_CS ensemble horizontal grid not associated')
  if (.not. associated(cs%GV)) call mom_error(fatal,'ODA_CS ensemble vertical grid not associated')
 
  !! switch to global pelist
  call set_current_pelist(cs%filter_pelist)
  call mom_mesg('Setting prior')
 
  isc=cs%Grid%isc;iec=cs%Grid%iec;jsc=cs%Grid%jsc;jec=cs%Grid%jec
  call mpp_get_compute_domain(cs%domains(cs%ensemble_id)%mpp_domain,is,ie,js,je)
  call mpp_get_data_domain(cs%domains(cs%ensemble_id)%mpp_domain,isd,ied,jsd,jed)
  allocate(t(isd:ied,jsd:jed,cs%nk))
  allocate(s(isd:ied,jsd:jed,cs%nk))
 
  do j=js,je; do i=is,ie
    call remapping_core_h(cs%remapCS, gv%ke, h(i,j,:), tv%T(i,j,:), &
         cs%nk, cs%h(i,j,:), t(i,j,:))
    call remapping_core_h(cs%remapCS, gv%ke, h(i,j,:), tv%S(i,j,:), &
         cs%nk, cs%h(i,j,:), s(i,j,:))
  enddo ; enddo
 
  do m=1,cs%ensemble_size
    call mpp_redistribute(cs%domains(m)%mpp_domain, t,&
         cs%mpp_domain, cs%Ocean_prior%T(:,:,:,m), complete=.true.)
    call mpp_redistribute(cs%domains(m)%mpp_domain, s,&
         cs%mpp_domain, cs%Ocean_prior%S(:,:,:,m), complete=.true.)
    if (cs%Ocean_prior%id_t(m)>0) &
      used=send_data(cs%Ocean_prior%id_t(m), cs%Ocean_prior%T(isc:iec,jsc:jec,:,m), cs%Time)
    if (cs%Ocean_prior%id_s(m)>0) &
      used=send_data(cs%Ocean_prior%id_s(m), cs%Ocean_prior%S(isc:iec,jsc:jec,:,m), cs%Time)
  enddo
  deallocate(t,s)
 
  !! switch back to ensemble member pelist
  call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
 
  return
 
end subroutine set_prior_tracer
 
!> Returns posterior adjustments or full state
!!Note that only those PEs associated with an ensemble member receive data
subroutine get_posterior_tracer(Time, CS, h, tv, increment)
  type(time_type), intent(in) :: time !< the current model time
  type(oda_cs), pointer :: cs !< ocean DA control structure
  real, dimension(:,:,:), pointer :: h    !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs), pointer :: tv   !< A structure pointing to various thermodynamic variables
  logical, optional, intent(in) :: increment !< True if returning increment only
 
  type(ocean_control_struct), pointer :: ocean_increment=>null()
  integer :: i, j, m
  logical :: used, get_inc
 
  ! return if not analysis time (retain pointers for h and tv)
  if (time < cs%Time) return
 
 
  !! switch to global pelist
  call set_current_pelist(cs%filter_pelist)
  call mom_mesg('Getting posterior')
 
  get_inc = .true.
  if (present(increment)) get_inc = increment
 
  if (get_inc) then
    allocate(ocean_increment)
    call init_ocean_ensemble(ocean_increment,cs%Grid,cs%GV,cs%ensemble_size)
    ocean_increment%T = cs%Ocean_posterior%T - cs%Ocean_prior%T
    ocean_increment%S = cs%Ocean_posterior%S - cs%Ocean_prior%S
  endif
  do m=1,cs%ensemble_size
    if (get_inc) then
      call mpp_redistribute(cs%mpp_domain, ocean_increment%T(:,:,:,m), &
              cs%domains(m)%mpp_domain, cs%tv%T, complete=.true.)
      call mpp_redistribute(cs%mpp_domain, ocean_increment%S(:,:,:,m), &
              cs%domains(m)%mpp_domain, cs%tv%S, complete=.true.)
    else
      call mpp_redistribute(cs%mpp_domain, cs%Ocean_posterior%T(:,:,:,m), &
              cs%domains(m)%mpp_domain, cs%tv%T, complete=.true.)
      call mpp_redistribute(cs%mpp_domain, cs%Ocean_posterior%S(:,:,:,m), &
              cs%domains(m)%mpp_domain, cs%tv%S, complete=.true.)
    endif
  enddo
 
  tv => cs%tv
  h => cs%h
  !! switch back to ensemble member pelist
  call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
 
end subroutine get_posterior_tracer
 
!> Gather observations and sall ODA routines
subroutine oda(Time, CS)
  type(time_type), intent(in) :: time !< the current model time
  type(oda_cs), intent(inout) :: cs !< the ocean DA control structure
 
  integer :: i, j
  integer :: m
  integer :: yr, mon, day, hr, min, sec
 
  if ( time >= cs%Time ) then
 
    !! switch to global pelist
    call set_current_pelist(cs%filter_pelist)
 
    call get_profiles(time, cs%Profiles, cs%CProfiles)
#ifdef ENABLE_ECDA
    call ensemble_filter(cs%Ocean_prior, cs%Ocean_posterior, cs%CProfiles, cs%kdroot, cs%mpp_domain, cs%oda_grid)
#endif
 
    !! switch back to ensemble member pelist
    call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
 
  endif
 
  return
end subroutine oda
 
!> Finalize DA module
subroutine oda_end(CS)
  type(oda_cs), intent(inout) :: cs !< the ocean DA control structure
 
end subroutine oda_end
 
!> Initialize DA module
subroutine init_ocean_ensemble(CS,Grid,GV,ens_size)
  type(ocean_control_struct), pointer :: CS !< Pointer to ODA control structure
  type(ocean_grid_type), pointer :: Grid !< Pointer to ocean analysis grid
  type(verticalgrid_type), pointer :: GV !< Pointer to DA vertical grid
  integer, intent(in) :: ens_size !< ensemble size
 
  integer :: n,is,ie,js,je,nk
 
  nk=gv%ke
  is=grid%isd;ie=grid%ied
  js=grid%jsd;je=grid%jed
  cs%ensemble_size=ens_size
  allocate(cs%T(is:ie,js:je,nk,ens_size))
  allocate(cs%S(is:ie,js:je,nk,ens_size))
  allocate(cs%SSH(is:ie,js:je,ens_size))
  allocate(cs%id_t(ens_size));cs%id_t(:)=-1
  allocate(cs%id_s(ens_size));cs%id_s(:)=-1
!  allocate(CS%U(is:ie,js:je,nk,ens_size))
!  allocate(CS%V(is:ie,js:je,nk,ens_size))
!  allocate(CS%id_u(ens_size));CS%id_u(:)=-1
!  allocate(CS%id_v(ens_size));CS%id_v(:)=-1
  allocate(cs%id_ssh(ens_size));cs%id_ssh(:)=-1
 
  return
end subroutine init_ocean_ensemble
 
!> Set the next analysis time
subroutine set_analysis_time(Time,CS)
  type(time_type), intent(in) :: time !< the current model time
  type(oda_cs), pointer, intent(inout) :: cs !< the DA control structure
 
  character(len=160) :: mesg  ! The text of an error message
  integer :: yr, mon, day, hr, min, sec
 
  if (time >= cs%Time) then
    cs%Time=increment_time(cs%Time,cs%assim_frequency*3600)
 
    call get_date(time, yr, mon, day, hr, min, sec)
    write(mesg,*) 'Model Time: ', yr, mon, day, hr, min, sec
    call mom_mesg("set_analysis_time: "//trim(mesg))
    call get_date(cs%time, yr, mon, day, hr, min, sec)
    write(mesg,*) 'Assimilation Time: ', yr, mon, day, hr, min, sec
    call mom_mesg("set_analysis_time: "//trim(mesg))
  endif
  if (cs%Time < time) then
    call mom_error(fatal, " set_analysis_time: " // &
         "assimilation interval appears to be shorter than " // &
         "the model timestep")
  endif
  return
 
end subroutine set_analysis_time
 
!> Write observation differences to a file
subroutine save_obs_diff(filename,CS)
  character(len=*), intent(in) :: filename !< name of output file
  type(oda_cs), pointer, intent(in) :: cs  !< pointer to DA control structure
 
  integer :: fid ! profile file handle
  type(ocean_profile_type), pointer :: prof=>null()
 
  fid = open_profile_file(trim(filename), nvar=2, thread=mpp_single, fset=mpp_single)
  prof=>cs%CProfiles
 
  !! switch to global pelist
  !call set_current_pelist(CS%filter_pelist)
 
  do while (associated(prof))
    call write_profile(fid,prof)
    prof=>prof%cnext
  enddo
  call close_profile_file(fid)
 
  !! switch back to ensemble member pelist
  !call set_current_pelist(CS%ensemble_pelist(CS%ensemble_id,:))
 
  return
end subroutine save_obs_diff
 
 
!> Apply increments to tracers
subroutine apply_oda_tracer_increments(dt,G,tv,h,CS)
  real,                     intent(in)    :: dt !< The tracer timestep [s]
  type(ocean_grid_type),    intent(in)    :: g  !< ocean grid structure
  type(thermo_var_ptrs),    intent(inout) :: tv !< A structure pointing to various thermodynamic variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                            intent(in)    :: h  !< layer thickness [H ~> m or kg m-2]
  type(oda_cs),             intent(inout) :: cs !< the data assimilation structure
 
end subroutine apply_oda_tracer_increments
 
!> \namespace MOM_oda_driver_mod
!!
!! \section section_ODA The Ocean data assimilation (DA) and Ensemble Framework
!!
!! The DA framework implements ensemble capability in MOM6.   Currently, this framework
!! is enabled using the cpp directive ENSEMBLE_OCEAN.  The ensembles need to be generated
!! at the level of the calling routine for oda_init or above. The ensemble instances may
!! exist on overlapping or non-overlapping processors. The ensemble information is accessed
!! via the FMS ensemble manager. An independent PE layout is used to gather (prior) ensemble
!! member information where this information is stored in the ODA control structure.  This
!! module was developed in collaboration with Feiyu Lu and Tony Rosati in the GFDL prediction
!! group for use in their coupled ensemble framework. These interfaces should be suitable for
!! interfacing MOM6 to other data assimilation packages as well.
 
end module mom_oda_driver_mod