oda_core_ecda.F90

Go to the documentation of this file.

! -*- f90 -*-
!***********************************************************************
!*                   GNU Lesser General Public License
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* FMS is free software: you can redistribute it and/or modify it under
!* the terms of the GNU Lesser General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or (at
!* your option) any later version.
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
!* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
!* for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************
module oda_core_ecda_mod
  ! FMS Shared modules
  use fms_mod, only : file_exist, read_data
  use fms_mod, only : open_namelist_file, check_nml_error, close_file
  use fms_mod, only : error_mesg, fatal, note
#ifdef INTERNAL_FILE_NML
  USE mpp_mod, ONLY: input_nml_file
#endif
  use mpp_mod, only : mpp_sum, stdout, stdlog, mpp_sync_self
  use mpp_mod, only : mpp_pe, mpp_root_pe
  use mpp_io_mod, only : mpp_open, mpp_close, mpp_ascii, mpp_rdonly, mpp_multi, mpp_single, mpp_netcdf
  use mpp_io_mod, only : mpp_get_atts, mpp_get_info, mpp_get_fields, mpp_read, axistype, fieldtype, mpp_get_axes
  use mpp_io_mod, only : mpp_get_axis_data, mpp_get_field_name
  use mpp_domains_mod, only : mpp_get_compute_domain, mpp_get_data_domain
  use mpp_domains_mod, only : domain2d, mpp_get_global_domain, mpp_update_domains
  use mpp_domains_mod, only : mpp_global_field
  use mpp_memutils_mod, only : mpp_print_memuse_stats
  use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time
  use time_manager_mod, only : operator( <= ), operator( - ), operator( > ), operator ( < )
  use get_cal_time_mod, only : get_cal_time
  use axis_utils_mod, only : frac_index
  use horiz_interp_type_mod, only: horiz_interp_type
  use horiz_interp_bilinear_mod, only : horiz_interp_bilinear_new
  use constants_mod, only : deg_to_rad

  ! ODA_tools modules
  use oda_types_mod, only : ocean_profile_type, ocn_obs_flag_type, grid_type, obs_clim_type
  use oda_types_mod, only : drop_profiler, mooring, satellite, drifter, ship, temp_id, salt_id, missing_value
  use oda_types_mod, only : unknown, tao
  use xbt_adjust, only : xbt_drop_rate_adjust

  implicit none

  private

  public :: copy_obs, oda_core_init, open_profile_dataset, &
       get_obs, get_obs_woa05t, get_obs_woa05s, get_obs_sst, get_obs_suv, &
       get_obs_eta, open_profile_dataset_sst, ocn_obs, ssh_td, max_profiles

  ! Parameters
  integer, parameter :: profile_file = 1
  integer, parameter :: sfc_file = 2

  ! oda_core_nml variables
  real :: max_misfit = 5.0 !< used to inflate observation errors where the difference from the first guess is large
  real :: ass_start_lat = -87.0 !< set obs domain
  real :: ass_end_lat = 87.0 !< set obs domain
  integer :: max_profiles = 50000
  namelist /oda_core_nml/ max_misfit, ass_start_lat, ass_end_lat, max_profiles

  ! Shared ocean_obs_nml namelist variables
  real :: eta_obs_start_lat = -80.0 !< set obs domain
  real :: eta_obs_end_lat = 85.0 !< set obs domain
  real :: sst_obs_start_lat = -82.0 !< set obs domain
  real :: sst_obs_end_lat = 89.0 !< set obs domain

  integer :: max_prflvs = 200 ! for vd test

  type(ocean_profile_type), target, dimension(:), allocatable :: profiles

  integer :: num_profiles, no_sst, no_prf, no_temp, no_salt, no_suv, no_eta ! total number of observations
  integer :: no_woa05

  integer :: isc, iec, jsc, jec, isd, ied, jsd, jed  ! indices for local domain on model grid
  integer :: isg, ieg, jsg, jeg
  integer :: isd_filt, ied_filt, jsd_filt, jed_filt
  integer :: isd_flt0, ied_flt0, jsd_flt0, jed_flt0
  integer :: nk

  real, dimension(:,:), allocatable, save :: mask_tao

  ! sst obs grid information
  real, allocatable :: woa05_lon(:), woa05_lat(:), woa05_z(:)
  real, allocatable :: sst_lon(:), sst_lat(:), obs_sst(:,:)
  real, allocatable, save :: obs_woa05t(:,:,:), obs_woa05s(:,:,:)
  integer ::  nlon, nlat, nlev
  integer ::  nlon_woa, nlat_woa, nlev_woa

  ! time window for DROP, MOORING and SATELLITE data respectively

  type(time_type) , dimension(0:100), public :: time_window

  type(grid_type), pointer :: grd

  type(horiz_interp_type), save :: interp

  real, allocatable, dimension(:, :) :: x_grid, y_grid, x_grid_uv, y_grid_uv
  real :: lon_out(1, 1), lat_out(1, 1)

  type(ocn_obs_flag_type) :: ocn_obs

  integer :: ssh_td

    type obs_entry_type
       character(len=128) :: filename
       character(len=16)  :: file_type
    end type obs_entry_type


contains

  subroutine init_observations(time_s, time_e, filt_domain, localize)
    type(time_type), intent(in) :: time_s, time_e
    type(domain2d), intent(in) :: filt_domain
    logical, intent(in), optional :: localize

    integer, parameter :: SUV_ID = 4, eta_id = 5, woat_id = 11, woas_id = 12

    ! ocean_obs_nml variables
    integer :: mooring_window = 5
    integer :: satellite_window = 10
    integer :: drop_window = 30
    integer :: drifter_window = 30
    integer :: ship_window = 30
    integer :: unknown_window = 30

    logical :: prfs_obs, salt_obs, sst_obs, eta_obs, suv_obs
    logical :: temp_obs_argo, salt_obs_argo, temp_obs_gtspp
    logical :: temp_obs_woa05, salt_obs_woa05
    integer :: eta_obs_td = 10
    integer :: max_files = 30
    integer :: max_files_argo = 10
    integer :: max_files_gtspp = 10
    namelist /ocean_obs_nml/ mooring_window, satellite_window, drop_window,&
         & drifter_window, ship_window, unknown_window,&
         & prfs_obs, salt_obs, sst_obs, eta_obs, suv_obs,&
         & temp_obs_argo, salt_obs_argo, temp_obs_gtspp,&
         & temp_obs_woa05, salt_obs_woa05, eta_obs_td,&
         & sst_obs_start_lat, sst_obs_end_lat, eta_obs_start_lat, eta_obs_end_lat,&
         & max_files, max_files_argo, max_files_gtspp

    integer :: i, j, n, obs_variable
    integer :: ioun, io_status, ierr
    integer :: stdout_unit, stdlog_unit
    integer :: nfiles, nrecs, unit
    integer :: nfiles_argo, nrecs_argo, unit_argo
    integer :: nfiles_gtspp, nrecs_gtspp, unit_gtspp
    integer, dimension(:), allocatable :: filetype
    integer, dimension(:), allocatable :: filetype_argo
    integer, dimension(:), allocatable :: filetype_gtspp

    character(len=128) :: input_files_woa05t, input_files_woa05s
    character(len=256) :: record
    character(len=128), dimension(:), allocatable :: input_files
    character(len=128), dimension(:), allocatable :: input_files_argo
    character(len=128), dimension(:), allocatable :: input_files_gtspp

    type(obs_entry_type) :: tbl_entry

    stdout_unit = stdout()
    stdlog_unit = stdlog()

#ifdef INTERNAL_FILE_NML
    read (input_nml_file, ocean_obs_nml, iostat=io_status)
#else
    ioun = open_namelist_file()
    read(unit=ioun, nml=ocean_obs_nml, iostat=io_status)
    ierr = check_nml_error(io_status,'ocean_obs_nml')
    call close_file(ioun)
#endif
    write (unit=stdlog_unit, nml=ocean_obs_nml)

    ! Allocate filetype* and input_files* variables
    allocate(filetype(max_files), input_files(max_files))
    allocate(filetype_argo(max_files_argo), input_files_argo(max_files_argo))
    allocate(filetype_gtspp(max_files_gtspp), input_files_gtspp(max_files_gtspp))

    filetype = -1
    filetype_argo = -1
    filetype_gtspp = -1
    input_files = ''
    input_files_argo = ''
    input_files_gtspp = ''

    if ( prfs_obs .or. salt_obs .or. temp_obs_argo .or. temp_obs_gtspp .or. salt_obs_argo ) then
       ocn_obs%use_prf_as_obs = .true.
    end if
    ocn_obs%use_sst_as_obs = sst_obs
    ocn_obs%use_ssh_as_obs = eta_obs
    ocn_obs%use_suv_as_obs = suv_obs
    ocn_obs%use_woa05_t = temp_obs_woa05
    ocn_obs%use_woa05_s = salt_obs_woa05
    ssh_td = eta_obs_td

    ! time window for DROP, MOORING and SATELLITE data respectively
    ! will be available from namelist
    time_window(:) = set_time(0,unknown_window)
    time_window(drop_profiler:drop_profiler+9) = set_time(0,drop_window)
    time_window(mooring:mooring+9) = set_time(0,mooring_window)
    time_window(satellite:satellite+9) = set_time(0,satellite_window)
    time_window(drifter:drifter+9) = set_time(0,drifter_window)
    time_window(ship:ship+9) = set_time(0,ship_window)

    nfiles = 0
    nrecs=0
    call mpp_open(unit, 'ocean_obs_table', action=mpp_rdonly)
    read_obs: do while ( nfiles <= max_files )
       read (unit=unit, fmt='(A)', iostat=io_status) record
       if ( io_status < 0 ) then
          exit read_obs
       else if ( io_status > 0 ) then
          cycle read_obs
       else
          nrecs = nrecs + 1
          if ( record(1:1) == '#' ) cycle read_obs
          read ( unit=record, fmt=*, iostat=io_status ) tbl_entry
          if ( io_status < 0 ) then
             exit read_obs
          else if ( io_status > 0 ) then
             cycle read_obs
          else
             nfiles = nfiles + 1
             input_files(nfiles) = tbl_entry%filename
             select case ( trim(tbl_entry%file_type) )
             case ('profiles')
                filetype(nfiles) = profile_file
             case ('sfc')
                filetype(nfiles) = sfc_file
             case default
                call error_mesg('oda_core_mod::init_observations', 'error in obs_table entry format', fatal)
             end select
          end if
       end if
    end do read_obs
    if ( nfiles > max_files ) then
       call error_mesg('oda_core_mod::init_observations', 'number of obs files exceeds max_files parameter', fatal)
    end if
    CALL mpp_close(unit)

    nfiles_argo = 0
    nrecs_argo = 0
    call mpp_open(unit_argo, 'ocean_obs_argo_table', action=mpp_rdonly)
    read_obs_argo: do while ( nfiles_argo <= max_files_argo )
       read (unit=unit_argo, fmt='(A)', iostat=io_status) record
       if ( io_status < 0 ) then
          exit read_obs_argo
       else if ( io_status > 0 ) then
          cycle read_obs_argo
       else
          nrecs_argo = nrecs_argo + 1
          if ( record(1:1) == '#' ) cycle read_obs_argo
          read (unit=record, fmt=*, iostat=io_status) tbl_entry
          if ( io_status < 0 ) then
             exit read_obs_argo
          else if ( io_status > 0 ) then
             cycle read_obs_argo
          else
             nfiles_argo = nfiles_argo + 1
             input_files_argo(nfiles_argo) = tbl_entry%filename
             select case ( trim(tbl_entry%file_type) )
             case ('profiles')
                filetype_argo(nfiles_argo) = profile_file
             case ('sfc')
                filetype_argo(nfiles_argo) = sfc_file
             case default
                call error_mesg('oda_core_mod::init_observations', 'error in obs_table entry format for argo', fatal)
             end select
          end if
       end if
    end do read_obs_argo
    if ( nfiles_argo > max_files_argo ) then
       call error_mesg('oda_core_mod::init_observations', 'number of obs files exceeds max_files_argo parameter', fatal)
    end if
    call mpp_close(unit_argo)

    nfiles_gtspp = 0
    nrecs_gtspp = 0
    call mpp_open(unit_gtspp, 'ocean_obs_gtspp_table', action=mpp_rdonly)
    read_obs_gtspp: do while ( nfiles_gtspp <= max_files_gtspp )
       read (unit=unit_gtspp, fmt='(A)', iostat=io_status) record
       if ( io_status < 0 ) then
          exit read_obs_gtspp
       else if ( io_status > 0 ) then
          cycle read_obs_gtspp
       else
          nrecs_gtspp = nrecs_gtspp + 1
          if ( record(1:1) == '#' ) cycle read_obs_gtspp
          read (unit=record, fmt=*, iostat=io_status) tbl_entry
          if ( io_status < 0 ) then
             exit read_obs_gtspp
          else if ( io_status > 0 ) then
             cycle read_obs_gtspp
          else
             nfiles_gtspp = nfiles_gtspp + 1
             input_files_gtspp(nfiles_gtspp) = tbl_entry%filename
             select case ( trim(tbl_entry%file_type) )
             case ('profiles')
                filetype_gtspp(nfiles_gtspp) = profile_file
             case ('sfc')
                filetype_gtspp(nfiles_gtspp) = sfc_file
             case default
                call error_mesg('oda_core_mod::init_observations', 'error in obs_table entry format for gtspp', fatal)
             end select
          end if
       end if
    end do read_obs_gtspp
    if ( nfiles_gtspp > max_files_gtspp ) then
       call error_mesg('oda_core_mod::init_observations', 'number of obs files exceeds max_files_gtspp parameter', fatal)
    end if
    CALL mpp_close(unit_gtspp)

    num_profiles = 0
    no_prf = 0
    no_sst = 0
    no_temp = 0
    no_salt = 0
    no_suv = 0
    no_eta = 0
    no_woa05 = 0

    ! get local indices for Model grid
    allocate(x_grid(isg:ieg,jsg:jeg), x_grid_uv(isg:ieg,jsg:jeg))
    allocate(y_grid(isg:ieg,jsg:jeg), y_grid_uv(isg:ieg,jsg:jeg))

    call mpp_global_field(filt_domain, grd%x(:,:), x_grid(:,:))
    call mpp_global_field(filt_domain, grd%y(:,:), y_grid(:,:))

    ! Allocate profiles
    allocate(profiles(max_profiles))

    do j=jsg, jeg
       do i=isg, ieg
          if ( x_grid(i,j) .lt. 80.0 ) x_grid(i,j) = x_grid(i,j) + 360.0
       end do
    end do

    ! uv grid may not be precise, need to be carefully checked
    x_grid_uv(:,:) = x_grid(:,:) + 0.5
    do j=jsg, jeg-1
       do i=isg, ieg
          y_grid_uv(i,j) = y_grid(i,j) + 0.5*(y_grid(i,j+1)-y_grid(i,j))
       end do
    end do
    do i=isg, ieg
       y_grid_uv(i,jeg) = 90.0
    end do

    if ( prfs_obs ) then
       obs_variable = temp_id
       if ( mpp_pe() == mpp_root_pe() ) then
          write (unit=stdout_unit, fmt='("TEMP_ID = ",I5)') temp_id
       end if
       do n=1, nfiles
          select case ( filetype(n) )
          case (profile_file)
             call open_profile_dataset(trim(input_files(n)), time_s, time_e, obs_variable, localize)
          case default
             call error_mesg('oda_core_mod::init_observations', 'filetype not currently supported for prfs_obs', fatal)
          end select
       end do
    end if

    if ( salt_obs ) then
       obs_variable = salt_id
       if ( mpp_pe() == mpp_root_pe() ) then
          write (unit=stdout_unit, fmt='("SALT_ID = ",I5)') salt_id
       end if
       do n=1, nfiles
          select case ( filetype(n) )
          case (profile_file)
             call open_profile_dataset(trim(input_files(n)), time_s, time_e, obs_variable, localize)
          case default
             call error_mesg('oda_core_mod::init_observations', 'filetype not currently supported for salt_obs', fatal)
          end select
       end do
    end if

    if ( temp_obs_gtspp ) then
       obs_variable = temp_id
       if ( mpp_pe() == mpp_root_pe() ) then
          write (unit=stdout_unit, fmt='("TEMP_ID = ",I5)') temp_id
       end if
       do n=1, nfiles_gtspp
          if ( mpp_pe() == mpp_root_pe() ) then
             write (unit=stdout_unit, fmt='("f_typ_gtspp = ",I8)') filetype_gtspp(n)
             write (unit=stdout_unit, fmt='("i_f_gtspp = ",A)') input_files_gtspp(n)
          end if
          select case ( filetype_gtspp(n) )
          case (profile_file)
             call open_profile_dataset_gtspp()
          case default
             call error_mesg('oda_core_mod::init_observations', 'filetype_gtspp not currently supported', fatal)
          end select
       end do
    end if

    if ( temp_obs_argo ) then
       obs_variable = temp_id
       if ( mpp_pe() == mpp_root_pe() ) then
          write (unit=stdout_unit, fmt='("TEMP_ID = ",I5)') temp_id
       end if
       do n=1, nfiles_argo
          if ( mpp_pe() == mpp_root_pe() ) then
             write (unit=stdout_unit, fmt='("f_typ_argo = ",I8)') filetype_argo(n)
             write (unit=stdout_unit, fmt='("i_f_argo = ",A)') input_files_argo(n)
          end if
          select case ( filetype_argo(n) )
          case (profile_file)
             call open_profile_dataset_argo(trim(input_files_argo(n)), time_s, time_e, obs_variable, localize)
          case default
             call error_mesg('oda_core_mod::init_observations', 'filetype_argo not currently supported', fatal)
          end select
       end do
    end if

    if ( salt_obs_argo ) then
       obs_variable = salt_id
       if ( mpp_pe() == mpp_root_pe() ) then
          write (unit=stdout_unit, fmt='("SALT_ID = ",I5)') salt_id
       end if
       do n=1, nfiles_argo
          select case ( filetype_argo(n) )
          case (profile_file)
             call open_profile_dataset_argo(trim(input_files_argo(n)), time_s, time_e, obs_variable, localize)
          case default
             call error_mesg('oda_core_mod::init_observations', 'filetype_argo not currently supported', fatal)
          end select
       end do
    end if

    if ( temp_obs_woa05 ) then
       obs_variable = woat_id
       if ( mpp_pe() == mpp_root_pe() ) then
          write (unit=stdout_unit, fmt='("WOAT_ID = ",I5)') woat_id
       end if
       input_files_woa05t = "INPUT/woa05_temp.nc"
       call open_profile_dataset_woa05t(trim(input_files_woa05t), obs_variable, localize)
    end if

    if ( salt_obs_woa05 ) then
       obs_variable = woas_id
       if ( mpp_pe() == mpp_root_pe() ) then
          write (unit=stdout_unit, fmt='("WOAS_ID = ",I5)') woas_id
       end if
       input_files_woa05s = "INPUT/woa05_salt.nc"
       call open_profile_dataset_woa05s(trim(input_files_woa05s), obs_variable, localize)
    end if

    if ( sst_obs ) then
       obs_variable = temp_id
       if ( mpp_pe() == mpp_root_pe() ) then
          write (unit=stdout_unit, fmt='("TEMP_ID for sst = ",I5)') temp_id
       end if
       nfiles = 1
       nrecs = 0
       input_files(nfiles) = "INPUT/sst_daily.nc"
       filetype(nfiles) = profile_file
