mom_oda_driver_mod Module Reference

Detailed Description

Interfaces for MOM6 ensembles and data assimilation.

Data Types
type	oda_cs
	Control structure that contains a transpose of the ocean state across ensemble members. More...
type	ptr_mpp_domain
	A structure with a pointer to a domain2d, to allow for the creation of arrays of pointers. More...

Functions/Subroutines
subroutine, public	init_oda (Time, G, GV, CS)
	initialize First_guess (prior) and Analysis grid information for all ensemble members More...
subroutine, public	set_prior_tracer (Time, G, GV, h, tv, CS)
	Copy ensemble member tracers to ensemble vector. More...
subroutine, public	get_posterior_tracer (Time, CS, h, tv, increment)
	Returns posterior adjustments or full state Note that only those PEs associated with an ensemble member receive data. More...
subroutine, public	oda (Time, CS)
	Gather observations and sall ODA routines. More...
subroutine, public	oda_end (CS)
	Finalize DA module. More...
subroutine	init_ocean_ensemble (CS, Grid, GV, ens_size)
	Initialize DA module. More...
subroutine, public	set_analysis_time (Time, CS)
	Set the next analysis time. More...
subroutine, public	save_obs_diff (filename, CS)
	Write observation differences to a file. More...
subroutine, public	apply_oda_tracer_increments (dt, G, tv, h, CS)
	Apply increments to tracers. More...

Function/Subroutine Documentation

apply_oda_tracer_increments()

subroutine, public mom_oda_driver_mod::apply_oda_tracer_increments	(	real, intent(in)	dt,
		type(ocean_grid_type), intent(in)	G,
		type(thermo_var_ptrs), intent(inout)	tv,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		type(oda_cs), intent(inout)	CS
	)

Apply increments to tracers.

Parameters

[in]	dt	The tracer timestep [s]
[in]	g	ocean grid structure
[in,out]	tv	A structure pointing to various thermodynamic variables
[in]	h	layer thickness [H ~> m or kg m-2]
[in,out]	cs	the data assimilation structure

Definition at line 529 of file MOM_oda_driver.F90.

  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
 

get_posterior_tracer()

subroutine, public mom_oda_driver_mod::get_posterior_tracer	(	type(time_type), intent(in)	Time,
		type(oda_cs), pointer	CS,
		real, dimension(:,:,:), pointer	h,
		type(thermo_var_ptrs), pointer	tv,
		logical, intent(in), optional	increment
	)

Returns posterior adjustments or full state Note that only those PEs associated with an ensemble member receive data.

Parameters

[in]	time	the current model time
	cs	ocean DA control structure
	h	Layer thicknesses [H ~> m or kg m-2]
	tv	A structure pointing to various thermodynamic variables
[in]	increment	True if returning increment only

Definition at line 365 of file MOM_oda_driver.F90.

  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,:))
 

init_ocean_ensemble()

subroutine mom_oda_driver_mod::init_ocean_ensemble ( type(ocean_control_struct), pointer  CS, type(ocean_grid_type), pointer  Grid, type(verticalgrid_type), pointer  GV, integer, intent(in)  ens_size  )

private

Initialize DA module.

Parameters

	cs	Pointer to ODA control structure
	grid	Pointer to ocean analysis grid
	gv	Pointer to DA vertical grid
[in]	ens_size	ensemble size

Definition at line 448 of file MOM_oda_driver.F90.

  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

init_oda()

subroutine, public mom_oda_driver_mod::init_oda	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), pointer	G,
		type(verticalgrid_type), intent(in)	GV,
		type(oda_cs), intent(inout), pointer	CS
	)

initialize First_guess (prior) and Analysis grid information for all ensemble members

Parameters

[in]	time	The current model time.
	g	domain and grid information for ocean model
[in]	gv	The ocean's vertical grid structure
[in,out]	cs	The DA control structure

Definition at line 115 of file MOM_oda_driver.F90.

 
  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%US)
  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.)
 
  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,:))

oda()

subroutine, public mom_oda_driver_mod::oda	(	type(time_type), intent(in)	Time,
		type(oda_cs), intent(inout)	CS
	)

Gather observations and sall ODA routines.

Parameters

[in]	time	the current model time
[in,out]	cs	the ocean DA control structure

Definition at line 415 of file MOM_oda_driver.F90.

  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

oda_end()

subroutine, public mom_oda_driver_mod::oda_end

(

type(oda_cs), intent(inout)

CS

)

Finalize DA module.

Parameters

[in,out]

cs

the ocean DA control structure

Definition at line 442 of file MOM_oda_driver.F90.

  type(ODA_CS), intent(inout) :: CS !< the ocean DA control structure
 

save_obs_diff()

subroutine, public mom_oda_driver_mod::save_obs_diff	(	character(len=*), intent(in)	filename,
		type(oda_cs), intent(in), pointer	CS
	)

Write observation differences to a file.

Parameters

[in]	filename	name of output file
[in]	cs	pointer to DA control structure

Definition at line 502 of file MOM_oda_driver.F90.

  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

set_analysis_time()

subroutine, public mom_oda_driver_mod::set_analysis_time	(	type(time_type), intent(in)	Time,
		type(oda_cs), intent(inout), pointer	CS
	)

Set the next analysis time.

Parameters

[in]	time	the current model time
[in,out]	cs	the DA control structure

Definition at line 475 of file MOM_oda_driver.F90.

  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
 

set_prior_tracer()

subroutine, public mom_oda_driver_mod::set_prior_tracer	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), pointer	G,
		type(verticalgrid_type), intent(in)	GV,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		type(thermo_var_ptrs), intent(in)	tv,
		type(oda_cs), pointer	CS
	)

Copy ensemble member tracers to ensemble vector.

Parameters

[in]	time	The current model time
	g	domain and grid information for ocean model
[in]	gv	The ocean's vertical grid structure
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	tv	A structure pointing to various thermodynamic variables
	cs	ocean DA control structure

Definition at line 308 of file MOM_oda_driver.F90.

  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.)
  enddo
  deallocate(t,s)
 
  !! switch back to ensemble member pelist
  call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
 
  return