!!$       call open_profile_dataset_sst(trim(input_files(nfiles)), obs_variable, localize)
    end if

    if ( eta_obs ) then
       obs_variable = eta_id
       nfiles = 1
       nrecs = 0
       input_files(nfiles) = "INPUT/ocean.19760101-20001231.eta_t.nc"
       filetype(nfiles) = profile_file
       call open_profile_dataset_eta(trim(input_files(nfiles)), obs_variable, localize)
    end if

    if ( suv_obs ) then
       obs_variable = suv_id
       nfiles = 1
       nrecs = 0
       input_files(nfiles) = "INPUT/sfc_current.197601-200012.nc"
       filetype(nfiles) = profile_file
       call open_profile_dataset_suv(trim(input_files(nfiles)), obs_variable, localize)
    end if

    ! Deallocate before exiting routine
    deallocate(filetype, input_files)
    deallocate(filetype_argo, input_files_argo)
    deallocate(filetype_gtspp, input_files_gtspp)
  end subroutine init_observations

  subroutine open_profile_dataset(filename, time_start, time_end, obs_variable, localize)
    character(len=*), intent(in) :: filename
    type(time_type), intent(in) :: time_start, time_end
    integer, intent(in) :: obs_variable
    logical, intent(in), optional :: localize

    integer, parameter :: max_levels = 1000
    integer, parameter :: max_lnks = 500

    real :: lon, lat, time, profile_error, rlink, flag_t, flag_s, fix_depth
    real :: ri0, rj0
    real, dimension(MAX_LEVELS) :: depth, data, t_flag, s_flag
    real, dimension(MAX_LNKS, MAX_LEVELS) :: data_bfr, depth_bfr, t_flag_bfr, s_flag_bfr

    integer :: unit, ndim, nvar, natt, nstation
    integer :: stdout_unit
    integer :: inst_type, var_id
    integer :: num_levs, k, kk, i, i0, j0, k0, nlevs, a, nn, ii, nlinks
    integer :: nprof_in_filt_domain
    integer :: bad_point, bad_point_g, out_bound_point

    logical :: data_is_local, localize_data, cont
    logical :: data_in_period
    logical :: prof_in_filt_domain
    logical, dimension(MAX_LEVELS) :: flag
    logical, dimension(MAX_LNKS, MAX_LEVELS) :: flag_bfr

    character(len=32) :: fldname, axisname, anal_fldname, time_units
    character(len=138) :: emsg_local

    type(time_type) :: profile_time
    type(axistype), pointer :: depth_axis, station_axis
    type(axistype), allocatable, dimension(:), target :: axes
    type(fieldtype), allocatable, dimension(:), target :: fields
    type(fieldtype), pointer :: field_lon, field_lat, field_flag, field_time, field_depth, field_data
    type(fieldtype), pointer :: field_error, field_link, field_t_flag, field_s_flag, field_fix_depth ! snz drop rate

    ! NOTE: fields are restricted to be in separate files

    if ( PRESENT(localize) ) then
       localize_data = localize
    else
       localize_data = .true.
    end if

    nprof_in_filt_domain = 0
    stdout_unit = stdout()

    anal_fldname = 'none'
    var_id=-1
    if ( obs_variable == temp_id ) then
       anal_fldname = 'temp'
       var_id = temp_id
    else if ( obs_variable == salt_id ) then
       anal_fldname = 'salt'
       var_id = salt_id
    end if

!    call mpp_print_memuse_stats('open_profile_dataset Start')

    call mpp_open(unit, filename, form=mpp_netcdf, fileset=mpp_single, threading=mpp_multi, action=mpp_rdonly)
    call mpp_get_info(unit, ndim, nvar, natt, nstation)

    if ( mpp_pe() .eq. mpp_root_pe() ) then
       write (unit=stdout_unit, fmt='("Opened profile dataset: ",A)') trim(filename)
    end if

    ! get axis information
    allocate(axes(ndim))
    call mpp_get_axes(unit, axes)
    do i=1, ndim
       call mpp_get_atts(axes(i), name=axisname)
       select case ( trim(axisname) )
       case ('depth_index')
          depth_axis => axes(i)
       case ('station_index')
          station_axis => axes(i)
       end select
    end do

    ! get field information
    allocate(fields(nvar))
    call mpp_get_fields(unit, fields)
    do i=1, nvar
       call mpp_get_atts(fields(i), name=fldname)
       if( var_id .eq. temp_id ) then
          select case (trim(fldname))
          case ('longitude')
             field_lon => fields(i)
          case ('latitude')
             field_lat => fields(i)
          case ('profile_flag')
             field_flag => fields(i)
          case ('time')
             field_time => fields(i)
          case ('temp')
             field_data => fields(i)
          case ('depth')
             field_depth => fields(i)
          case ('link')
             field_link => fields(i)
          case ('temp_error')
             field_error => fields(i)
          case ('temp_flag')
             field_t_flag => fields(i)
          case ('fix_depth') ! snz drop rate
             field_fix_depth => fields(i)
          end select
       else if( var_id .eq. salt_id ) then
          select case (trim(fldname))
          case ('longitude')
             field_lon => fields(i)
          case ('latitude')
             field_lat => fields(i)
          case ('profile_flag_s')
             field_flag => fields(i)
          case ('time')
             field_time => fields(i)
          case ('salt')
             field_data => fields(i)
          case ('depth')
             field_depth => fields(i)
          case ('link')
             field_link => fields(i)
          case ('salt_error')
             field_error => fields(i)
          case ('salt_flag')
             field_s_flag => fields(i)
          case ('fix_depth') ! snz drop rate
             field_fix_depth => fields(i)
          end select
       end if
    end do

    call mpp_get_atts(depth_axis, len=nlevs)

    if ( nlevs > max_levels ) then
       call error_mesg('oda_core_mod::open_profile_dataset', 'increase parameter MAX_LEVELS', fatal)
    else if (nlevs < 1) then
       call error_mesg('oda_core_mod::open_profile_dataset', 'Value of nlevs is less than 1.', fatal)
    end if

    if ( .NOT.ASSOCIATED(field_data) ) then
       call error_mesg('oda_core_mod::open_profile_dataset',&
            & 'profile dataset not used because data not needed for Analysis', note)
       return
    end if

    write(unit=stdout_unit, fmt='("There are ",I8," records in this dataset.")') nstation
    write(unit=stdout_unit, fmt='("Searching for profiles . . .")')

    call mpp_get_atts(field_time, units=time_units)

    bad_point = 0
    out_bound_point = 0

    i = 1
    cont = .true.

    do while ( cont )
       prof_in_filt_domain = .false.
       depth = missing_value  ! snz add
       data = missing_value   ! snz add

       call mpp_read(unit, field_lon, lon, tindex=i)
       call mpp_read(unit, field_lat, lat, tindex=i)
       call mpp_read(unit, field_time, time, tindex=i)
       call mpp_read(unit, field_depth, depth(1:nlevs), tindex=i)
       call mpp_read(unit, field_data, data(1:nlevs), tindex=i)
       call mpp_read(unit, field_error, profile_error, tindex=i)
       call mpp_read(unit, field_fix_depth, fix_depth, tindex=i) ! snz drop rate
       if ( var_id == temp_id ) then
          call mpp_read(unit, field_t_flag, t_flag(1:nlevs), tindex=i)
          call mpp_read(unit, field_flag, flag_t, tindex=i)
       else if ( var_id == salt_id ) then
          call mpp_read(unit, field_s_flag, s_flag(1:nlevs), tindex=i)
          call mpp_read(unit, field_flag, flag_s, tindex=i)
       end if
       call mpp_read(unit, field_link, rlink, tindex=i)

       inst_type = 20 ! snz change one line
!!$       inst_type = DRIFTER + ARGO

       data_is_local = .false.
       data_in_period = .false.

       if ( lon .lt. 0.0 ) lon = lon + 360.0
       if ( lon .gt. 360.0 ) lon = lon - 360.0
       if ( lon .lt. 80.0 ) lon = lon + 360.0

       if ( lat > ass_start_lat .and. lat < ass_end_lat ) data_is_local = .true.

       profile_time = get_cal_time(time, time_units, 'julian')
       if ( profile_time > time_start .and. profile_time < time_end ) data_in_period = .true.
       if ( (data_in_period .and. data_is_local) .and. (.NOT.localize_data) ) then ! localize

          if (isd_filt < 1 .and. ied_filt > ieg) then
             ! filter domain is a full x band
             if (lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  lat >= y_grid(1,jsd_flt0) .and. lat <= y_grid(ieg-1,jsd_flt0)) then
                prof_in_filt_domain = .true.
             end if
          else if (isd_filt >= 1 .and. ied_filt <= ieg) then
             ! Interior filter domain
            if (lon >= x_grid(isd_filt,jsd_flt0) .and. lon <= x_grid(ied_filt-1,jsd_flt0) .and.&
               & lat >= y_grid(isd_filt,jsd_flt0) .and. lat <=  y_grid(ied_filt-1,jed_flt0-1)) then
                prof_in_filt_domain = .true.
            end if
          else if (isd_filt < 1 .and. ied_filt <= ieg) then
             ! lhs filter domain
             isd_flt0 = isd_filt + ieg
             if ((lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ied_filt-1,jsd_flt0) .and.&
                  & lat >= y_grid(1,jsd_flt0) .and. lat <=  y_grid(ied_filt-1,jed_flt0-1)).or.&
                  & (lon >= x_grid(isd_flt0,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  & lat >= y_grid(isd_flt0,jsd_flt0) .and. lat <=  y_grid(ieg-1,jed_flt0-1))) then
                prof_in_filt_domain = .true.
             end if
          else if (isd_filt >= 1 .and. ied_filt > ieg) then
             ! rhs filter domain
             ied_flt0 = ied_filt - ieg
             if ( lon >= x_grid(isd_filt,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  & lat >= y_grid(isd_filt,jsd_flt0) .and. lat <=  y_grid(ieg-1,jed_flt0-1) ) then
                prof_in_filt_domain = .true.
             end if
             if (ied_flt0-1 > 1) then
                if ( lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ied_flt0-1,jsd_flt0) .and.&
                     & lat >= y_grid(1,jsd_flt0) .and. lat <=  y_grid(ied_flt0-1,jed_flt0-1) ) then
                   prof_in_filt_domain = .true.
                end if
             end if
          end if

          if ( var_id == temp_id .and. flag_t == 0.0 ) then
             num_profiles = num_profiles + 1
             no_temp = no_temp + 1
             no_prf = no_prf + 1
          end if
          if ( var_id == salt_id .and. flag_s == 0.0 ) then
             num_profiles = num_profiles + 1
             no_salt = no_salt + 1
             no_prf = no_prf + 1
          end if

          if ( num_profiles > max_profiles ) then
             call error_mesg('oda_core_mod::open_profile_dataset',&
                  & 'maximum number of profiles exceeded, increase max_profiles in oda_core_nml', fatal)
          end if

          num_levs = 0
          do k=1, max_levels
             flag(k) = .true.

             if ( depth(k) > 2000.0 ) depth(k) = missing_value  ! snz add for rdat-hybn

             if ( var_id == temp_id ) then
                if ( data(k) .eq. missing_value .or.&
                     & depth(k) .eq. missing_value .or. t_flag(k) .ne. 0.0 ) then
                   flag(k) = .false.
                else
                   num_levs = num_levs + 1
                end if
             else if ( var_id == salt_id ) then
                if ( data(k) .eq. missing_value .or.&
                     & depth(k) .eq. missing_value .or. s_flag(k) .ne. 0.0 ) then
                   flag(k) = .false.
                else
                   num_levs = num_levs+1
                end if
             end if
          end do

          ! large profile are stored externally in separate records
          ! read linked records and combine profile
          ii = i + 1
          nlinks = 0
          do while ( rlink > 0.0 )
             nlinks = nlinks + 1

             if ( nlinks > max_lnks ) then
                write (emsg_local, '("nlinks (",I6,") > MAX_LNKS (",I6,")")')&
                     & nlinks, max_lnks
                call error_mesg('oda_core_mod::open_profile_dataset',&
                     & trim(emsg_local)//' in file "'//trim(filename)//&
                     & '".  Increase parameter MAX_LNKS', fatal)
             end if

             depth_bfr(nlinks,:) = missing_value
             data_bfr(nlinks,:) = missing_value
             call mpp_read(unit, field_depth, depth_bfr(nlinks,1:nlevs), tindex=ii)
             call mpp_read(unit, field_data, data_bfr(nlinks,1:nlevs), tindex=ii)
             if ( var_id == temp_id ) then
                call mpp_read(unit, field_t_flag, t_flag_bfr(nlinks,1:nlevs), tindex=ii)
             else if ( var_id == salt_id ) then
                call mpp_read(unit, field_s_flag, s_flag_bfr(nlinks,1:nlevs), tindex=ii)
             end if
             call mpp_read(unit, field_link, rlink, tindex=ii)
             ii = ii + 1
          end do
          i = ii ! set record counter to start of next profile

          if ( nlinks > 0 ) then
             do nn=1, nlinks
                do k=1, max_levels
                   flag_bfr(nn,k) = .true.

                   if ( depth_bfr(nn,k) > 2000.0 ) depth_bfr(nn,k) = missing_value  ! snz add for rdat-hybn

                   if ( var_id == temp_id ) then
                      if ( data_bfr(nn,k) .eq. missing_value .or.&
                           & depth_bfr(nn,k) .eq. missing_value .or.&
                           & t_flag_bfr(nn,k) .ne. 0.0 ) then
                         flag_bfr(nn,k) = .false.
                      else
                         num_levs = num_levs+1
                      end if
                   else if (var_id == salt_id) then
                      if ( data_bfr(nn,k) .eq. missing_value .or.&
                           & depth_bfr(nn,k) .eq. missing_value .or.&
                           & s_flag_bfr(nn,k) .ne. 0.0 ) then
                         flag_bfr(nn,k) = .false.
                      else
                         num_levs = num_levs+1
                      end if
                   end if
                end do
             end do
          end if

          ! mh2 asks to change from [if (num_levs == 0) cycle]
          if ( num_levs == 0 ) then
             if ( i .gt. nstation ) cont = .false.
             cycle
          end if

          if ( num_profiles > 0 .and. prof_in_filt_domain ) then ! snz - 05 Nov 2012

             allocate(profiles(num_profiles)%depth(num_levs))
             allocate(profiles(num_profiles)%data(num_levs))
             allocate(profiles(num_profiles)%flag(num_levs))
             profiles(num_profiles)%variable = var_id
             if ( inst_type < 1 ) inst_type = unknown
             profiles(num_profiles)%inst_type = inst_type
             profiles(num_profiles)%levels = num_levs
             profiles(num_profiles)%lat = lat
             profiles(num_profiles)%lon = lon
!             allocate(profiles(num_profiles)%ms(num_levs))
!             allocate(profiles(num_profiles)%ms_inv(num_levs))
!             profiles(num_profiles)%ms(:) = 0.5
             kk = 1
             do k=1, max_levels
                if ( flag(k) ) then
                   if ( kk > profiles(num_profiles)%levels ) then
                      call error_mesg('oda_core_mod::open_profile_dataset',&
                           & 'Loop value "kk" is greater than profile levels', fatal)
                   end if
                   profiles(num_profiles)%depth(kk) = depth(k)
                   profiles(num_profiles)%data(kk) = data(k)
!                   profiles(num_profiles)%ms_inv(kk) = 1./profiles(num_profiles)%ms(kk)
                   kk = kk + 1
                end if
             end do

             do nn=1, nlinks
                do k=1, max_levels
                   if ( flag_bfr(nn,k) ) then
                      if ( kk > profiles(num_profiles)%levels ) then
                         call error_mesg('oda_core_mod::open_profile_dataset',&
                              & 'Loop value "kk" is greater than profile levels (bfr loop)', fatal)
                      end if
                      profiles(num_profiles)%depth(kk) = depth_bfr(nn,k)
                      profiles(num_profiles)%data(kk) = data_bfr(nn,k)
!                      profiles(num_profiles)%ms_inv(kk) = 1./profiles(num_profiles)%ms(kk)
                      kk = kk + 1
                   end if
                end do
             end do

             profiles(num_profiles)%time = profile_time

             ! calculate interpolation coefficients (make sure to account for grid offsets here!)
             if ( lat < 65.0 ) then ! regular grids
                ri0 = frac_index(lon, x_grid(:,1))
                rj0 = frac_index(lat, y_grid(90,:))
                i0 = floor(ri0)
                j0 = floor(rj0)
                if ( i0 > ieg .or. j0 > jeg ) then
                   write (unit=emsg_local, fmt='("i0 = ",I8,", j0 = ",I8)') mpp_pe(), i0, j0
                   call error_mesg('oda_core_mod::open_profile_dataset',&
                        & 'For regular grids, either i0 > ieg or j0 > jeg.  '//trim(emsg_local), fatal)
                end if
                if ( isd_filt >= 1 .and. ied_filt <= ieg ) then
                  if ( i0 < isd_filt .or. i0 > ied_filt .or. j0 < jsd_filt .or. j0 > jed_filt ) then
                     write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                          & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                     call error_mesg('oda_core_mod::open_profile_dataset',&
                          & 'i0,j0 out of bounds in prfs01. '//trim(emsg_local), fatal)
                  end if
                end if
                if ( isd_filt < 1 .and. i0 > ied_filt-1 .and. i0 < isd_filt + ieg ) then
                   write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                        & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                   call error_mesg('oda_core_mod::open_profile_dataset',&
                        & 'i0,j0 out of bounds in prfs02. '//trim(emsg_local), fatal)
                end if
                if ( ied_filt > ieg .and. i0 > ied_filt-ieg-1 .and. ied_filt < isd_filt ) then
                   write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                        & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                   call error_mesg('oda_core_mod::open_profile_dataset',&
                        & 'i0,j0 out of bounds in prfs03. '//trim(emsg_local), fatal)
                end if
                profiles(num_profiles)%i_index = ri0
                profiles(num_profiles)%j_index = rj0
             else ! tripolar grids
                lon_out(1,1) = (lon-360.0)*deg_to_rad
                lat_out(1,1) = lat*deg_to_rad
                call horiz_interp_bilinear_new (interp, (x_grid-360.0)*deg_to_rad, y_grid*deg_to_rad,&
                     & lon_out, lat_out, new_search=.true., no_crash_when_not_found=.true.)

                if ( interp%i_lon(1,1,1) == -999. ) bad_point = bad_point + 1
                if ( interp%wti(1,1,2) < 1.0 ) then
                   i0 = interp%i_lon(1,1,1)
                else
                   i0 = interp%i_lon(1,1,2)
                end if
                if ( interp%wtj(1,1,2) < 1.0 ) then
                   j0 = interp%j_lat(1,1,1)
                else
                   j0 = interp%j_lat(1,1,2)
                end if
                if ( i0 > ieg .or. j0 > jeg ) then
                   write (unit=emsg_local, fmt='("i0 = ",I6,", j0 = ",I6)') mpp_pe(), i0, j0
                   call error_mesg('oda_core_mod::open_profile_dataset',&
                        & 'For tripolar grids, either i0 > ieg or j0 > jeg', fatal)
                end if
                if( i0 < isd_filt .or. i0 > ied_filt .or. j0 < jsd_filt .or. j0 > jed_filt ) then
!!$                   print*,'prfs.pe,i0,j0= ',mpp_pe(), i0, j0,&
!!$                        & 'isd_filt,ied_filt,jsd_filt,jed_filt= ',isd_filt,ied_filt,jsd_filt,jed_filt
!!$                   print*,'pe,lon,lat=',mpp_pe(),lon,lat,'x_grid(i0+-1)',x_grid(i0-1:i0+1,j0),&
!!$                        & 'y_grid(i0,j0+-1)=',y_grid(i0,j0-1:j0+1)
!!$                   print*,'lono11,lato11=',x_grid(i0,j0),y_grid(i0,j0),'lono21,lato21=',x_grid(i0+1,j0),y_grid(i0+1,j0)
!!$                   print*,'lono12,lato12=',x_grid(i0,j0+1),y_grid(i0,j0+1),'lono22,lato22=',x_grid(i0+1,j0+1),y_grid(i0+1,j0+1)
!!$                   print*,'lonm11,latm11=',x_grid(isd_filt,jsd_filt),y_grid(isd_filt,jsd_filt),&
!!$                        & 'lonm21,latm21=',x_grid(ied_filt,jsd_filt),y_grid(ied_filt,jsd_filt)
!!$                   print*,'lonm12,latm12=',x_grid(isd_filt,jed_filt),y_grid(isd_filt,jed_filt),&
!!$                        & 'lonm22,latm22=',x_grid(ied_filt,jed_filt),y_grid(ied_filt,jed_filt)
!!$                   print*,'wti(1:2)=',Interp%wti(1,1,:),'wtj(1:2)=',Interp%wtj(1,1,:)

                   out_bound_point = out_bound_point + 1
                end if
                if ( interp%wti(1,1,2) < 1.0 ) then
                   profiles(num_profiles)%i_index =interp%i_lon(1,1,1) + interp%wti(1,1,2)
                else
                   profiles(num_profiles)%i_index =interp%i_lon(1,1,2)
                end if
                if (interp%wtj(1,1,2) < 1.0) then
                   profiles(num_profiles)%j_index =interp%j_lat(1,1,1) + interp%wtj(1,1,2)
                else
                   profiles(num_profiles)%j_index =interp%j_lat(1,1,2)
                end if
             end if ! grids

             profiles(num_profiles)%accepted = .true.

             if ( var_id == temp_id .and. flag_t /= 0.0 ) profiles(num_profiles)%accepted = .false.
             if ( var_id == salt_id .and. flag_s /= 0.0 ) profiles(num_profiles)%accepted = .false.

             if (i0 < 1 .or. j0 < 1) then
                profiles(num_profiles)%accepted = .false.
             end if
             if( i0 < isd_filt .or. i0 >= ied_filt .or. j0 < jsd_filt .or. j0 >= jed_filt ) then
                profiles(num_profiles)%accepted = .false.
             end if

             if ( profiles(num_profiles)%accepted ) then ! here
                if ( i0 /= ieg .and. j0 /= jeg ) then
                   if (grd%mask(i0,j0,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0,1) == 0.0 .or.&
                        & grd%mask(i0,j0+1,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0+1,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                else if ( i0 == ieg .and. j0 /= jeg ) then
                   if (grd%mask(i0,j0,1) == 0.0 .or.&
                        & grd%mask(1,j0,1) == 0.0 .or.&
                        & grd%mask(i0,j0+1,1) == 0.0 .or.&
                        & grd%mask(1,j0+1,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                else if ( i0 /= ieg .and. j0 == jeg ) then
                   if ( grd%mask(i0,j0,1) == 0.0 .or. grd%mask(i0+1,j0,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                else
                   if ( grd%mask(i0,j0,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                end if
             end if ! here

             if ( profiles(num_profiles)%accepted .and.&
                  & profiles(num_profiles)%inst_type == mooring+tao ) then
                if ( allocated(mask_tao) ) then
                   if ( mask_tao(i0,j0) < 1.0 ) then
                      profiles(num_profiles)%accepted = .false.
                      write (unit=stdout_unit,&
                           & fmt='("Rejecting tao mooring at (lat,lon) = (",F10.5,",",F10.5,") based on user-specified mask.")')&
                           & profiles(num_profiles)%lat,&
                           & profiles(num_profiles)%lon
                   end if
                end if
             end if

             if ( profiles(num_profiles)%accepted ) then ! accepted
                profiles(num_profiles)%flag(:) = .true.
                allocate(profiles(num_profiles)%k_index(profiles(num_profiles)%levels))
                do k=1, profiles(num_profiles)%levels
                   if (profiles(num_profiles)%depth(k) < grd%z(1)) then
                     profiles(num_profiles)%k_index(k) = 1.0
                   else
                     profiles(num_profiles)%k_index(k) = frac_index(profiles(num_profiles)%depth(k), (/0.,grd%z(:)/))! - 1 snz modify to v3.2 JAN3012
                   end if
                   if ( profiles(num_profiles)%k_index(k) < 1.0 ) then
                      if ( profiles(num_profiles)%depth(k) < 0.0 ) then
                         profiles(num_profiles)%k_index(k) = 0.0
                      else if ( profiles(num_profiles)%depth(k) > grd%z(size(grd%z,1)) ) then
                         profiles(num_profiles)%k_index(k) = real(nk)
                      end if
                   else
                      profiles(num_profiles)%k_index(k) = profiles(num_profiles)%k_index(k) - 1.0
                   end if
                   if ( profiles(num_profiles)%k_index(k) > real(nk) ) then
                      call error_mesg('oda_core_mod::open_profile_dataset', 'Profile k_index is greater than nk', fatal)
                   else if ( profiles(num_profiles)%k_index(k) < 0.0 ) then
                      call error_mesg('oda_core_mod::open_profile_dataset', 'Profile k_index is less than 0', fatal)
                   end if
                   k0 = floor(profiles(num_profiles)%k_index(k))

                   IF ( k0 >= 1 ) THEN ! snz add
                      if ( profiles(num_profiles)%flag(k) ) then ! flag
                         if ( i0 /= ieg .and. j0 /= jeg ) then
                            if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0) == 0.0 .or.&
                                 & grd%mask(i0+1,j0+1,k0) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else if ( i0 == ieg .and. j0 /= jeg ) then
                            if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                                 & grd%mask(1,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0) == 0.0 .or.&
                                 & grd%mask(1,j0+1,k0) == 0.0) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else if ( i0 /= ieg .and. j0 == jeg ) then
                            if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0) == 0.0) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else
                            if ( grd%mask(i0,j0,k0) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         end if

                         if ( i0 /= ieg .and. j0 /= jeg) then
                            if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0+1) == 0.0 .or.&
                                 & grd%mask(i0+1,j0+1,k0+1) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else if ( i0 == ieg .and. j0 /= jeg ) then
                            if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(1,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0+1) == 0.0 .or.&
                                 & grd%mask(1,j0+1,k0+1) == 0.0) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else if ( i0 /= ieg .and. j0 == jeg ) then
                            if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0+1) == 0.0) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else
                            if ( grd%mask(i0,j0,k0+1) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         end if

                         if ( abs(profiles(num_profiles)%data(k)) > 1.e4 &
                              & .or. abs(profiles(num_profiles)%depth(k)) > 1.e4 ) then
                            profiles(num_profiles)%flag(k) = .false.
                         end if
                      end if ! flag
                   end if ! snz add
                end do
             end if ! accepted
          endif ! 05 Nov 2012
       else ! localize
          i = i+1
       end if ! localize

       if ( var_id == temp_id .and. num_profiles > 0 ) call xbt_drop_rate_adjust(profiles(num_profiles))

       if ( i .gt. nstation ) cont = .false.
    end do

    a = nprof_in_filt_domain
    bad_point_g = bad_point
    call mpp_sum(a)
    call mpp_sum(bad_point_g)
    call mpp_sum(out_bound_point)

    if ( no_prf /= num_profiles ) then
       write(unit=stdout_unit, fmt='("PE: ",I6," no_prf = ",I8,", num_profiles = ",I8)') mpp_pe(), no_prf, num_profiles
    end if
    if ( var_id == temp_id ) then
       write(unit=stdout_unit, fmt='("A grand total of ",I8," temp prfs within global domain")') no_temp
       write(unit=stdout_unit, fmt='("The total of bad_point temp ",I8," within global domain")') bad_point_g
       write(unit=stdout_unit, fmt='("The total out_bound_point temp ",I8)') out_bound_point
    else if ( var_id == salt_id ) then
       write(unit=stdout_unit, fmt='("A grand total of ",I8," salt prfs within global domain")') no_salt
       write(unit=stdout_unit, fmt='("The total of bad_point salt",I8," within global domain")') bad_point_g
       write(unit=stdout_unit, fmt='("A grand total of ",I8," prfs within global domain")') no_prf
       write(unit=stdout_unit, fmt='("A grand total of ",I8," prfs within current PEs computer domain")') a
       write(unit=stdout_unit, fmt='("The total out_bound_point salt ",I8)') out_bound_point
    end if

    call mpp_sync_self()
    call mpp_close(unit)
    deallocate(axes)
    deallocate(fields)

!    call mpp_print_memuse_stats('open_profile_dataset End')

  end subroutine open_profile_dataset

  subroutine open_profile_dataset_gtspp()
    return
  end subroutine open_profile_dataset_gtspp

  subroutine open_profile_dataset_argo(filename, time_start, time_end, obs_variable, localize)
    character(len=*), intent(in) :: filename
    type(time_type), intent(in) :: time_start, time_end
    integer, intent(in) :: obs_variable
    logical, intent(in), optional :: localize

    integer, parameter :: MAX_LEVELS = 1000
    integer, parameter :: MAX_LNKS = 500

    real :: lon, lat, time, rlink, prf_type
    real :: ri0, rj0
    real, dimension(MAX_LEVELS) :: depth, data
    real, dimension(MAX_LNKS, MAX_LEVELS) :: data_bfr, depth_bfr

    integer :: unit, ndim, nvar, natt, nstation
    integer :: stdout_unit
    integer :: inst_type, var_id
    integer :: num_levs, k, kk, i, i0, j0, k0, nlevs, a, nn, ii, nlinks
    integer :: nprof_in_filt_domain, out_bound_point

    character(len=32) :: fldname, axisname, anal_fldname, time_units
    character(len=128) :: emsg_local

    logical :: data_is_local, localize_data, cont
    logical :: data_in_period
    logical :: prof_in_filt_domain
    logical, dimension(MAX_LEVELS) :: flag
    logical, dimension(MAX_LNKS, MAX_LEVELS) :: flag_bfr

    type(time_type) :: profile_time
    type(axistype), pointer :: depth_axis, station_axis
    type(axistype), allocatable, dimension(:), target :: axes
    type(fieldtype), allocatable, dimension(:), target :: fields
    type(fieldtype), pointer :: field_lon, field_lat, field_flag, field_time
    type(fieldtype), pointer :: field_depth, field_data, field_link, field_var_type

    ! NOTE: fields are restricted to be in separate files

    stdout_unit = stdout()

    if ( PRESENT(localize) ) then
       localize_data = localize
    else
       localize_data = .true.
    end if

    nprof_in_filt_domain = 0

    anal_fldname = 'none'
    var_id=-1
    if ( obs_variable == temp_id ) then
       anal_fldname = 'temp'
       var_id = temp_id
    else if ( obs_variable == salt_id ) then
       anal_fldname = 'salt'
       var_id = salt_id
    end if

!    call mpp_print_memuse_stats('open_profile_dataset_argo Start')

    call mpp_open(unit, filename, form=mpp_netcdf, fileset=mpp_single, threading=mpp_multi, action=mpp_rdonly)
    call mpp_get_info(unit, ndim, nvar, natt, nstation)

    write (unit=stdout_unit, fmt='("Opened profile dataset: ",A)') trim(filename)

    ! get axis information

    allocate(axes(ndim))
    call mpp_get_axes(unit, axes)
    do i=1, ndim
       call mpp_get_atts(axes(i), name=axisname)
       select case (trim(axisname))
       case ('depth_index')
          depth_axis => axes(i)
       case ('station_index')
          station_axis => axes(i)
       end select
    end do

    ! get field information
    allocate(fields(nvar))
    call mpp_get_fields(unit, fields)
    do i=1, nvar
       call mpp_get_atts(fields(i), name=fldname)
       if( var_id .eq. temp_id ) then
          select case (trim(fldname))
          case ('longitude')
             field_lon => fields(i)
          case ('latitude')
             field_lat => fields(i)
          case ('dens_flag')
             field_flag => fields(i)
          case ('time')
             field_time => fields(i)
          case ('temp')
             field_data => fields(i)
          case ('depth')
             field_depth => fields(i)
          case ('link')
             field_link => fields(i)
          case ('var_type')
             field_var_type => fields(i)
          end select
       else if( var_id .eq. salt_id ) then
          select case (trim(fldname))
          case ('longitude')
             field_lon => fields(i)
          case ('latitude')
             field_lat => fields(i)
          case ('dens_flag')
             field_flag => fields(i)
          case ('time')
             field_time => fields(i)
          case ('salt')
             field_data => fields(i)
          case ('depth')
             field_depth => fields(i)
          case ('link')
             field_link => fields(i)
          case ('var_type')
             field_var_type => fields(i)
          end select
       end if
    end do

    call mpp_get_atts(depth_axis, len=nlevs)

    if ( nlevs > max_levels ) then
       call error_mesg('oda_core_mod::open_profile_dataset_argo', 'increase parameter MAX_LEVELS', fatal)
    else if (nlevs < 1) then
       call error_mesg('oda_core_mod::open_profile_dataset_argo', 'nlevs less than 1.', fatal)
    end if

    if ( .NOT.ASSOCIATED(field_data) ) then
       call error_mesg('oda_core_mod::open_profile_dataset_argo',&
            & 'profile dataset not used because data not needed for Analysis', note)
       return
    end if

    write (unit=stdout_unit, fmt='("There are ",I8," records in this dataset")') nstation
    write (unit=stdout_unit, fmt='("Searching for profiles . . .")')

    call mpp_get_atts(field_time, units=time_units)

    out_bound_point = 0

    i=1
    cont=.true.

    do while (cont)
       prof_in_filt_domain = .false.
       depth = missing_value  ! snz add
       data = missing_value   ! snz add

       call mpp_read(unit, field_lon, lon, tindex=i)
       call mpp_read(unit, field_lat, lat, tindex=i)
       call mpp_read(unit, field_time, time, tindex=i)
       call mpp_read(unit, field_depth, depth(1:nlevs), tindex=i)
       call mpp_read(unit, field_data, data(1:nlevs), tindex=i)
       call mpp_read(unit, field_var_type, prf_type, tindex=i)
       call mpp_read(unit, field_link, rlink, tindex=i)

!!$       inst_type = DRIFTER + ARGO
       inst_type = 20 ! snz change one line
       data_is_local = .false.
       data_in_period = .false.

       if ( lon .lt. 0.0 ) lon = lon + 360.0
       if ( lon .gt. 360.0 ) lon = lon - 360.0
       if ( lon .lt. 80.0 ) lon = lon + 360.0

       if ( lat > ass_start_lat .and. lat < ass_end_lat ) data_is_local = .true.

       profile_time = get_cal_time(time, time_units, 'NOLEAP')
       if ( profile_time > time_start .and. profile_time < time_end ) data_in_period = .true.
       if ( (data_in_period .and. data_is_local) .and. (.NOT.localize_data) ) then

          if (isd_filt < 1 .and. ied_filt > ieg) then
             ! filter domain is a full x band
             if (lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  lat >= y_grid(1,jsd_flt0) .and. lat <= y_grid(ieg-1,jsd_flt0)) then
                prof_in_filt_domain = .true.
             end if
          else if (isd_filt >= 1 .and. ied_filt <= ieg) then
             ! Interior filter domain
            if (lon >= x_grid(isd_filt,jsd_flt0) .and. lon <= x_grid(ied_filt-1,jsd_flt0) .and.&
               & lat >= y_grid(isd_filt,jsd_flt0) .and. lat <=  y_grid(ied_filt-1,jed_flt0-1)) then
                prof_in_filt_domain = .true.
            end if
          else if (isd_filt < 1 .and. ied_filt <= ieg) then
             ! lhs filter domain
             isd_flt0 = isd_filt + ieg
             if ((lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ied_filt-1,jsd_flt0) .and.&
                  & lat >= y_grid(1,jsd_flt0) .and. lat <=  y_grid(ied_filt-1,jed_flt0-1)).or.&
                  & (lon >= x_grid(isd_flt0,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  & lat >= y_grid(isd_flt0,jsd_flt0) .and. lat <=  y_grid(ieg-1,jed_flt0-1))) then
                prof_in_filt_domain = .true.
             end if
          else if (isd_filt >= 1 .and. ied_filt > ieg) then
             ! rhs filter domain
             ied_flt0 = ied_filt - ieg
             if ( lon >= x_grid(isd_filt,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  & lat >= y_grid(isd_filt,jsd_flt0) .and. lat <=  y_grid(ieg-1,jed_flt0-1) ) then
                prof_in_filt_domain = .true.
             end if
             if (ied_flt0-1 > 1) then
                if ( lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ied_flt0-1,jsd_flt0) .and.&
                     & lat >= y_grid(1,jsd_flt0) .and. lat <=  y_grid(ied_flt0-1,jed_flt0-1) ) then
                   prof_in_filt_domain = .true.
                end if
             end if
          end if

          if ( var_id == temp_id ) then
             num_profiles = num_profiles + 1
             no_temp = no_temp + 1
             no_prf = no_prf + 1
          else if ( var_id == salt_id .and. prf_type == 2.0 ) then
             num_profiles =num_profiles + 1
             no_salt = no_salt + 1
             no_prf = no_prf + 1
          end if

          if ( num_profiles > max_profiles ) then
             call error_mesg('oda_core_mod::open_profile_dataset_argo',&
                  & 'maximum number of profiles exceeded, increase max_profiles in oda_core_nml', fatal)
          end if

          num_levs = 0
          do k=1, max_levels
             flag(k) = .true.
             if ( depth(k) > 2000.0 ) depth(k) = missing_value  ! snz add for rdat-hybn
             if ( var_id == temp_id ) then
                if ( data(k) .eq. missing_value .or. depth(k) .eq. missing_value ) then
                   flag(k) = .false.
                else
                   num_levs = num_levs+1
                end if
             else if ( var_id == salt_id ) then
                if ( data(k) .eq. missing_value .or. depth(k) .eq. missing_value ) then
                   flag(k) = .false.
                else
                   num_levs = num_levs+1
                end if
             end if
          end do

          ! large profile are stored externally in separate records
          ! read linked records and combine profile
          ii=i+1
          nlinks = 0
          do while ( rlink > 0.0 )
             nlinks = nlinks + 1

             if ( nlinks > max_lnks ) then
                write (emsg_local, '("nlinks (",I6,") > MAX_LNKS (",I6,")")')&
                     & nlinks, max_lnks
                call error_mesg('oda_core_mod::open_profile_dataset_argo',&
                     & trim(emsg_local)//' in file "'//trim(filename)//&
                     & '".  Increase parameter MAX_LNKS', fatal)
             end if

             depth_bfr(nlinks,:) = missing_value
             data_bfr(nlinks,:) = missing_value
             call mpp_read(unit,field_depth,depth_bfr(nlinks,1:nlevs),tindex=ii)
             call mpp_read(unit,field_data,data_bfr(nlinks,1:nlevs),tindex=ii)
             call mpp_read(unit,field_link,rlink,tindex=ii)
             ii=ii+1
          end do
          i=ii ! set record counter to start of next profile

          if ( nlinks > 0 ) then
             do nn=1, nlinks
                do k=1, max_levels
                   flag_bfr(nn,k) = .true.

                   if ( depth_bfr(nn,k) > 2000.0 ) depth_bfr(nn,k) = missing_value  ! snz add for rdat-hybn

                   if ( var_id == temp_id ) then
                      if ( data_bfr(nn,k) .eq. missing_value .or. depth_bfr(nn,k) .eq. missing_value ) then
                         flag_bfr(nn,k) = .false.
                      else
                         num_levs = num_levs+1
                      end if
                   else if ( var_id == salt_id ) then
                      if ( data_bfr(nn,k) .eq. missing_value .or. depth_bfr(nn,k) .eq. missing_value ) then
                         flag_bfr(nn,k) = .false.
                      else
                         num_levs = num_levs+1
                      end if
                   end if
                end do
             end do
          end if

          ! mh2 asks to change from [if (num_levs == 0) cycle]
          if ( num_levs == 0 ) then
             if ( i .gt. nstation ) cont = .false.
             cycle
          end if

          if (nprof_in_filt_domain > 0 .and. prof_in_filt_domain) then ! snz 05 Nov 2012

          allocate(profiles(num_profiles)%depth(num_levs))
          allocate(profiles(num_profiles)%data(num_levs))
          allocate(profiles(num_profiles)%flag(num_levs))
          profiles(num_profiles)%variable = var_id
          if ( inst_type < 1 ) inst_type = unknown
          profiles(num_profiles)%inst_type = inst_type
          profiles(num_profiles)%levels = num_levs
          profiles(num_profiles)%lat = lat
          profiles(num_profiles)%lon = lon
!          allocate(profiles(num_profiles)%ms(num_levs))
!          allocate(profiles(num_profiles)%ms_inv(num_levs))
!          profiles(num_profiles)%ms(:) = 0.5

          kk= 1
          do k=1, max_levels
             if ( flag(k) ) then
                if ( kk > profiles(num_profiles)%levels ) then
                   call error_mesg('oda_core_mod::open_profile_dataset_argo',&
                        & 'Loop variable "kk" is greater than profile levels', fatal)
                end if
                profiles(num_profiles)%depth(kk) = depth(k)
                profiles(num_profiles)%data(kk) = data(k)
!                profiles(num_profiles)%ms_inv(kk) = 1./profiles(num_profiles)%ms(kk)
                kk = kk + 1
             end if
          end do

          do nn=1, nlinks
             do k=1, max_levels
                if ( flag_bfr(nn,k) ) then
                   if ( kk > profiles(num_profiles)%levels ) then
                      call error_mesg('oda_core_mod::open_profile_dataset_argo',&
                           & 'Loop variable "kk" is greater than profile levels (bfr loop)', fatal)
                   end if
                   profiles(num_profiles)%depth(kk) = depth_bfr(nn,k)
                   profiles(num_profiles)%data(kk) = data_bfr(nn,k)
!                   profiles(num_profiles)%ms_inv(kk) = 1./profiles(num_profiles)%ms(kk)
                   kk = kk + 1
                end if
             end do
          end do

          profiles(num_profiles)%time = profile_time

! snz uses the following to test excluding the coast area salt profiles
!          if (profiles(num_profiles)%variable == SALT_ID .and. &
!              profiles(num_profiles)%depth(num_levs) < 900.0) profiles(num_profiles)%accepted = .false.

          ! calculate interpolation coefficients (make sure to account for grid offsets here!)
          ! note that this only works for lat/lon grids
          if ( lat < 65.0 ) then ! regular grids
             ri0 = frac_index(lon, x_grid(:,1))
             rj0 = frac_index(lat, y_grid(90,:))
             i0 = floor(ri0)
             j0 = floor(rj0)
             if ( i0 > ieg .or. j0 > jeg ) then
                write (unit=emsg_local, fmt='("i0 = ",I6,", j0 = ",I6)') mpp_pe(), i0, j0
                call error_mesg('oda_core_mod::open_profile_dataset_argo',&
                     & 'For regular grids, either i0 > ieg or j0 > jeg.  '//trim(emsg_local), fatal)
             end if
             if ( isd_filt >= 1 .and. ied_filt <= ieg ) then
               if ( i0 < isd_filt .or. i0 > ied_filt .or. j0 < jsd_filt .or. j0 > jed_filt ) then
                  write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                       & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                  call error_mesg('oda_core_mod::open_profile_dataset',&
                       & 'i0,j0 out of bounds in argo01. '//trim(emsg_local), fatal)
               end if
             end if
             if ( isd_filt < 1 .and. i0 > ied_filt-1 .and. i0 < isd_filt + ieg ) then
                write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                     & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                call error_mesg('oda_core_mod::open_profile_dataset',&
                     & 'i0,j0 out of bounds in argo02. '//trim(emsg_local), fatal)
             end if
             if ( ied_filt > ieg .and. i0 > ied_filt-ieg-1 .and. ied_filt < isd_filt ) then
                write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                     & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                call error_mesg('oda_core_mod::open_profile_dataset',&
                     & 'i0,j0 out of bounds in argo03. '//trim(emsg_local), fatal)
             end if
             profiles(num_profiles)%i_index = ri0
             profiles(num_profiles)%j_index = rj0
          else ! tripolar grids
             lon_out(1,1) = lon*deg_to_rad
             lat_out(1,1) = lat*deg_to_rad
             call horiz_interp_bilinear_new (interp, x_grid*deg_to_rad, y_grid*deg_to_rad, lon_out, lat_out)
             if(interp%wti(1,1,2) < 1.0) then
                i0 = interp%i_lon(1,1,1)
             else
                i0 = interp%i_lon(1,1,2)
             end if
             if ( interp%wtj(1,1,2) < 1.0 ) then
                j0 = interp%j_lat(1,1,1)
             else
                j0 = interp%j_lat(1,1,2)
             end if
             if ( i0 > ieg .or. j0 > jeg ) then
                write (unit=emsg_local, fmt='("i0 = ",I6,", j0 = ",I6)') mpp_pe(), i0, j0
                call error_mesg('oda_core_mod::open_profile_dataset_argo',&
                     & 'For tirpolar grids, either i0 > ieg or j0 > jeg.  '//trim(emsg_local), fatal)
             end if
             if ( i0 < isd_filt .or. i0 > ied_filt .or. j0 < jsd_filt .or. j0 > jed_filt ) then
!!$                print*,'argo.pe,i0,j0= ',mpp_pe(), i0, j0,&
!!$                     & 'isd_filt,ied_filt,jsd_filt,jed_filt= ',isd_filt,ied_filt,jsd_filt,jed_filt
!!$                print*,'pe,lon,lat=',mpp_pe(),lon,lat,'x_grid(i0+-1)',x_grid(i0-1:i0+1,j0),&
!!$                     & 'y_grid(i0,j0+-1)=',y_grid(i0,j0-1:j0+1)
!!$                print*,'lono11,lato11=',x_grid(i0,j0),y_grid(i0,j0),'lono21,lato21=',x_grid(i0+1,j0),y_grid(i0+1,j0)
!!$                print*,'lono12,lato12=',x_grid(i0,j0+1),y_grid(i0,j0+1),'lono22,lato22=',x_grid(i0+1,j0+1),y_grid(i0+1,j0+1)
!!$                print*,'lonm11,latm11=',x_grid(isd_filt,jsd_filt),y_grid(isd_filt,jsd_filt),&
!!$                     & 'lonm21,latm21=',x_grid(ied_filt,jsd_filt),y_grid(ied_filt,jsd_filt)
!!$                print*,'lonm12,latm12=',x_grid(isd_filt,jed_filt),y_grid(isd_filt,jed_filt),&
!!$                     & 'lonm22,latm22=',x_grid(ied_filt,jed_filt),y_grid(ied_filt,jed_filt)
!!$                print*,'wti(1:2)=',Interp%wti(1,1,:),'wtj(1:2)=',Interp%wtj(1,1,:)

                out_bound_point = out_bound_point + 1
             end if
             if ( interp%wti(1,1,2) < 1.0 ) then
                profiles(num_profiles)%i_index =interp%i_lon(1,1,1) + interp%wti(1,1,2)
             else
                profiles(num_profiles)%i_index =interp%i_lon(1,1,2)
             end if
             if ( interp%wtj(1,1,2) < 1.0 ) then
                profiles(num_profiles)%j_index =interp%j_lat(1,1,1) + interp%wtj(1,1,2)
             else
                profiles(num_profiles)%j_index =interp%j_lat(1,1,2)
             end if
          end if ! grids

           profiles(num_profiles)%accepted = .true.
           if ( i0 < 1 .or. j0 < 1 ) then
              profiles(num_profiles)%accepted = .false.
           end if
           if ( i0 < isd_filt .or. i0 >= ied_filt .or. j0 < jsd_filt .or. j0 >= jed_filt ) then
              profiles(num_profiles)%accepted = .false.
           end if

           if ( profiles(num_profiles)%accepted ) then ! here
              if ( i0 /= ieg .and. j0 /= jeg ) then
                 if ( grd%mask(i0,j0,1) == 0.0 .or.&
                      & grd%mask(i0+1,j0,1) == 0.0 .or.&
                      & grd%mask(i0,j0+1,1) == 0.0 .or.&
                      & grd%mask(i0+1,j0+1,1) == 0.0) then
                    profiles(num_profiles)%accepted = .false.
                 end if
              else if ( i0 == ieg .and. j0 /= jeg ) then
                 if (grd%mask(i0,j0,1) == 0.0 .or.&
                      & grd%mask(1,j0,1) == 0.0 .or.&
                      & grd%mask(i0,j0+1,1) == 0.0 .or.&
                      & grd%mask(1,j0+1,1) == 0.0) then
                    profiles(num_profiles)%accepted = .false.
                 end if
              else if ( i0 /= ieg .and. j0 == jeg ) then
                 if (grd%mask(i0,j0,1) == 0.0 .or.&
                      & grd%mask(i0+1,j0,1) == 0.0) then
                    profiles(num_profiles)%accepted = .false.
                 end if
              else
                 if ( grd%mask(i0,j0,1) == 0.0 ) then
                    profiles(num_profiles)%accepted = .false.
                 end if
              end if
           end if ! here

           if ( profiles(num_profiles)%accepted .and. profiles(num_profiles)%inst_type == mooring+tao) then
              if ( allocated(mask_tao) ) then
                 if ( mask_tao(i0,j0) < 1.0 ) then
                    profiles(num_profiles)%accepted = .false.
                    write (unit=stdout_unit,&
                         & fmt='("Rejecting tao mooring at (lat,lon) = (",F10.5,",",F10.5,") based on user-specified mask.")')&
                         & profiles(num_profiles)%lat,&
                         & profiles(num_profiles)%lon
                 end if
              end if
           end if

           if ( profiles(num_profiles)%accepted ) then
              profiles(num_profiles)%flag(:) = .true.
              allocate(profiles(num_profiles)%k_index(profiles(num_profiles)%levels))
              do k=1, profiles(num_profiles)%levels
                 if (profiles(num_profiles)%depth(k) < grd%z(1)) then
                   profiles(num_profiles)%k_index(k) = 1.0
                 else
                   profiles(num_profiles)%k_index(k) = frac_index(profiles(num_profiles)%depth(k), (/0.,grd%z(:)/))! - 1 snz modify to v3.2 JAN3012
                 end if
                 if ( profiles(num_profiles)%k_index(k) < 1 ) then
                    if ( profiles(num_profiles)%depth(k) < 0 ) then
                       profiles(num_profiles)%k_index(k) = 0
                    else if ( profiles(num_profiles)%depth(k) > grd%z(size(grd%z,1)) ) then
                       profiles(num_profiles)%k_index(k) = nk
                    end if
                 else
                    profiles(num_profiles)%k_index(k) = profiles(num_profiles)%k_index(k) - 1
                 end if
                 if ( profiles(num_profiles)%k_index(k) > nk ) then
                    call error_mesg('oda_core_mod::open_profile_dataset_argo', 'Profile k_index is greater than nk', fatal)
                 else if ( profiles(num_profiles)%k_index(k) < 0 ) then
                    call error_mesg('oda_core_mod::open_profile_dataset_argo', 'Profile k_index is less than 0', fatal)
                 end if
                 k0 = floor(profiles(num_profiles)%k_index(k))

                 if ( k0 >= 1 ) then ! snz add
                    if ( profiles(num_profiles)%flag(k) ) then ! flag
                       if ( i0 /= ieg .and. j0 /= jeg ) then
                          if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                               & grd%mask(i0+1,j0,k0) == 0.0 .or.&
                               & grd%mask(i0,j0+1,k0) == 0.0 .or.&
                               & grd%mask(i0+1,j0+1,k0) == 0.0 ) then
                             profiles(num_profiles)%flag(k) = .false.
                          end if
                       else if ( i0 == ieg .and. j0 /= jeg ) then
                          if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                               & grd%mask(1,j0,k0) == 0.0 .or.&
                               & grd%mask(i0,j0+1,k0) == 0.0 .or.&
                               & grd%mask(1,j0+1,k0) == 0.0 ) then
                             profiles(num_profiles)%flag(k) = .false.
                          end if
                       else if ( i0 /= ieg .and. j0 == jeg ) then
                          if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                               & grd%mask(i0+1,j0,k0) == 0.0 ) then
                             profiles(num_profiles)%flag(k) = .false.
                          end if
                       else
                          if ( grd%mask(i0,j0,k0) == 0.0 ) then
                             profiles(num_profiles)%flag(k) = .false.
                          end if
                       end if

                       if ( i0 /= ieg .and. j0 /= jeg ) then
                          if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                               & grd%mask(i0+1,j0,k0+1) == 0.0 .or.&
                               & grd%mask(i0,j0+1,k0+1) == 0.0 .or.&
                               & grd%mask(i0+1,j0+1,k0+1) == 0.0 ) then
                             profiles(num_profiles)%flag(k) = .false.
                          end if
                       else if ( i0 == ieg .and. j0 /= jeg ) then
                          if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                               & grd%mask(1,j0,k0+1) == 0.0 .or.&
                               & grd%mask(i0,j0+1,k0+1) == 0.0 .or.&
                               & grd%mask(1,j0+1,k0+1) == 0.0) then
                             profiles(num_profiles)%flag(k) = .false.
                          end if
                       else if ( i0 /= ieg .and. j0 == jeg ) then
                          if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                               & grd%mask(i0+1,j0,k0+1) == 0.0 ) then
                             profiles(num_profiles)%flag(k) = .false.
                          end if
                       else
                          if ( grd%mask(i0,j0,k0+1) == 0.0 ) then
                             profiles(num_profiles)%flag(k) = .false.
                          end if
                       end if

                       if ( abs(profiles(num_profiles)%data(k)) > 1.e4 &
                            & .or. abs(profiles(num_profiles)%depth(k)) > 1.e4 ) then
                          profiles(num_profiles)%flag(k) = .false.
                       end if
                    end if ! flag
                 end if ! snz add
              end do
           end if  ! accepted

           end if ! 05 Nov 2012

        else ! localize
           i = i+1
        end if ! localize

        if ( i .gt. nstation ) cont = .false.
     end do

     a = nprof_in_filt_domain
     call mpp_sum(a)
     call mpp_sum(out_bound_point)

    if ( no_prf /= num_profiles ) then
       write(unit=stdout_unit, fmt='("PE: ",I6," no_prf = ",I8,", num_profiles = ",I8)') mpp_pe(), no_prf, num_profiles
    end if
    if ( var_id == temp_id ) then
       write(unit=stdout_unit, fmt='("A grand total of ",I8," argo temp prfs within global domain")') no_temp
       write(unit=stdout_unit, fmt='("A total out of bound points",I8," argo temp within global domain")') out_bound_point
    else if ( var_id == salt_id ) then
       write(unit=stdout_unit, fmt='("A grand total of ",I8," argo salt prfs within global domain")') no_salt
       write(unit=stdout_unit, fmt='("A grand total of ",I8," argo prfs within global domain")') no_prf
       write(unit=stdout_unit, fmt='("A grand total of ",I8," argo prfs within current PEs computer domain")') a
       write(unit=stdout_unit, fmt='("A total out of bound points",I8," argo salt within global domain")') out_bound_point
    end if

    call mpp_sync_self()
    call mpp_close(unit)
    deallocate(axes)
    deallocate(fields)

!    call mpp_print_memuse_stats('open_profile_dataset_argo End')

  end subroutine open_profile_dataset_argo

  ! get profiles and sfc
  ! obs relevant to current analysis interval
  subroutine get_obs(model_time, Prof, nprof)
    type(time_type), intent(in) :: model_time
    type(ocean_profile_type), dimension(:), intent(inout) :: prof
    integer, intent(inout) :: nprof

    integer :: i, k, kk, k_interval
    integer :: yr, mon, day, hr, min, sec
    integer :: stdout_unit

    type(time_type) :: tdiff

    nprof = 0
    stdout_unit = stdout()

    write (unit=stdout_unit, fmt='("Gathering profiles for current analysis time")')
    call get_date(model_time, yr, mon, day, hr, min, sec)
    write (unit=stdout_unit, fmt='("Current YYYY/MM/DD = ",I4,"/",I2,"/",I2)') yr, mon, day

    do i=1, no_prf
       if ( profiles(i)%time <= model_time ) then
          tdiff = model_time - profiles(i)%time
       else
          tdiff = profiles(i)%time - model_time
       end if

       ! no tdiff criteria for monthly mean data like
       ! but tdiff criteria has to be set for daily data
       if ( tdiff <= time_window(profiles(i)%inst_type) .and. profiles(i)%accepted ) then
          ! for single profile test

          nprof = nprof + 1
          if ( nprof > size(prof,1) ) then
             call error_mesg('oda_core_mod::get_obs',&
                  & 'Passed in array "Prof" is smaller than number of profiles, increase size of Prof before call.',&
                  & fatal)
          end if
          call copy_obs(profiles(i:i), prof(nprof:nprof))

          prof(nprof)%tdiff = tdiff

          ! snz add the following few lines for increasing deep water data
          if ( prof(nprof)%levels > max_prflvs ) then
             k_interval = (prof(nprof)%levels-max_prflvs+50)/50 + 1
             kk = max_prflvs - 50
             do k=max_prflvs-50+1, prof(nprof)%levels, k_interval
                kk = kk + 1
                prof(nprof)%depth(kk) = prof(nprof)%depth(k)
                prof(nprof)%k_index(kk) = prof(nprof)%k_index(k)
                prof(nprof)%data(kk) = prof(nprof)%data(k)
                prof(nprof)%flag(kk) = prof(nprof)%flag(k)
             end do
             prof(nprof)%levels = kk
          end if
          ! snz end the adding lines
       end if
    end do

    write (unit=stdout_unit,&
         & fmt='("A total of ",I8," profiles are being used for the current analysis step.")') nprof

    return
  end subroutine get_obs

  subroutine oda_core_init(Domain, Grid, time_s, time_e, filt_domain, localize)
    type(domain2d), intent(inout) :: domain
    type(grid_type), target, intent(in) :: grid
    logical, intent(in), optional :: localize
    type(time_type), intent(in) :: time_s, time_e
    type(domain2d), intent(in) :: filt_domain

    integer :: ioun, ierr, io_status
    integer :: stdlog_unit

    stdlog_unit = stdlog()

    ! Read in the namelist file
#ifdef INTERNAL_FILE_NML
    read (input_nml_file, oda_core_nml, iostat=io_status)
#else
    ioun = open_namelist_file()
    read(ioun, nml=oda_core_nml, iostat=io_status)
    ierr = check_nml_error(io_status, 'oda_core_nml')
    call close_file(ioun)
#endif

    write(stdlog_unit, nml=oda_core_nml)

    grd => grid

    call mpp_get_compute_domain(domain, isc, iec, jsc, jec)
    call mpp_get_data_domain(domain, isd, ied, jsd, jed)
    call mpp_get_global_domain(domain, isg, ieg, jsg, jeg)
    call mpp_get_data_domain(filt_domain, isd_filt, ied_filt, jsd_filt, jed_filt)

    jsd_flt0 = jsd_filt
    jed_flt0 = jed_filt
    if (jsd_filt < 1) jsd_flt0 = 1
    if (jed_filt > jeg) jed_flt0 = jeg

    nk = size(grid%z)

    call init_observations(time_s, time_e, filt_domain, localize)
  end subroutine oda_core_init

  subroutine copy_obs(obs_in, obs_out)
    type(ocean_profile_type), dimension(:), intent(in) :: obs_in
    type(ocean_profile_type), dimension(:), intent(inout) :: obs_out

    integer :: n

    if ( size(obs_in) .ne. size(obs_out) ) then
       call error_mesg('oda_core_mod::copy_obs', 'Size of in and out obs variables are not equal.', fatal)
    end if

    do n=1, size(obs_in)
       obs_out(n)%variable = obs_in(n)%variable
       obs_out(n)%inst_type = obs_in(n)%inst_type
       obs_out(n)%levels = obs_in(n)%levels
       obs_out(n)%lon = obs_in(n)%lon
       obs_out(n)%lat = obs_in(n)%lat
       obs_out(n)%accepted = obs_in(n)%accepted
       if ( associated(obs_out(n)%depth) ) then
          deallocate(obs_out(n)%depth)
          nullify(obs_out(n)%depth)
       end if
       allocate(obs_out(n)%depth(obs_in(n)%levels))
       obs_out(n)%depth(:) = obs_in(n)%depth(:)
       if ( associated(obs_out(n)%data) ) then
          deallocate(obs_out(n)%data)
          nullify(obs_out(n)%data)
       end if
       allocate(obs_out(n)%data(obs_in(n)%levels))
       obs_out(n)%data(:) = obs_in(n)%data(:)
       if ( associated(obs_out(n)%flag) ) then
          deallocate(obs_out(n)%flag)
          nullify(obs_out(n)%flag)
       end if
       allocate(obs_out(n)%flag(obs_in(n)%levels))
       obs_out(n)%flag(:) = obs_in(n)%flag(:)
       obs_out(n)%time = obs_in(n)%time
       obs_out(n)%yyyy = obs_in(n)%yyyy
       obs_out(n)%mmdd = obs_in(n)%mmdd
       obs_out(n)%i_index = obs_in(n)%i_index
       obs_out(n)%j_index = obs_in(n)%j_index
       if ( associated(obs_out(n)%k_index) ) then
          deallocate(obs_out(n)%k_index)
          nullify(obs_out(n)%k_index)
       end if
       allocate(obs_out(n)%k_index(obs_in(n)%levels))
       obs_out(n)%k_index = obs_in(n)%k_index

!       if ( associated(Obs_out(n)%ms) ) then
!          deallocate(Obs_out(n)%ms)
!          nullify(Obs_out(n)%ms)
!       end if
!       allocate(Obs_out(n)%ms(Obs_in(n)%levels))
!       Obs_out(n)%ms = Obs_in(n)%ms
!       if ( associated(Obs_out(n)%ms_inv) ) then
!          deallocate(Obs_out(n)%ms_inv)
!          nullify(Obs_out(n)%ms_inv)
!       end if
!       allocate(Obs_out(n)%ms_inv(Obs_in(n)%levels))
!       Obs_out(n)%ms_inv = 1./Obs_in(n)%ms

       obs_out(n)%tdiff = obs_in(n)%tdiff
       if ( associated(obs_out(n)%Forward_model%wgt) ) then
          deallocate(obs_out(n)%Forward_model%wgt)
          nullify(obs_out(n)%Forward_model%wgt)
       end if
    end do
  end subroutine copy_obs

  subroutine open_profile_dataset_sst(filename, obs_variable, localize)
    character(len=*), intent(in) :: filename
    integer, intent(in) :: obs_variable
    logical, intent(in), optional :: localize

    integer, parameter :: max_levels = 1

    real :: lon, lat, rms_err
    real :: ri0, rj0
    real, dimension(MAX_LEVELS) :: depth, data

    integer :: unit, ndim, nvar, natt, ntime
    integer :: var_id, inst_type
    integer :: num_levs, k, kk, i, j, i0, j0
    integer :: stdout_unit

    logical :: data_is_local, localize_data
    logical, dimension(MAX_LEVELS) :: flag

    character(len=32) :: axisname, anal_fldname
    character(len=128) :: emsg_local

    type(axistype), dimension(:), allocatable, target :: axes
    type(axistype), pointer :: lon_axis, lat_axis

    if ( PRESENT(localize) ) then
       localize_data = localize
    else
       localize_data = .true.
    end if

    stdout_unit = stdout()

    anal_fldname = 'temp'
    var_id = obs_variable

    call mpp_open(unit, trim(filename), mpp_rdonly, mpp_netcdf, threading=mpp_multi, fileset=mpp_single)
    call mpp_get_info(unit, ndim, nvar, natt, ntime)
    allocate(axes(ndim))
    call mpp_get_axes(unit, axes)

    do i=1, ndim
       call mpp_get_atts(axes(i),name=axisname)
       select case (trim(axisname))
!!$       case ('GRIDLON_T') ! cm2 grids
!!$       case ('gridlon_t') ! cm2 grids
       case ('XT_OCEAN') ! cm2.5 grids
!!$       case ('lon') ! after 2008
          lon_axis => axes(i)
!!$       case ('GRIDLAT_T') ! cm2 grids
!!$       case ('gridlat_t') ! cm2 grids
       case ('YT_OCEAN') ! cm2 grids
!!$       case ('lat') ! for after 2008
          lat_axis => axes(i)
       end select
    end do

    call mpp_get_atts(lon_axis,len=nlon)
    call mpp_get_atts(lat_axis,len=nlat)
    if ( nlon /= 1440 .or. nlat /= 1070 ) then ! after 2008
       write (unit=emsg_local, fmt='("sst obs dim is not same as in file. nlon = ",I5,", nlat = ",I5)') nlon, nlat
       call error_mesg('oda_core_mod::open_profile_dataset_sst', trim(emsg_local), fatal)
    end if

    ! idealized
    do j=1, nlat
       do i=1, nlon
          lon = x_grid(i,j)
          lat = y_grid(i,j)
          rms_err = 0.5
          inst_type = 20
          data_is_local = .true.

          if ( lon .lt. 0.0 ) lon = lon + 360.0
          if ( lon .gt. 360.0 ) lon = lon - 360.0
          if ( lon .lt. 80.0 ) lon = lon + 360.0

          if ( lat < sst_obs_start_lat .or. lat > sst_obs_end_lat ) data_is_local = .false. ! at the final test

          if ( grd%mask(i,j,1) == 0 ) data_is_local = .false.

          if ( abs(lat) < 40.0 ) then
             if ( i/4*4 /= i .or. j/4*4 /= j ) data_is_local = .false.
          else if ( abs(lat) < 60.0 ) then
             if ( i/8*8 /= i .or. j/6*6 /= j ) data_is_local = .false.
          else
             if ( i/16*16 /= i .or. j/8*8 /= j ) data_is_local = .false.
          end if

          if ( data_is_local .and. (.NOT.localize_data) ) then
             if ( lat < 60.0 ) then ! regular grids
                ri0 = frac_index(lon, x_grid(:,nlat/2))
                rj0 = frac_index(lat, y_grid(nlon/4,:))
                i0 = floor(ri0)
                j0 = floor(rj0)
             else ! tripolar grids
                lon_out(1,1) = lon*deg_to_rad
                lat_out(1,1) = lat*deg_to_rad
                call horiz_interp_bilinear_new (interp, x_grid*deg_to_rad, y_grid*deg_to_rad, lon_out, lat_out)
                if ( interp%wti(1,1,2) < 1.0 ) then
                   i0 = interp%i_lon(1,1,1)
                else
                   i0 = interp%i_lon(1,1,2)
                end if
                if ( interp%wtj(1,1,2) < 1.0 ) then
                   j0 = interp%j_lat(1,1,1)
                else
                   j0 = interp%j_lat(1,1,2)
                end if

                if ( i0 > ieg .or. j0 > jeg ) then
                   write (unit=emsg_local, fmt='("i0 = ",I8,", j0 = ",I8)') i0, j0
                   call error_mesg('oda_core_mod::open_profile_dataset_sst',&
                        & 'For tripolar grids, either i0 > ieg or j0 > jeg.  '//trim(emsg_local), fatal)
                end if
             end if

             if ( i0 /= ieg .and. j0 /= jeg ) then ! exclude SSTs at ieg and jeg
                if ( grd%mask(i0,j0,1) /= 0.0 .and. grd%mask(i0+1,j0,1) /= 0.0 .and.&
                     & grd%mask(i0,j0+1,1) /= 0.0 .and. grd%mask(i0+1,j0+1,1) /= 0.0 ) then
                   no_sst = no_sst+1
                   num_profiles=num_profiles+1

                   if ( num_profiles > max_profiles ) then
                      call error_mesg('oda_core_mod::open_profile_dataset_sst',&
                           & 'Maximum number of profiles exceeded, increase max_profiles in oda_core_nml', fatal)
                   end if

                   num_levs = 0
                   flag = .false.
                   do k=1, 1
                      flag(k) = .true.
                      data(k) = 0.0
                      depth(k) = 0.0
                      num_levs = num_levs + 1
                   end do
                   if ( num_levs == 0 ) cycle
                   allocate(profiles(num_profiles)%depth(num_levs))
                   allocate(profiles(num_profiles)%data(num_levs))
                   allocate(profiles(num_profiles)%flag(num_levs))
                   allocate(profiles(num_profiles)%ms(num_levs))
                   allocate(profiles(num_profiles)%ms_inv(num_levs))
                   profiles(num_profiles)%variable = var_id
                   profiles(num_profiles)%inst_type = inst_type
                   profiles(num_profiles)%levels = num_levs
                   profiles(num_profiles)%lat = lat
                   profiles(num_profiles)%lon = lon

                   kk = 1
                   do k=1, 1
                      if ( flag(k) ) then
                         profiles(num_profiles)%depth(kk) = depth(k)
                         profiles(num_profiles)%data(kk) = data(k)
                         profiles(num_profiles)%ms(kk) = 1.0
                         profiles(num_profiles)%ms_inv(kk) = 1.0
                         kk=kk+1
                      end if
                   end do

                   ! calculate interpolation coefficients (make sure to account for grid offsets here!)
                   if ( lat < 60.0 ) then ! for regular grids
                      profiles(num_profiles)%i_index = ri0
                      profiles(num_profiles)%j_index = rj0
                   else ! for tripolar grids
                      lon_out(1,1) = lon*deg_to_rad
                      lat_out(1,1) = lat*deg_to_rad
                      call horiz_interp_bilinear_new (interp, x_grid*deg_to_rad, y_grid*deg_to_rad, lon_out, lat_out)
                      if ( interp%wti(1,1,2) < 1.0 ) then
                         profiles(num_profiles)%i_index =interp%i_lon(1,1,1) + interp%wti(1,1,2)
                      else
                         profiles(num_profiles)%i_index =interp%i_lon(1,1,2)
                      end if
                      if ( interp%wtj(1,1,2) < 1.0 ) then
                         profiles(num_profiles)%j_index =interp%j_lat(1,1,1) + interp%wtj(1,1,2)
                      else
                         profiles(num_profiles)%j_index =interp%j_lat(1,1,2)
                      end if
                   end if

                   profiles(num_profiles)%accepted = .true.
                   if ( i0 < 1 .or. j0 < 1 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                   if ( profiles(num_profiles)%accepted ) then
                      if ( grd%mask(i0,j0,1) == 0.0 .or.&
                           & grd%mask(i0+1,j0,1) == 0.0 .or.&
                           & grd%mask(i0,j0+1,1) == 0.0 .or.&
                           & grd%mask(i0+1,j0+1,1) == 0.0 ) then
                         profiles(num_profiles)%accepted = .false.
                      end if
                   end if
                   if ( profiles(num_profiles)%accepted ) then
                      profiles(num_profiles)%flag(:) = .true.
                      allocate(profiles(num_profiles)%k_index(profiles(num_profiles)%levels))
                      do k=1, profiles(num_profiles)%levels
                         profiles(num_profiles)%k_index(k) = frac_index(depth(k), (/0.,grd%z(:)/)) - 1
                         !::sdu:: Do we need the same out-of-range check here?
                      end do
                   end if
                end if ! exclude SSTs at ieg and jeg
             end if
          end if
       end do
    end do

    call mpp_close(unit)
    deallocate(axes)
    write (unit=stdout_unit, fmt='("A grand total of ",I8," sst points within global domain")') no_sst
    write (unit=stdout_unit, fmt='("A final total @sst of ",I8," prfs within global domain")') num_profiles
  end subroutine open_profile_dataset_sst

  subroutine open_profile_dataset_woa05t(filename, obs_variable, localize)
    character(len=*), intent(in) :: filename
    integer, intent(in) :: obs_variable
    logical, intent(in), optional :: localize

    integer, parameter :: MAX_LEVELS = 24

    real :: lon, lat, rms_err
    real :: ri0, rj0
    real, dimension(MAX_LEVELS) :: depth, data

    integer :: unit, ndim, nvar, natt, ntime
    integer :: var_id, inst_type
    integer :: num_levs, k, kk, i, j, i0, j0, k0
    integer :: stdout_unit, istat
    integer :: out_bound_point

    logical :: data_is_local, localize_data
    logical :: prof_in_filt_domain
    logical, dimension(MAX_LEVELS) :: flag

    character(len=32) :: axisname, anal_fldname
    character(len=128) :: emsg_local

    type(axistype), dimension(:), allocatable, target :: axes
    type(axistype), pointer :: lon_axis, lat_axis, z_axis, t_axis

    if ( PRESENT(localize) ) then
       localize_data = localize
    else
       localize_data = .true.
    end if

    stdout_unit = stdout()
    anal_fldname = 'temp'
    var_id = obs_variable

!    call mpp_print_memuse_stats('open_profile_dataset_woa05t Start')

    call mpp_open(unit, trim(filename), mpp_rdonly, mpp_netcdf, threading=mpp_multi, fileset=mpp_single)

    write (unit=stdout_unit, fmt='("Opened profile woa05t dataset: ",A)') trim(filename)

    call mpp_get_info(unit, ndim, nvar, natt, ntime)
    allocate(axes(ndim))
    call mpp_get_axes(unit,axes)
    do i=1, ndim
       call mpp_get_atts(axes(i), name=axisname)
       select case (trim(axisname))
       case ('lon')
          lon_axis => axes(i)
       case ('lat')
          lat_axis => axes(i)
       case ('depth')
          z_axis => axes(i)
       end select
    end do

    call mpp_get_atts(lon_axis, len=nlon_woa)
    call mpp_get_atts(lat_axis, len=nlat_woa)
    call mpp_get_atts(z_axis, len=nlev_woa)

    if ( nlon_woa /= 360 .or. nlat_woa /= 180 ) then
       write (unit=emsg_local, fmt='("woa05 obs dim is not same as in file. nlon_woa = ",I8,", nlat_woa = ",I8)') nlon_woa, nlat_woa
       call error_mesg('oda_core_mod::open_profile_dataset_woa05t', trim(emsg_local), fatal)
    end if

    allocate(woa05_lon(nlon_woa), woa05_lat(nlat_woa), woa05_z(nlev_woa))

    call mpp_get_axis_data(lon_axis, woa05_lon)
    call mpp_get_axis_data(lat_axis, woa05_lat)
    call mpp_get_axis_data(z_axis, woa05_z)

    out_bound_point = 0
    ! idealized
    do j=1, nlat_woa
       do i=1, nlon_woa
          lon = woa05_lon(i)
          lat = woa05_lat(j)
          rms_err = 5
          inst_type = 20
          data_is_local = .true.
          prof_in_filt_domain = .false.

          if ( lon .lt. 0.0 ) lon = lon + 360.0
          if ( lon .gt. 360.0 ) lon = lon - 360.0
          if ( lon .lt. 80.0 ) lon = lon + 360.0

          if ( lat < -80.0 .or. lat > 80.0 ) data_is_local = .false.
          if ( abs(lat) < 20.0 .and. (mod(i,2) /= 0 .or. mod(j,2) /= 0) ) data_is_local = .false.
          if ( abs(lat) >= 20.0 .and. (mod(i,4) /= 0 .or. mod(j,4) /= 0) ) data_is_local = .false.
          if ( abs(lat) >= 60.0 .and. (mod(i,6) /= 0 .or. mod(j,6) /= 0) ) data_is_local = .false.

          if (isd_filt < 1 .and. ied_filt > ieg) then
             ! filter domain is a full x band
             if (lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  lat >= y_grid(1,jsd_flt0) .and. lat <= y_grid(ieg-1,jsd_flt0)) then
                prof_in_filt_domain = .true.
             end if
          else if (isd_filt >= 1 .and. ied_filt <= ieg) then
             ! Interior filter domain
            if (lon >= x_grid(isd_filt,jsd_flt0) .and. lon <= x_grid(ied_filt-1,jsd_flt0) .and.&
               & lat >= y_grid(isd_filt,jsd_flt0) .and. lat <=  y_grid(ied_filt-1,jed_flt0-1)) then
                prof_in_filt_domain = .true.
            end if
          else if (isd_filt < 1 .and. ied_filt <= ieg) then
             ! lhs filter domain
             isd_flt0 = isd_filt + ieg
             if ((lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ied_filt-1,jsd_flt0) .and.&
                  & lat >= y_grid(1,jsd_flt0) .and. lat <=  y_grid(ied_filt-1,jed_flt0-1)).or.&
                  & (lon >= x_grid(isd_flt0,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  & lat >= y_grid(isd_flt0,jsd_flt0) .and. lat <=  y_grid(ieg-1,jed_flt0-1))) then
                prof_in_filt_domain = .true.
             end if
          else if (isd_filt >= 1 .and. ied_filt > ieg) then
             ! rhs filter domain
             ied_flt0 = ied_filt - ieg
             if ( lon >= x_grid(isd_filt,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  & lat >= y_grid(isd_filt,jsd_flt0) .and. lat <=  y_grid(ieg-1,jed_flt0-1) ) then
                prof_in_filt_domain = .true.
             end if
             if (ied_flt0-1 > 1) then
                if ( lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ied_flt0-1,jsd_flt0) .and.&
                     & lat >= y_grid(1,jsd_flt0) .and. lat <=  y_grid(ied_flt0-1,jed_flt0-1) ) then
                   prof_in_filt_domain = .true.
                end if
             end if
          end if

          if ( data_is_local .and. (.NOT.localize_data) ) then ! global index

            no_woa05 = no_woa05 + 1
            num_profiles = num_profiles + 1
             if ( num_profiles > max_profiles ) then
                call error_mesg('oda_core_mod::open_profile_dataset_woa05t',&
                     & 'Maximum number of profiles exceeded, increase max_profiles in oda_core_nml', fatal)
             end if

            num_levs = 0
            flag = .false.
            do k=1, nlev
              flag(k) = .true.
              data(k) = 0.0
              depth(k) = woa05_z(k)
              num_levs = num_levs+1
            end do
            if ( num_levs == 0 ) cycle

            if ( prof_in_filt_domain ) then ! localize

             allocate(profiles(num_profiles)%depth(num_levs))
             allocate(profiles(num_profiles)%data(num_levs))
             allocate(profiles(num_profiles)%flag(num_levs))
!             allocate(profiles(num_profiles)%ms(num_levs))
!             allocate(profiles(num_profiles)%ms_inv(num_levs))
             profiles(num_profiles)%variable = var_id
             profiles(num_profiles)%inst_type = inst_type
             profiles(num_profiles)%levels = num_levs
             profiles(num_profiles)%lat = lat
             profiles(num_profiles)%lon = lon

             kk = 1
             do k=1, nlev
                if ( flag(k) ) then
                   profiles(num_profiles)%depth(kk) = depth(k)
                   profiles(num_profiles)%data(kk) = data(k)
!                   profiles(num_profiles)%ms(kk) = 1.0
!                   profiles(num_profiles)%ms_inv(kk) = 1.0
                   kk = kk + 1
                end if
             end do

             ! calculate interpolation coefficients (make sure to account for grid offsets here!)
             if ( lat < 65.0 ) then ! regular grids
                ri0 = frac_index(lon, x_grid(:,1))
                rj0 = frac_index(lat, y_grid(90,:))
                i0 = floor(ri0)
                j0 = floor(rj0)
                if ( i0 > ieg .or. j0 > jeg ) then
                   write (unit=emsg_local, fmt='("i0 = ",I8,", j0 = ",I8)') i0, j0
                   call error_mesg('oda_core_mod::open_profile_dataset_woa05t',&
                        & 'For regular grids, either i0 > ieg or j0 > jeg.  '//trim(emsg_local), fatal)
                end if
                if ( isd_filt >= 1 .and. ied_filt <= ieg ) then
                  if ( i0 < isd_filt .or. i0 > ied_filt .or. j0 < jsd_filt .or. j0 > jed_filt ) then
                     write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                          & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                     call error_mesg('oda_core_mod::open_profile_dataset',&
                          & 'i0,j0 out of bounds in woat01. '//trim(emsg_local), fatal)
                  end if
                end if
                if ( isd_filt < 1 .and. i0 > ied_filt-1 .and. i0 < isd_filt + ieg ) then
                   write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                        & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                   call error_mesg('oda_core_mod::open_profile_dataset',&
                        & 'i0,j0 out of bounds in woat02. '//trim(emsg_local), fatal)
                end if
                if ( ied_filt > ieg .and. i0 > ied_filt-ieg-1 .and. ied_filt < isd_filt ) then
                   write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                        & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                   call error_mesg('oda_core_mod::open_profile_dataset',&
                        & 'i0,j0 out of bounds in woat03. '//trim(emsg_local), fatal)
                end if
                profiles(num_profiles)%i_index = ri0
                profiles(num_profiles)%j_index = rj0
             else ! tripolar grids
                lon_out(1,1) = lon*deg_to_rad
                lat_out(1,1) = lat*deg_to_rad
                call horiz_interp_bilinear_new (interp, x_grid*deg_to_rad, y_grid*deg_to_rad, lon_out, lat_out)
                if ( interp%wti(1,1,2) < 1.0 ) then
                   i0 = interp%i_lon(1,1,1)
                else
                   i0 = interp%i_lon(1,1,2)
                end if
                if ( interp%wtj(1,1,2) < 1.0 ) then
                   j0 = interp%j_lat(1,1,1)
                else
                   j0 = interp%j_lat(1,1,2)
                end if
                if ( i0 > ieg .or. j0 > jeg ) then
                   write (unit=emsg_local, fmt='("i0 = ",I8,", j0 = ",I8)') i0, j0
                   call error_mesg('oda_core_mod::open_profile_dataset_woa05t',&
                        & 'For tripolar grids, either i0 > ieg or j0 > jeg.  '//trim(emsg_local), fatal)
                end if
                if ( i0 < isd_filt .or. i0 > ied_filt .or. j0 < jsd_filt .or. j0 > jed_filt ) then
                   write (unit=stdout_unit, fmt='("woat.pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                        & mpp_pe(),i0, j0,isd_filt,ied_filt,jsd_filt,jed_filt
                   out_bound_point = out_bound_point + 1
                end if
                if ( interp%wti(1,1,2) < 1.0 ) then
                   profiles(num_profiles)%i_index =interp%i_lon(1,1,1) + interp%wti(1,1,2)
                else
                   profiles(num_profiles)%i_index =interp%i_lon(1,1,2)
                end if
                if ( interp%wtj(1,1,2) < 1.0 ) then
                   profiles(num_profiles)%j_index =interp%j_lat(1,1,1) + interp%wtj(1,1,2)
                else
                   profiles(num_profiles)%j_index =interp%j_lat(1,1,2)
                end if
             end if ! grids

             profiles(num_profiles)%accepted = .true.
             if ( i0 < 1 .or. j0 < 1 ) then
                profiles(num_profiles)%accepted = .false.
             end if
             if ( i0 < isd_filt .or. i0 >= ied_filt .or. j0 < jsd_filt .or. j0 >= jed_filt ) then
                profiles(num_profiles)%accepted = .false.
             end if

             if ( profiles(num_profiles)%accepted ) then ! here
                if ( i0 /= ieg .and. j0 /= jeg ) then
                   if ( grd%mask(i0,j0,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0,1) == 0.0 .or.&
                        & grd%mask(i0,j0+1,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0+1,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                else if ( i0 == ieg .and. j0 /= jeg ) then
                   if ( grd%mask(i0,j0,1) == 0.0 .or.&
                        & grd%mask(1,j0,1) == 0.0 .or.&
                        & grd%mask(i0,j0+1,1) == 0.0 .or.&
                        & grd%mask(1,j0+1,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                else if ( i0 /= ieg .and. j0 == jeg ) then
                   if ( grd%mask(i0,j0,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                else
                   if ( grd%mask(i0,j0,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                end if
             end if ! here

             if ( profiles(num_profiles)%accepted ) then ! accepted
                profiles(num_profiles)%flag(:) = .true.
                allocate(profiles(num_profiles)%k_index(profiles(num_profiles)%levels))
                do k=1, profiles(num_profiles)%levels
                   if (depth(k) < grd%z(1)) then
                     profiles(num_profiles)%k_index(k) = 0.0
                   else
                     profiles(num_profiles)%k_index(k) = frac_index(depth(k), (/0.,grd%z(:)/)) - 1.0 ! snz modify to v3.2 JAN3012
                   end if
                   if ( profiles(num_profiles)%k_index(k) > nk ) then
                      call error_mesg('oda_core_mod::open_profile_dataset_woa05t',&
                           & 'Profile k_index is greater than nk', fatal)
                   end if
                   k0 = floor(profiles(num_profiles)%k_index(k))

                   if ( k0 >= 1 ) then ! snz add
                      if ( profiles(num_profiles)%flag(k) ) then ! flag
                         if ( i0 /= ieg .and. j0 /= jeg ) then
                            if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0) == 0.0 .or.&
                                 & grd%mask(i0+1,j0+1,k0) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else if ( i0 == ieg .and. j0 /= jeg ) then
                            if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                                 & grd%mask(1,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0) == 0.0 .or.&
                                 & grd%mask(1,j0+1,k0) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else if ( i0 /= ieg .and. j0 == jeg ) then
                            if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else
                            if ( grd%mask(i0,j0,k0) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         end if

                         if ( i0 /= ieg .and. j0 /= jeg ) then
                            if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0+1) == 0.0 .or.&
                                 & grd%mask(i0+1,j0+1,k0+1) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else if ( i0 == ieg .and. j0 /= jeg ) then
                            if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(1,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0+1) == 0.0 .or.&
                                 & grd%mask(1,j0+1,k0+1) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else if ( i0 /= ieg .and. j0 == jeg ) then
                            if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0+1) == 0.0) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else
                            if ( grd%mask(i0,j0,k0+1) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         end if

                         if ( abs(profiles(num_profiles)%data(k)) > 1.e4 &
                              & .or. abs(profiles(num_profiles)%depth(k)) > 1.e4 ) then
                            profiles(num_profiles)%flag(k) = .false.
                         end if
                      end if ! flag
                   end if ! snz add
                end do
             end if ! accepted
           end if ! localize
         end if ! global index
       end do
    end do

    call mpp_close(unit)
    deallocate(axes)
    write (unit=stdout_unit, fmt='("A grand total of ",I8," woa05t points within global domain")') no_woa05
    write (unit=stdout_unit, fmt='("A final total @woa05t of ",I8," prfs within global domain")') num_profiles

!    call mpp_print_memuse_stats('open_profile_dataset_woa05t End')

  end subroutine open_profile_dataset_woa05t

  subroutine open_profile_dataset_woa05s(filename, obs_variable, localize)
    character(len=*), intent(in) :: filename
    integer, intent(in) :: obs_variable
    logical, intent(in), optional :: localize

    integer, parameter :: MAX_LEVELS = 24

    real :: lon, lat, rms_err
    real :: ri0, rj0
    real, dimension(MAX_LEVELS) :: depth, data

    integer :: unit, ndim, nvar, natt, ntime
    integer :: var_id, inst_type
    integer :: num_levs, k, kk, i, j, i0, j0, k0
    integer :: stdout_unit
    integer :: out_bound_point

    logical :: data_is_local, localize_data
    logical :: prof_in_filt_domain
    logical, dimension(MAX_LEVELS) :: flag

    character(len=32) :: axisname, anal_fldname
    character(len=128) :: emsg_local

    type(axistype), dimension(:), allocatable, target :: axes
    type(axistype), pointer :: lon_axis, lat_axis, z_axis

    if ( present(localize) ) then
       localize_data = localize
    else
       localize_data = .true.
    end if

    stdout_unit = stdout()

    anal_fldname = 'salt'
    var_id = obs_variable

!    call mpp_print_memuse_stats('open_profile_dataset_woa05s Start')

    call mpp_open(unit, trim(filename), mpp_rdonly, mpp_netcdf, threading=mpp_multi, fileset=mpp_single)

    write (unit=stdout_unit, fmt='("Opened profile woa05s dataset: ",A)') trim(filename)

    call mpp_get_info(unit, ndim, nvar, natt, ntime)
    allocate(axes(ndim))
    call mpp_get_axes(unit, axes)
    do i=1, ndim
       call mpp_get_atts(axes(i), name=axisname)
       select case ( trim(axisname) )
       case ('lon')
          lon_axis => axes(i)
       case ('lat')
          lat_axis => axes(i)
       case ('depth')
          z_axis => axes(i)
       end select
    end do
    call mpp_get_atts(lon_axis, len=nlon_woa)
    call mpp_get_atts(lat_axis, len=nlat_woa)
    call mpp_get_atts(z_axis, len=nlev_woa)
    if ( nlon_woa /= 360 .or. nlat_woa /= 180 ) then
       write (unit=emsg_local, fmt='("woa05 obs dim is not same as in file nlon_woa = ",I8,", nlat_woa = ",I8)') nlon_woa, nlat_woa
       call error_mesg('oda_core_mod::open_profile_dataset_woa05s', trim(emsg_local), fatal)
    end if

    out_bound_point = 0
    ! idealized
    do j=1, nlat_woa
       do i=1, nlon_woa
          lon = woa05_lon(i)
          lat = woa05_lat(j)
          rms_err = 5
          inst_type = 20
          data_is_local = .true.
          prof_in_filt_domain = .false.

          if ( lon .lt. 0.0 ) lon = lon + 360.0
          if ( lon .gt. 360.0 ) lon = lon - 360.0
          if ( lon .lt. 80.0 ) lon = lon + 360.0

          if ( lat < -80.0 .or. lat > 80.0 ) data_is_local = .false.
          if ( abs(lat) < 20.0 .and. (mod(i,2) /= 0 .or. mod(j,2) /= 0) ) data_is_local = .false.
          if ( abs(lat) >= 20.0 .and. (mod(i,4) /= 0 .or. mod(j,4) /= 0) ) data_is_local = .false.
          if ( abs(lat) >= 60.0 .and. (mod(i,6) /= 0 .or. mod(j,6) /= 0) ) data_is_local = .false.

          if (isd_filt < 1 .and. ied_filt > ieg) then
             ! filter domain is a full x band
             if (lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  lat >= y_grid(1,jsd_flt0) .and. lat <= y_grid(ieg-1,jsd_flt0)) then
                prof_in_filt_domain = .true.
             end if
          else if (isd_filt >= 1 .and. ied_filt <= ieg) then
             ! Interior filter domain
            if (lon >= x_grid(isd_filt,jsd_flt0) .and. lon <= x_grid(ied_filt-1,jsd_flt0) .and.&
               & lat >= y_grid(isd_filt,jsd_flt0) .and. lat <=  y_grid(ied_filt-1,jed_flt0-1)) then
                prof_in_filt_domain = .true.
            end if
          else if (isd_filt < 1 .and. ied_filt <= ieg) then
             ! lhs filter domain
             isd_flt0 = isd_filt + ieg
             if ((lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ied_filt-1,jsd_flt0) .and.&
                  & lat >= y_grid(1,jsd_flt0) .and. lat <=  y_grid(ied_filt-1,jed_flt0-1)).or.&
                  & (lon >= x_grid(isd_flt0,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  & lat >= y_grid(isd_flt0,jsd_flt0) .and. lat <=  y_grid(ieg-1,jed_flt0-1))) then
                prof_in_filt_domain = .true.
             end if
          else if (isd_filt >= 1 .and. ied_filt > ieg) then
             ! rhs filter domain
             ied_flt0 = ied_filt - ieg
             if ( lon >= x_grid(isd_filt,jsd_flt0) .and. lon <= x_grid(ieg-1,jsd_flt0) .and.&
                  & lat >= y_grid(isd_filt,jsd_flt0) .and. lat <=  y_grid(ieg-1,jed_flt0-1) ) then
                prof_in_filt_domain = .true.
             end if
             if (ied_flt0-1 > 1) then
                if ( lon >= x_grid(1,jsd_flt0) .and. lon <= x_grid(ied_flt0-1,jsd_flt0) .and.&
                     & lat >= y_grid(1,jsd_flt0) .and. lat <=  y_grid(ied_flt0-1,jed_flt0-1) ) then
                   prof_in_filt_domain = .true.
                end if
             end if
          end if

          if ( data_is_local .and. (.NOT.localize_data) ) then ! global index

            no_woa05 = no_woa05 + 1
            num_profiles=num_profiles + 1
             if ( num_profiles > max_profiles ) then
                call error_mesg('oda_core_mod::open_profile_dataset_woa05s',&
                     & 'Maximum number of profiles exceeded, increase max_profiles in oda_core_nml.', fatal)
             end if

            num_levs = 0
            flag = .false.
            do k=1, nlev_woa
              flag(k) = .true.
              data(k) = 0.0
              depth(k) = woa05_z(k)
              num_levs = num_levs + 1
            end do
            if ( num_levs == 0 ) cycle

           if ( prof_in_filt_domain ) then ! localize

             allocate(profiles(num_profiles)%depth(num_levs))
             allocate(profiles(num_profiles)%data(num_levs))
             allocate(profiles(num_profiles)%flag(num_levs))
!             allocate(profiles(num_profiles)%ms(num_levs))
!             allocate(profiles(num_profiles)%ms_inv(num_levs))
             profiles(num_profiles)%variable = var_id
             profiles(num_profiles)%inst_type = inst_type
             profiles(num_profiles)%levels = num_levs
             profiles(num_profiles)%lat = lat
             profiles(num_profiles)%lon = lon

             kk = 1
             do k=1, nlev_woa
                if ( flag(k) ) then
                   profiles(num_profiles)%depth(kk) = depth(k)
                   profiles(num_profiles)%data(kk) = data(k)
!                   profiles(num_profiles)%ms(kk) = 1.0
!                   profiles(num_profiles)%ms_inv(kk) = 1.0
                   kk = kk + 1
                end if
             end do

             ! calculate interpolation coefficients (make sure to account for grid offsets here!)
             if ( lat < 65.0 ) then ! regular grids
                ri0 = frac_index(lon, x_grid(:,1))
                rj0 = frac_index(lat, y_grid(90,:))
                i0 = floor(ri0)
                j0 = floor(rj0)
                if ( i0 > ieg .or. j0 > jeg ) then
                   write (unit=emsg_local, fmt='("i0 = ",I8,", j0 = ",I8)') i0, j0
                   call error_mesg('oda_core_mod::open_profile_dataset_woa05s',&
                        & 'For regular grids, either i0 > ieg or j0 > jeg.  '//trim(emsg_local), fatal)
                end if
                if ( isd_filt >= 1 .and. ied_filt <= ieg ) then
                  if ( i0 < isd_filt .or. i0 > ied_filt .or. j0 < jsd_filt .or. j0 > jed_filt ) then
                     write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                          & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                     call error_mesg('oda_core_mod::open_profile_dataset',&
                          & 'i0,j0 out of bounds in woas01. '//trim(emsg_local), fatal)
                  end if
                end if
                if ( isd_filt < 1 .and. i0 > ied_filt-1 .and. i0 < isd_filt + ieg ) then
                   write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                        & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                   call error_mesg('oda_core_mod::open_profile_dataset',&
                        & 'i0,j0 out of bounds in woas02. '//trim(emsg_local), fatal)
                end if
                if ( ied_filt > ieg .and. i0 > ied_filt-ieg-1 .and. ied_filt < isd_filt ) then
                   write (unit=emsg_local, fmt='("pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                        & mpp_pe(), i0, j0, isd_filt,ied_filt,jsd_filt,jed_filt
                   call error_mesg('oda_core_mod::open_profile_dataset',&
                        & 'i0,j0 out of bounds in woas03. '//trim(emsg_local), fatal)
                end if
                profiles(num_profiles)%i_index = ri0
                profiles(num_profiles)%j_index = rj0
             else ! tripolar grids
                lon_out(1,1) = lon*deg_to_rad
                lat_out(1,1) = lat*deg_to_rad
                call horiz_interp_bilinear_new (interp, x_grid*deg_to_rad, y_grid*deg_to_rad, lon_out, lat_out)
                if ( interp%wti(1,1,2) < 1.0 ) then
                   i0 = interp%i_lon(1,1,1)
                else
                   i0 = interp%i_lon(1,1,2)
                end if
                if ( interp%wtj(1,1,2) < 1.0 ) then
                   j0 = interp%j_lat(1,1,1)
                else
                   j0 = interp%j_lat(1,1,2)
                end if
                if ( i0 > ieg .or. j0 > jeg ) then
                   write (unit=emsg_local, fmt='("i0 = ",I8,", j0 = ",I8)') i0, j0
                   call error_mesg('oda_core_mod::open_profile_dataset_woa05s',&
                        & 'For tripolar grids, either i0 > ieg or j0 > jeg.  '//trim(emsg_local), fatal)
                end if
                if ( i0 < isd_filt .or. i0 > ied_filt .or. j0 < jsd_filt .or. j0 > jed_filt ) then
                   write (unit=stdout_unit, fmt='("woas.pe,i0,j0= ",3I8,"isd_filt,ied_filt,jsd_filt,jed_filt= ",4I8)')&
                        & mpp_pe(),i0, j0,isd_filt,ied_filt,jsd_filt,jed_filt
                   out_bound_point = out_bound_point + 1
                end if
                if ( interp%wti(1,1,2) < 1.0 ) then
                   profiles(num_profiles)%i_index =interp%i_lon(1,1,1) + interp%wti(1,1,2)
                else
                   profiles(num_profiles)%i_index =interp%i_lon(1,1,2)
                end if
                if ( interp%wtj(1,1,2) < 1.0 ) then
                   profiles(num_profiles)%j_index =interp%j_lat(1,1,1) + interp%wtj(1,1,2)
                else
                   profiles(num_profiles)%j_index =interp%j_lat(1,1,2)
                end if
             end if ! grids

             profiles(num_profiles)%accepted = .true.
             if ( i0 < 1 .or. j0 < 1 ) then
                profiles(num_profiles)%accepted = .false.
             end if
             if ( i0 < isd_filt .or. i0 >= ied_filt .or. j0 < jsd_filt .or. j0 >= jed_filt ) then
                profiles(num_profiles)%accepted = .false.
             end if

             if ( profiles(num_profiles)%accepted ) then ! here
                if ( i0 /= ieg .and. j0 /= jeg ) then
                   if ( grd%mask(i0,j0,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0,1) == 0.0 .or.&
                        & grd%mask(i0,j0+1,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0+1,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                else if ( i0 == ieg .and. j0 /= jeg ) then
                   if ( grd%mask(i0,j0,1) == 0.0 .or.&
                        & grd%mask(1,j0,1) == 0.0 .or.&
                        & grd%mask(i0,j0+1,1) == 0.0 .or.&
                        & grd%mask(1,j0+1,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                else if ( i0 /= ieg .and. j0 == jeg ) then
                   if ( grd%mask(i0,j0,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                else
                   if ( grd%mask(i0,j0,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                end if
             end if ! here

             if ( profiles(num_profiles)%accepted ) then ! accepted
                profiles(num_profiles)%flag(:) = .true.
                allocate(profiles(num_profiles)%k_index(profiles(num_profiles)%levels))
                do k=1, profiles(num_profiles)%levels
                   if (depth(k) < grd%z(1)) then
                     profiles(num_profiles)%k_index(k) = 0.0
                   else
                     profiles(num_profiles)%k_index(k) = frac_index(depth(k), (/0.,grd%z(:)/)) - 1.0 ! snz modify to v3.2 JAN3012
                   end if
                   if ( profiles(num_profiles)%k_index(k) > nk ) then
                      call error_mesg('oda_core_mod::open_profile_dataset_woa05s',&
                           & 'Profile k_index is greater than nk', fatal)
                   end if
                   k0 = floor(profiles(num_profiles)%k_index(k))

                   if ( k0 >= 1 ) then ! snz add
                      if ( profiles(num_profiles)%flag(k) ) then ! flag
                         if ( i0 /= ieg .and. j0 /= jeg ) then
                            if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0) == 0.0 .or.&
                                 & grd%mask(i0+1,j0+1,k0) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else if ( i0 == ieg .and. j0 /= jeg ) then
                            if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                                 & grd%mask(1,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0) == 0.0 .or.&
                                 & grd%mask(1,j0+1,k0) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else if ( i0 /= ieg .and. j0 == jeg ) then
                            if ( grd%mask(i0,j0,k0) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else
                            if ( grd%mask(i0,j0,k0) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         end if

                         if ( i0 /= ieg .and. j0 /= jeg ) then
                            if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0+1) == 0.0 .or.&
                                 & grd%mask(i0+1,j0+1,k0+1) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         elseif ( i0 == ieg .and. j0 /= jeg ) then
                            if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(1,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0,j0+1,k0+1) == 0.0 .or.&
                                 & grd%mask(1,j0+1,k0+1) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else if ( i0 /= ieg .and. j0 == jeg ) then
                            if ( grd%mask(i0,j0,k0+1) == 0.0 .or.&
                                 & grd%mask(i0+1,j0,k0+1) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         else
                            if ( grd%mask(i0,j0,k0+1) == 0.0 ) then
                               profiles(num_profiles)%flag(k) = .false.
                            end if
                         end if

                         if ( abs(profiles(num_profiles)%data(k)) > 1.e4 &
                              & .or. abs(profiles(num_profiles)%depth(k)) > 1.e4 ) then
                            profiles(num_profiles)%flag(k) = .false.
                         end if
                      end if ! flag
                   end if ! snz add
                end do
             end if ! accepted
           end if ! localize
         end if ! global index
       end do
    end do

    call mpp_close(unit)
    deallocate(axes)
    write (unit=stdout_unit, fmt='("A grand total of ",I8," woa05s points within global domain")') no_woa05
    write (unit=stdout_unit, fmt='("A final total @woa05s of ",I8," prfs within global domain")') num_profiles

!    call mpp_print_memuse_stats('open_profile_dataset_woa05s Ens')

  end subroutine open_profile_dataset_woa05s

  subroutine get_obs_sst(model_time, Prof, nprof, no_prf0, sst_climo, Filter_domain)
    type(time_type), intent(in) :: model_time
    type(ocean_profile_type), dimension(:), intent(inout) :: prof
    integer, intent(inout) :: nprof
    integer, intent(in) :: no_prf0
    type(obs_clim_type), intent(inout) :: sst_climo
    type(domain2d), intent(in) :: filter_domain

    integer :: i0, j0, i, days, seconds, days1, seconds1
    integer :: unit, ndim, nvar, natt, ntime, time_idx
    integer :: iy0, in0, id0, ih0, im0, is0, i_m
    integer :: stdout_unit, istat
    integer, dimension(12) :: n_days
    integer, save :: year_on_first_read = 0 !< Year on first read of file during
                                            !! module run
    character(len=128) :: sst_filename, emsg_local

    type(fieldtype), dimension(:), allocatable :: fields
    type(time_type) :: tdiff, sst_time0, time1

    n_days  = (/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)

    nprof = 0

    stdout_unit = stdout()

    call get_time(model_time, seconds, days)
    call get_date(model_time, iy0, in0, id0, ih0, im0, is0)
    if ( year_on_first_read == 0 ) then
       year_on_first_read = iy0
    end if
    time1=set_date(year_on_first_read, 1, 1, 0, 0, 0)
    call get_time(time1, seconds1, days1)
    time_idx = days-days1+1

    ! daily data

    if ( mpp_pe() == mpp_root_pe() ) then
       write (unit=stdout_unit, fmt='("time_idx = ",I8)') time_idx
    end if

    sst_filename = "INPUT/sst_daily.nc"

    call mpp_open(unit, trim(sst_filename), mpp_rdonly, mpp_netcdf, threading=mpp_multi, fileset=mpp_single)
    call mpp_get_info(unit, ndim, nvar, natt, ntime)

    allocate(fields(nvar), stat=istat)
    if ( istat .ne. 0 ) then
       call error_mesg('oda_core_mod::get_obs_sst', 'Unable to allocate fields', fatal)
    end if

    call mpp_get_fields(unit, fields)
    do i=1, nvar
       select case ( mpp_get_field_name(fields(i)) )
       case ('SST1') ! for AVHRR daily SST
          call mpp_read(unit, fields(i), filter_domain, sst_climo%sst_obs, tindex=time_idx)
       end select
    end do

    ! get profiles and sst
    ! obs relevant to current analysis interval
    sst_time0 = set_date(iy0, in0, id0, ih0, im0, is0)

    if ( no_sst > 1 ) then

    do i=no_prf+no_woa05+1, no_prf+no_woa05+no_sst
       profiles(i)%time = sst_time0

       tdiff = model_time - profiles(i)%time

       i0 = floor(profiles(i)%i_index)
       if ( i0 < 1 .or. i0 > 1440 ) then
          write (unit=emsg_local, fmt='("Profiles(",I8,")%lon = ",I4,", i0 = ",I8)') i, profiles(i)%lon, i0
          call error_mesg('oda_core_mod::get_obs_sst',&
               & 'Profile longitude index outside range [1,1440].  '//trim(emsg_local), fatal)
       end if

       j0 = floor(profiles(i)%j_index)
       if ( j0 < 1 .or. j0 > 1070 ) then
          write (unit=emsg_local, fmt='("Profiles(",I8,")%lat = ",I4,", j0 = ",I8)') i, profiles(i)%lon, j0
          call error_mesg('oda_core_mod::get_obs_sst',&
               & 'Profile latitude index outside range [1,1070].  '//trim(emsg_local), fatal)
       end if

       if ( profiles(i)%accepted ) then
          nprof = nprof + 1
          if ( nprof > size(prof,1) ) then
             call error_mesg('oda_core_mod::get_obs_sst',&
                  & 'Passed in array "Prof" is smaller than number of profiles, increase size of Prof before call.',&
                  & fatal)
          end if

          call copy_obs(profiles(i:i), prof(no_prf0+nprof:no_prf0+nprof))

          prof(nprof+no_prf0)%tdiff = tdiff
       end if
    end do

    end if ! for no_sst > 1

    if ( mpp_pe() == mpp_root_pe() ) then
       write (unit=stdout_unit, fmt='("no of sst records: ",I8)') nprof
    end if

    deallocate(fields)
    call mpp_close(unit)
  end subroutine get_obs_sst

  subroutine get_obs_woa05t(model_time, Prof, nprof, no_prf0)
    type(time_type), intent(in) :: model_time
    type(ocean_profile_type), dimension(:), intent(inout) :: prof
    integer, intent(inout) :: nprof
    integer, intent(in) :: no_prf0

    real :: ri0, rj0, lon_woa05

    integer :: i0, j0
    integer :: i, k, unit, time_idx
    integer :: iy0, in0, id0, ih0, im0, is0
    integer :: ndim, nvar, natt, ntime
    integer :: stdout_unit, istat

    character(len=32) :: axisname
    character(len=128) :: woa05t_filename

    type(fieldtype), dimension(:), allocatable :: fields
    type(time_type) :: tdiff, woa05_time0

    nprof = 0
    stdout_unit = stdout()

    call get_date(model_time, iy0, in0, id0, ih0, im0, is0)

    time_idx = in0

    ! daily data
    woa05t_filename = "INPUT/woa05_temp.nc"

    call mpp_open(unit, trim(woa05t_filename), mpp_rdonly, mpp_netcdf, threading=mpp_multi, fileset=mpp_single)
    call mpp_get_info(unit, ndim, nvar, natt, ntime)

    allocate(fields(nvar), stat=istat)
    if ( istat .ne. 0 ) then
       call error_mesg('oda_core_mod::get_obs_woa05t', 'Unable to allocate fields', fatal)
    end if

    call mpp_get_fields(unit, fields)

    allocate(obs_woa05t(nlon_woa,nlat_woa,nlev_woa))

    do i=1, nvar
       select case ( mpp_get_field_name(fields(i)) )
       case ('t0112an1')
          call mpp_read(unit, fields(i), obs_woa05t, tindex=time_idx)
       end select
    end do

    woa05_time0 = set_date(iy0, in0, id0, ih0, im0, is0)

    do i=no_prf+1, no_prf+no_woa05/2
       profiles(i)%time = woa05_time0

       tdiff = model_time - profiles(i)%time

       lon_woa05 = profiles(i)%lon
       if ( lon_woa05 < 0.0 ) lon_woa05 = lon_woa05 + 360.0
       if ( lon_woa05 > 360.0 ) lon_woa05 = lon_woa05 - 360.0
       ri0 = frac_index(lon_woa05, woa05_lon)
       i0 = floor(ri0)
       if ( i0 < 1 ) i0 = 1
       if ( i0 > nlon_woa ) i0 = nlon_woa

       rj0 = frac_index(profiles(i)%lat, woa05_lat)
       j0 = floor(rj0)
       if(j0 < 1 ) j0 = 1
       if(j0 > nlat_woa) j0 = nlat_woa

       if ( profiles(i)%accepted ) then
          nprof = nprof + 1
          if ( nprof > size(prof,1) ) then
             call error_mesg('oda_core_mod::get_obs_woa05t',&
                  & 'Passed in array "Prof" is smaller than number of profiles, increase size of Prof before call.',&
                  & fatal)
          end if
          profiles(i)%data(1:nlev_woa) = obs_woa05t(i0,j0,1:nlev_woa)
          do k=1, nlev_woa
             if ( abs(profiles(i)%data(k)) > 1.e3 .or.&
                  & abs(profiles(i)%depth(k)) > 1.e5 ) then
                profiles(i)%flag(k) = .false.
             end if
          end do
          call copy_obs(profiles(i:i), prof(no_prf0+nprof:no_prf0+nprof))
          prof(no_prf0+nprof)%tdiff = tdiff
       end if
    end do

    if ( mpp_pe() == mpp_root_pe() ) then
       write (unit=stdout_unit, fmt='("no of woa05t records: ",I8)') nprof
    end if

    deallocate(fields, obs_woa05t)
    call mpp_close(unit)
  end subroutine get_obs_woa05t

  subroutine get_obs_woa05s(model_time, Prof, nprof, no_prf0)
    type(time_type), intent(in) :: model_time
    type(ocean_profile_type), dimension(:), intent(inout) :: prof
    integer, intent(inout) :: nprof
    integer, intent(in) :: no_prf0

    real :: ri0, rj0, lon_woa05

    integer :: i0, j0
    integer :: i, k, unit, time_idx
    integer :: iy0, in0, id0, ih0, im0, is0
    integer :: ndim, nvar, natt, ntime
    integer :: stdout_unit, istat

    character(len=128) :: woa05s_filename

    type(fieldtype), dimension(:), allocatable :: fields
    type(time_type) :: tdiff, woa05_time0

    nprof = 0
    stdout_unit = stdout()

    call get_date(model_time, iy0,in0,id0,ih0,im0,is0)

    time_idx = in0

    ! climatological data
    woa05s_filename = "INPUT/woa05_salt.nc"

    call mpp_open(unit, trim(woa05s_filename), mpp_rdonly, mpp_netcdf, threading=mpp_multi, fileset=mpp_single)
    call mpp_get_info(unit, ndim, nvar, natt, ntime)

    allocate(fields(nvar), stat=istat)
    if ( istat .ne. 0 ) then
       call error_mesg('oda_core_mod::get_obs_woa05s', 'Unable to allocate fields', fatal)
    end if

    call mpp_get_fields(unit, fields)

    allocate(obs_woa05s(nlon_woa,nlat_woa,nlev_woa))

    do i=1, nvar
       select case ( mpp_get_field_name(fields(i)) )
       case ('s0112an1')
          call mpp_read(unit, fields(i), obs_woa05s, tindex=time_idx)
       end select
    end do

    woa05_time0 = set_date(iy0, in0, id0, ih0, im0, is0)

    do i=no_prf+no_woa05/2+1, no_prf+no_woa05
       profiles(i)%time = woa05_time0

       tdiff = model_time - profiles(i)%time

       lon_woa05 = profiles(i)%lon
       if ( lon_woa05 < 0.0 ) lon_woa05 = lon_woa05 + 360.0
       if ( lon_woa05 > 360.0 ) lon_woa05 = lon_woa05 - 360.0
       ri0 = frac_index(lon_woa05, woa05_lon)
       i0 = floor(ri0)
       if ( i0 < 1 ) i0 = 1
       if ( i0 > nlon_woa ) i0 = nlon_woa

       rj0 = frac_index(profiles(i)%lat, woa05_lat)
       j0 = floor(rj0)
       if ( j0 < 1 ) j0 = 1
       if ( j0 > nlat_woa ) j0 = nlat_woa

       if ( profiles(i)%accepted ) then
          nprof = nprof + 1
          if ( nprof > size(prof,1) ) then
             call error_mesg('oda_core_mod::get_obs_woa05s',&
                  & 'Passed in array "Prof" is smaller than number of profiles, increase size of Prof before call',&
                  & fatal)
          end if
          profiles(i)%data(1:nlev_woa) = obs_woa05s(i0,j0,1:nlev_woa)
          do k=1, nlev_woa
             if ( abs(profiles(i)%data(k)) > 1.e3 .or.&
                  & abs(profiles(i)%depth(k)) > 1.e5 ) then
                profiles(i)%flag(k) = .false.
             end if
          end do
          call copy_obs(profiles(i:i),prof(no_prf0+nprof:no_prf0+nprof))
          prof(no_prf0+nprof)%tdiff = tdiff
       end if
    end do

    if ( mpp_pe() == mpp_root_pe() ) then
       write (unit=stdout_unit, fmt='("no of woa05t records: ",I8)') nprof
    end if

    deallocate(fields, obs_woa05s)
    call mpp_close(unit)
  end subroutine get_obs_woa05s

  subroutine open_profile_dataset_eta(filename, obs_variable, localize)
    character(len=*), intent(in) :: filename
    integer, intent(in) :: obs_variable
    logical, intent(in), optional :: localize

    integer, parameter :: MAX_LEVELS = 1000

    real :: lon, lat, rms_err
    real :: ri0, rj0
    real, dimension(MAX_LEVELS) :: depth, data

    integer :: inst_type, var_id
    integer :: num_levs, k, kk, i, j, i0, j0
    integer :: stdout_unit

    logical :: data_is_local, localize_data
    logical, dimension(MAX_LEVELS) :: flag

    character(len=32) :: anal_fldname

    if ( PRESENT(localize) ) then
       localize_data = localize
    else
       localize_data = .true.
    end if

    stdout_unit = stdout()

    anal_fldname = 'eta'
    var_id = obs_variable

    !snz idealized
    do j=1, size(x_grid, dim=1)
       do i=1, size(x_grid, dim=2)
          lon = x_grid(i,j)
          lat = y_grid(i,j)
          rms_err = 0.5

          inst_type = 20
          data_is_local = .true.

          if ( lon .lt. 0.0 ) lon = lon + 360.0
          if ( lon .gt. 360.0 ) lon = lon - 360.0
          if ( lon .lt. 80.0 ) lon = lon + 360.0

          if ( lat < eta_obs_start_lat .or. lat > eta_obs_end_lat ) data_is_local = .false.
          if ( abs(lat) < 20.0 .and.&
               & (mod(floor(lon),2) /= 0 .or. mod(floor(lat),2) /= 0) ) data_is_local = .false.
          if ( (abs(lat) >= 20.0 .and. abs(lat) < 40.0) .and.&
               & (mod(floor(lon),4) /= 0 .or. mod(floor(lat),4) /= 0) ) data_is_local = .false.
          if ( (abs(lat) >= 40.0) .and.&
               & (mod(floor(lon),6) /= 0 .or. mod(floor(lat),6) /= 0) ) data_is_local = .false.

          if ( data_is_local .and. (.NOT.localize_data) ) then
             ri0 = frac_index(lon, x_grid(:,1))
             rj0 = frac_index(lat, y_grid(90,:))
             i0 = floor(ri0)
             j0 = floor(rj0)
             if ( grd%mask(i0,j0,1) /= 0.0 .and. grd%mask(i0+1,j0,1) /= 0.0 .and.&
                  & grd%mask(i0,j0+1,1) /= 0.0 .and. grd%mask(i0+1,j0+1,1) /= 0.0 ) then
                no_eta = no_eta+1
                num_profiles=num_profiles+1
                if ( num_profiles > max_profiles ) then
                   call error_mesg('oda_core_mod::open_profile_dataset_eta',&
                        & 'Maximum number of profiles exceeded, increase max_profiles in oda_core_nml.', fatal)
                end if


                num_levs = 0
                flag = .false.
                do k=1, 1
                   flag(k) = .true.
                   data(k) = 0.0
                   depth(k) = 0.0
                   num_levs = num_levs+1
                end do
                if ( num_levs == 0 ) cycle
                allocate(profiles(num_profiles)%depth(num_levs))
                allocate(profiles(num_profiles)%data(num_levs))
                allocate(profiles(num_profiles)%flag(num_levs))
                allocate(profiles(num_profiles)%ms(num_levs))
                allocate(profiles(num_profiles)%ms_inv(num_levs))
                profiles(num_profiles)%variable = var_id
                profiles(num_profiles)%inst_type = inst_type
                profiles(num_profiles)%levels = num_levs
                profiles(num_profiles)%lat = lat
                profiles(num_profiles)%lon = lon
                kk = 1
                do k=1, 1
                   if ( flag(k) ) then
                      profiles(num_profiles)%depth(kk) = depth(k)
                      profiles(num_profiles)%data(kk) = data(k)
                      profiles(num_profiles)%ms(kk) = 1.0
                      profiles(num_profiles)%ms_inv(kk) = 1.0
                      kk = kk + 1
                   end if
                end do

                ! calculate interpolation coefficients (make sure to account for grid offsets here!)
                profiles(num_profiles)%i_index = ri0
                profiles(num_profiles)%j_index = rj0
                profiles(num_profiles)%accepted = .true.
                if ( i0 < 1 .or. j0 < 1 ) then
                   profiles(num_profiles)%accepted = .false.
                end if
                if ( profiles(num_profiles)%accepted ) then
                   if ( grd%mask(i0,j0,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0,1) == 0.0 .or.&
                        & grd%mask(i0,j0+1,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0+1,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                end if
                if ( profiles(num_profiles)%accepted ) then
                   profiles(num_profiles)%flag(:) = .true.
                   allocate(profiles(num_profiles)%k_index(profiles(num_profiles)%levels))
                   do k=1, profiles(num_profiles)%levels
                      if (depth(k) < grd%z(1)) then
                        profiles(num_profiles)%k_index(k) = 0.0
                      else
                        profiles(num_profiles)%k_index(k) = frac_index(depth(k), (/0.,grd%z(:)/)) - 1
                      end if
                      if ( profiles(num_profiles)%k_index(k) < 1.0 ) profiles(num_profiles)%flag(k) = .false.
                   end do
                end if
             end if
          end if
       end do
    end do

    write (unit=stdout_unit, fmt='("A grand total of ",I8," eta points within global domain")') no_eta
    write (unit=stdout_unit, fmt='("A final total @eta of ",I8," prfs within global domain")') num_profiles
  end subroutine open_profile_dataset_eta

  subroutine open_profile_dataset_suv(filename, obs_variable, localize)
    character(len=*), intent(in) :: filename
    integer, intent(in) :: obs_variable
    logical, intent(in), optional :: localize

    integer, parameter :: MAX_LEVELS = 1000

    real :: ri0, rj0
    real :: lon, lat, rms_err
    real, dimension(MAX_LEVELS) :: depth, data

    integer :: inst_type, var_id
    integer :: num_levs, k, kk, i, j, i0, j0
    integer :: stdout_unit

    logical :: data_is_local, localize_data
    logical, dimension(MAX_LEVELS) :: flag

    character(len=32) :: anal_fldname

    if ( PRESENT(localize) ) then
       localize_data = localize
    else
       localize_data = .true.
    end if

    stdout_unit = stdout()

    anal_fldname = 'suv'
    var_id = obs_variable

    !snz idealized
    do j=1, size(x_grid_uv,dim=1)
       do i=1, size(x_grid_uv,dim=2)
          lon = x_grid_uv(i,j)
          lat = y_grid_uv(i,j)
          rms_err = 0.5

          inst_type = 20
          data_is_local = .true.

          if ( lon .lt. 0.0 ) lon = lon + 360.0
          if ( lon .gt. 360.0 ) lon = lon - 360.0
          if ( lon .lt. 80.0 ) lon = lon + 360.0

          if ( lat < -40.0 .or. lat > 40.0 ) data_is_local = .false.
          if ( abs(lat) < 20.0 .and.&
               & (mod(floor(lon),2) /= 0 .or. mod(floor(lat),2) /= 0) ) data_is_local = .false.
          if ( (abs(lat) >= 20.0 .and. abs(lat) < 40.0) .and.&
               & (mod(floor(lon),4) /= 0 .or. mod(floor(lat),4) /= 0) ) data_is_local = .false.
          if ( (abs(lat) >= 40.0) .and.&
               & (mod(floor(lon),6) /= 0 .or. mod(floor(lat),6) /= 0) ) data_is_local = .false.

          if ( data_is_local .and. (.NOT.localize_data) ) then
             ri0 = frac_index(lon, x_grid_uv(:,1))
             rj0 = frac_index(lat, y_grid_uv(90,:))
             i0 = floor(ri0)
             j0 = floor(rj0)
             if ( grd%mask(i0,j0,1) /= 0.0 .and. grd%mask(i0+1,j0,1) /= 0.0 .and.&
                  & grd%mask(i0,j0+1,1) /= 0.0 .and. grd%mask(i0+1,j0+1,1) /= 0.0 ) then
                no_suv = no_suv+1
                num_profiles = num_profiles + 1
                if ( num_profiles > max_profiles ) then
                   call error_mesg('oda_core_mod::open_profile_dataset_suv',&
                        & 'Maximum number of profiles exceeded, increase max_profiles in oda_core_nml.', fatal)
                end if


                num_levs = 0
                flag = .false.
                do k=1, 2
                   flag(k) = .true.
                   data(k) = 0.0
                   depth(k) = 0.0
                   num_levs = num_levs + 1
                end do
                if ( num_levs == 0 ) cycle
                allocate(profiles(num_profiles)%depth(num_levs))
                allocate(profiles(num_profiles)%data(num_levs))
                allocate(profiles(num_profiles)%flag(num_levs))
                allocate(profiles(num_profiles)%ms(num_levs))
                allocate(profiles(num_profiles)%ms_inv(num_levs))
                profiles(num_profiles)%variable = var_id
                profiles(num_profiles)%inst_type = inst_type
                profiles(num_profiles)%levels = num_levs
                profiles(num_profiles)%lat = lat
                profiles(num_profiles)%lon = lon
                kk = 1
                do k=1, 2
                   if ( flag(k) ) then
                      profiles(num_profiles)%depth(kk) = depth(k)
                      profiles(num_profiles)%data(kk) = data(k)
                      profiles(num_profiles)%ms(kk) = 1.0
                      profiles(num_profiles)%ms_inv(kk) = 1.0
                      kk = kk + 1
                   end if
                end do

                ! calculate interpolation coefficients (make sure to account for grid offsets here!)
                profiles(num_profiles)%i_index = ri0
                profiles(num_profiles)%j_index = rj0
                profiles(num_profiles)%accepted = .true.
                if ( i0 < 1 .or. j0 < 1 ) then
                   profiles(num_profiles)%accepted = .false.
                end if
                if ( profiles(num_profiles)%accepted ) then
                   if ( grd%mask(i0,j0,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0,1) == 0.0 .or.&
                        & grd%mask(i0,j0+1,1) == 0.0 .or.&
                        & grd%mask(i0+1,j0+1,1) == 0.0 ) then
                      profiles(num_profiles)%accepted = .false.
                   end if
                end if
                if ( profiles(num_profiles)%accepted ) then
                   profiles(num_profiles)%flag(:) = .true.
                   allocate(profiles(num_profiles)%k_index(profiles(num_profiles)%levels))
                   do k=1, profiles(num_profiles)%levels
                      if (depth(k) < grd%z(1)) then
                        profiles(num_profiles)%k_index(k) = 0.0
                      else
                        profiles(num_profiles)%k_index(k) = frac_index(depth(k), (/0.,grd%z(:)/)) - 1.0 ! snz modify to v3.2 JAN3012
                      end if
                      if ( profiles(num_profiles)%k_index(k) < 1.0 ) profiles(num_profiles)%flag(k) = .false.
                   end do
                end if
             end if
          end if
       end do
    end do

    write (unit=stdout_unit, fmt='("A grand total of ",I8," suv points within global domain")') no_suv
    write (unit=stdout_unit, fmt='("A final total @suv of ",I8," prfs within global domain")') num_profiles
  end subroutine open_profile_dataset_suv

  subroutine get_obs_suv(model_time, Prof, nprof, no_prf0)
    type(time_type), intent(in) :: model_time
    type(ocean_profile_type), dimension(:), intent(inout) :: prof
    integer, intent(inout) :: nprof
    integer, intent(in) :: no_prf0

    ! get sst data and put into profiles
    ! only current day

    real :: sfc_lon, sfc_lat, ri0, rj0
    real, dimension(1440,1070,1) :: sfc_u, sfc_v

    integer :: i, k, i_m, i0, j0
    integer :: unit, time_idx, ndim, nvar, natt, ntime
    integer :: iy0, in0, id0, ih0, im0, is0
    integer :: stdout_unit, istat
    integer, dimension(12) :: n_days

    character(len=80) :: sfc_filename
    character(len=256) :: emsg_local

    type(fieldtype), dimension(:), allocatable :: fields
    type(time_type) :: tdiff, sfc_time0

    n_days  = (/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
    nprof = 0
    stdout_unit = stdout()

    call get_date(model_time, iy0, in0, id0, ih0, im0, is0)

    !monthly
!!$    time_idx = (iy0-1984)*12+in0
    ! daily
    if ( in0 == 1 ) then
       time_idx = (iy0-1984)*365 + id0
    else
       time_idx = 0
       do i_m=1, in0-1
          time_idx = time_idx + (iy0-1984)*365 + n_days(i_m)
       end do
       time_idx = time_idx + id0
    end if

    if ( mpp_pe() == mpp_root_pe() ) then
       write (unit=stdout_unit, fmt='("time_idx = ",I8)') time_idx
    end if

    sfc_time0 = set_date(iy0, in0, id0, ih0, im0, is0)

    sfc_filename = "INPUT/sfc_current.198401-198412.nc"

    call mpp_open(unit, trim(sfc_filename), mpp_rdonly, mpp_netcdf, threading=mpp_multi, fileset=mpp_single)
    call mpp_get_info(unit, ndim, nvar, natt, ntime)

    allocate(fields(nvar), stat=istat)
    if ( istat .ne. 0 ) then
       call error_mesg('oda_core_mod::get_obs_sst', 'Unable to allocate fields', fatal)
    end if

    call mpp_get_fields(unit, fields)
    do i=1, nvar
       select case ( mpp_get_field_name(fields(i)) )
       case ('U_SFC')
          call mpp_read(unit, fields(i), sfc_u, tindex=time_idx)
       case ('V_SFC')
          call mpp_read(unit, fields(i), sfc_v, tindex=time_idx)
       end select
    end do

    do i=no_prf+no_sst+no_eta+1, no_prf+no_sst+no_eta+no_suv
       profiles(i)%time = sfc_time0

       tdiff = model_time - profiles(i)%time

       sfc_lon = profiles(i)%lon
       ri0 = frac_index(sfc_lon, x_grid_uv(:,1))
       i0 = floor(ri0)
       if ( i0 < 1 .or. i0 > 1440 ) then
          write (unit=emsg_local, fmt='("Profiles(",I8,")%lon = ",I4,", i0 = ",I8)') i, profiles(i)%lon, i0
          call error_mesg('oda_core_mod::get_obs_suv',&
               & 'Profile longitude index outside range [1,1440].  '//trim(emsg_local), fatal)
       end if

       sfc_lat = profiles(i)%lat
       rj0 = frac_index(sfc_lat, y_grid_uv(90,:))
       j0 = floor(rj0)
       if ( j0 < 1 .or. j0 > 1070 ) then
          write (unit=emsg_local, fmt='("Profiles(",I8,")%lat = ",I4,", j0 = ",I8)') i, profiles(i)%lon, j0
          call error_mesg('oda_core_mod::get_obs_suv',&
               & 'Profile latitude index outside range [1,1070].  '//trim(emsg_local), fatal)
       end if

       if ( profiles(i)%accepted ) then
          nprof = nprof + 1
          if ( nprof > size(prof,1) ) then
             call error_mesg('oda_core_mod::get_obs_suv',&
                  & 'Passed in array "Prof" is smaller than number of profiles, increase size of Prof before call',&
                  & fatal)
          end if
          profiles(i)%data(1) = sfc_u(i0,j0,1)
          profiles(i)%data(2) = sfc_v(i0,j0,1)

          call copy_obs(profiles(i:i),prof(nprof+no_prf0:nprof+no_prf0))

          prof(nprof+no_prf0)%tdiff = tdiff
       end if
    end do

    deallocate(fields)
    call mpp_close(unit)
  end subroutine get_obs_suv

  subroutine get_obs_eta(model_time, Prof, nprof, no_prf0)
    type(time_type), intent(in) :: model_time
    type(ocean_profile_type), dimension(:), intent(inout) :: prof
    integer, intent(inout) :: nprof
    integer, intent(in) :: no_prf0

    ! get sst data and put into profiles
    ! only current day

    real :: eta_lon, eta_lat, ri0, rj0
    real, dimension(1440,1070) :: eta_t

    integer :: i, i0, j0, i_m
    integer :: iy0, in0, id0, ih0, im0, is0
    integer :: unit, time_idx, ndim, nvar, natt, ntime
    integer :: stdout_unit, istat
    integer, dimension(12) :: n_days

    character(len=80) :: eta_filename
    character(len=256) :: emsg_local

    type(fieldtype), dimension(:), allocatable :: fields
    type(time_type) :: tdiff, sfc_time0

    n_days  = (/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
    nprof = 0
    stdout_unit = stdout()

    call get_date(model_time, iy0, in0, id0, ih0, im0, is0)

    !monthly
!!$    time_idx = (iy0-1984)*12+in0
    ! daily
    if ( in0 == 1 ) then
       time_idx = (iy0-1976)*365 + id0 - 1
    else
       time_idx = 0
       do i_m=1, in0-1
          time_idx = time_idx + n_days(i_m)
       end do
       time_idx = (iy0-1976)*365 + time_idx + id0 - 1
    end if

    if ( mpp_pe() == mpp_root_pe() ) then
       write (unit=stdout_unit, fmt='("time_idx = ",I8)') time_idx
    end if

    sfc_time0 = set_date(iy0, in0, id0, ih0, im0, is0)

    eta_filename='INPUT/ocean.19760101-20001231.eta_t.nc'

    call mpp_open(unit, trim(eta_filename), mpp_rdonly, mpp_netcdf, threading=mpp_multi, fileset=mpp_single)
    call mpp_get_info(unit, ndim, nvar, natt, ntime)

    allocate(fields(nvar), stat=istat)
    if ( istat .ne. 0 ) then
       call error_mesg('oda_core_mod::get_obs_sst', 'Unable to allocate fields', fatal)
    end if

    call mpp_get_fields(unit, fields)
    do i=1, nvar
       select case ( mpp_get_field_name(fields(i)) )
       case ('eta_t')
          call mpp_read(unit, fields(i), eta_t, tindex=time_idx)
       end select
    end do

    do i=no_prf+no_sst+1, no_prf+no_sst+no_eta
       profiles(i)%time = sfc_time0

       tdiff = model_time - profiles(i)%time

       eta_lon = profiles(i)%lon
       ri0 = frac_index(eta_lon, x_grid(:,1))
       i0 = floor(ri0)
       if ( i0 < 1 .or. i0 > 1440 ) then
          write (unit=emsg_local, fmt='("Profiles(",I8,")%lon = ",I4,", i0 = ",I8)') i, profiles(i)%lon, i0
          call error_mesg('oda_core_mod::get_obs_eta',&
               & 'Profile longitude index outside range [1,1440].  '//trim(emsg_local), fatal)
       end if

       eta_lat = profiles(i)%lat
       rj0 = frac_index(eta_lat, y_grid(90,:))
       j0 = floor(rj0)
       if ( j0 < 1 .or. j0 > 1070 ) then
          write (unit=emsg_local, fmt='("Profiles(",I8,")%lat = ",I4,", j0 = ",I8)') i, profiles(i)%lon, j0
          call error_mesg('oda_core_mod::get_obs_suv',&
               & 'Profile latitude index outside range [1,1070].  '//trim(emsg_local), fatal)
       end if

       if ( profiles(i)%accepted ) then
          if ( eta_t(i0,j0) > -9.9 ) then !!! excluding missing values
             nprof = nprof + 1
             if ( nprof > size(prof,1) ) then
                call error_mesg('oda_core_mod::get_obs_eta',&
                     & 'Passed in array "Prof" is smaller than number of profiles, increase size of Prof before call.',&
                     & fatal)
             end if
             profiles(i)%data(1) = eta_t(i0,j0)

             call copy_obs(profiles(i:i),prof(nprof+no_prf0:nprof+no_prf0))

             prof(nprof+no_prf0)%tdiff = tdiff
          end if
       end if
    end do

    deallocate(fields)
    call mpp_close(unit)
  end subroutine get_obs_eta
end module oda_core_ecda_mod