!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! 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 General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module oda_core_mod
!
!<CONTACT EMAIL="matthew.harrison@noaa.gov"> Matthew Harrison
!</CONTACT>
!
!<OVERVIEW>
! core ocean data assimilation subroutines for ocean models. This module
! includes interfaces to :
!  
! (i)   initialize the ODA core with observations and model
! (ii)  request ocean state increments using observed or idealized data 
! (iii) calculate model guess differences and output to file
! (iv)  terminate the DA core.  
!  
! NOTE: Profile files conform to v0.1a profile metadata standard.
! This metadata standard will evolve in time and will maintain
! backward compability.  
!  
! NOTE: This code is still under development. ODA data structures should be 
! opaque in future releases.
!  
!</OVERVIEW>
!
!<NAMELIST NAME="oda_core_nml">
! <DATA NAME="max_misfit" TYPE="real">
!  flag measurement if abs(omf) exceeds max_misfit
!</DATA>  
! <DATA NAME="min_obs_err_t" TYPE="real">
! minumum rms temperature observation error (degC) for
! variable obs error, else nominal temp error  
!</DATA>    
! <DATA NAME="min_obs_err_s" TYPE="real">
! minimum rms salinity observation error (g/kg) for
! variable obs error,else nominal salt error  
!</DATA>    
! <DATA NAME="eta_tide_const" TYPE="real">
! Tidal internal displacement amplitude (meters) for observational error estimation.  
!</DATA>    
! <DATA NAME="min_prof_depth" TYPE="real">
! Minimum profile depth (meters)  
!</DATA>    
! <DATA NAME="max_prof_spacing" TYPE="real">
! Data must be contiguous to this resolution (meters). Otherwise flag is activated.
! </DATA>
! <DATA NAME="data_window" TYPE="integer">
! Half-width of profile time window (days)
! </DATA>
! <DATA NAME="add_tidal_aliasing" TYPE="logical">
! Add tidal aliasing to observation error. Use eta_tide_const * dT/dz to estimate
! observation.  
! </DATA>
! <DATA NAME="max_profiles" TYPE="integer">
! Allocation size of profile array
! </DATA>
! <DATA NAME="max_sfc_obs" TYPE="integer">
! Allocation size of surface data array
! </DATA>
! <DATA NAME="temp_obs_rmse" TYPE="real">
! nominal temperature error rms error
! </DATA>
! <DATA NAME="salt_obs_rmse" TYPE="real">
! nominal salinity error rms error
! </DATA>      
!</NAMELIST>  
  use fms_mod, only : file_exist,read_data
  use mpp_mod, only : mpp_error, FATAL, NOTE, mpp_sum, stdout,&
                      mpp_sync_self, mpp_pe,mpp_npes,mpp_root_pe,&
                      mpp_broadcast
  use mpp_domains_mod, only : mpp_get_compute_domain, mpp_get_data_domain, &
       domain2d, mpp_get_global_domain, mpp_update_domains
  use time_manager_mod, only : time_type, operator( <= ), operator( - ), &
       operator( > ), operator ( < ),  set_time, set_date, &
       get_date, get_time
  use get_cal_time_mod, only : get_cal_time
  use axis_utils_mod, only : frac_index  
  use constants_mod, only : radian, pi
  use oda_types_mod
  use write_ocean_data_mod, only : write_ocean_data_init
  
  implicit none

  private
  
  real, private :: max_misfit = 5.0 ! reject fg diff gt max_misfit


  real, parameter, private :: depth_min=0.0, depth_max=10000.
  real, parameter, private :: temp_min=-3.0, temp_max=40.
  real, parameter, private :: salt_min=0.0, salt_max=45.
  
  integer :: max_profiles = 250000, max_sfc_obs = 1 
  integer, parameter, private :: max_files=100 
  integer, parameter, private :: PROFILE_FILE = 1,SFC_FILE= 2,&
                                 IDEALIZED_PROFILES=3
  
  
! parameters for obs error inflation due to internal tidal displacements
  
  real, private :: min_obs_err_t = 0.5, min_obs_err_s=0.1, eta_tide_const = 7.0
  
  type(ocean_profile_type), target, save, private, allocatable  :: profiles(:)
  type(ocean_surface_type), target, save, private, allocatable  :: sfc_obs(:)  

  integer, private, save :: num_profiles, num_sfc_obs ! total number of observations 

  integer, private, save :: isc, iec, jsc, jec, isd, ied, jsd, jed  ! indices for local and global domain
  integer, private, save :: isg, ieg, jsg, jeg
  integer, private, save :: nk

! parameters used to restrict profiles
  real, private, save :: min_prof_depth = 200.0 ! profile data ending
                                                   ! above this level are
                                                   ! discarded.  
  real, private, save :: max_prof_spacing = 1.e5 ! reject profile data
                                                    ! if not contiguous
                                                    ! at this resolution

! internal arrays for forward and backward observations
  real, dimension(:,:,:), allocatable, private, save :: sum_wgt, nobs

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

! the DA module maintains a unique grid,domain association
  
  type(grid_type), pointer :: Grd

! DA grid (1-d) coordinates.  Generalized horizontal coordinates are not supporte.
  
  real, allocatable, dimension(:) :: x_grid, y_grid

  real :: temp_obs_rmse = 0.7071
  real :: salt_obs_rmse = 0.1
  
  logical :: add_tidal_aliasing=.false.

  logical :: localize_data = .true.

  logical :: debug=.false.

  integer :: ndebug=10

  integer, allocatable, dimension(:) :: nprof_in_comp_domain
  
  namelist /oda_core_nml/ max_misfit,add_tidal_aliasing,min_obs_err_t,&
                          min_obs_err_s, eta_tide_const, debug, max_profiles,&
                          max_sfc_obs, temp_obs_rmse, salt_obs_rmse, ndebug
                          

! NOTE: Surface observation type is not yet implemented.

! transform from Grd => Obs
  
  interface forward_obs
     module procedure forward_obs_profile
     module procedure forward_obs_sfc
  end interface

! transform from Obs => Grd
  
  interface backward_obs
     module procedure backward_obs_profile
     module procedure backward_obs_sfc
  end interface


! one-time association between profile and grid, used for storing
! interpolation weights
  
  interface assign_forward_model
     module procedure assign_forward_model_profile
     module procedure assign_forward_model_sfc
  end interface

! duplicate observation array
  
  interface copy_obs
     module procedure copy_obs_prof
     module procedure copy_obs_sfc
  end interface

! multiply observation data by inverse error variance
 
  interface mult_obs_I_mse
     module procedure mult_obs_I_mse_profile
     module procedure mult_obs_I_mse_sfc
  end interface

! difference between two observations (e.g. model first guess and observations)
  
  interface diff_obs
     module procedure diff_obs_profile
     module procedure diff_obs_sfc
  end interface

! inflate observational error based on first guess misfit 
! and time window.
  
  interface adjust_obs_error
     module procedure adjust_obs_error_profile
     module procedure adjust_obs_error_sfc
  end interface


  interface nullify_obs
     module procedure nullify_obs_prof
  end interface
  
  public :: forward_obs, backward_obs, &
            copy_obs, adjust_obs_error, &
            oda_core_init, open_profile_dataset, get_obs, &
            assign_forward_model, diff_obs, mult_obs_I_mse, &
            purge_obs, nullify_obs
  
  contains


  subroutine open_profile_dataset(filename, localize)  
!    <DESCRIPTION>
!    open dataset containing profile information in fms station format.
!    store internally.
!    </DESCRIPTION>    
    
    use mpp_io_mod, only : mpp_open, mpp_get_atts, mpp_get_info, &
         mpp_get_fields, mpp_read, MPP_SINGLE, MPP_MULTI, MPP_NETCDF,&
         axistype, atttype, fieldtype, MPP_RDONLY, mpp_get_axes, mpp_close,&
         mpp_get_att_char
    
    character(len=*), intent(in) :: filename
    logical, intent(in), optional :: localize
    logical :: found_neighbor,continue_looking
    integer, parameter :: max_levels=1000
    integer :: unit, ndim, nvar, natt, nstation
    character(len=32) :: fldname, axisname
    type(atttype), allocatable, dimension(:), target :: global_atts
    type(atttype), pointer :: version => NULL()
    type(axistype), pointer :: depth_axis => NULL(), station_axis => NULL()
    type(axistype), allocatable, dimension(:), target :: axes
    type(fieldtype), allocatable, dimension(:), target :: fields
    
    type(fieldtype), pointer :: field_lon => NULL(), field_lat => NULL(), field_probe => NULL(),&
         field_time => NULL(), field_depth => NULL()
    type(fieldtype), pointer :: field_temp => NULL(), field_salt => NULL(), &
                                field_error => NULL(), field_link => NULL()
    ! NOTE: fields are restricted to be in separate files
    real :: lon, lat, time, depth(max_levels), temp(max_levels), salt(max_levels),&
         error(max_levels), rprobe, profile_error, rlink
    integer :: nlev, probe, yr, mon, day, hr, minu, sec, kl, outunit
    integer :: num_levs, num_levs_temp, num_levs_salt,&
               k, kk, ll, i, i0, j0, k0, nlevs, a, nn, ii, nlinks
    real :: ri0, rj0, rk0, dx1, dx2, dy1, dy2
    character(len=128) :: time_units, attname, catt
    type(time_type) :: profile_time
    integer :: flag_t(max_levels), flag_s(max_levels), cpe
    logical :: data_is_local, &
               continue
    logical :: found_temp=.false., found_salt=.false.
    real, dimension(max_links,max_levels) :: temp_bfr, salt_bfr, depth_bfr
    integer, dimension(max_links,max_levels) :: flag_t_bfr, flag_s_bfr
    real :: temp_missing=missing_value,salt_missing=missing_value,&
         depth_missing=missing_value
    real :: max_prof_depth, zdist
    
    
    ! read timeseries of local observational data from NetCDF files
    ! and allocate ocean_obs arrays.  

    ! File structure:
    !    dimensions:
    !       depth_index;station=UNLIMITED;
    !    variables:
    !       depth_index(depth_index);
    !       station(station);
    !       longitude(station);
    !       latitude(station);
    !       time(station); 
    !       data(station,depth_index);
    !       depth(station,depth_index);
    !       probe(station);
    !       err(station, depth_index);

    cpe = mpp_pe()

    dx1 = (x_grid(isc)-x_grid(isc-1))/2.0
    dx2 = (x_grid(iec+1)-x_grid(iec))/2.0
    dy1 = (y_grid(jsc)-y_grid(jsc-1))/2.0
    dy2 = (y_grid(jec+1)-y_grid(jec))/2.0
    
    localize_data = .true.

    if (PRESENT(localize)) localize_data = localize

    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)

    outunit = stdout()
    write(outunit,*) 'Opened profile dataset :',trim(filename)

    ! get version number of profiles

    allocate(global_atts(natt))
    call mpp_get_atts(unit,global_atts)

    do i=1,natt
       catt = mpp_get_att_char(global_atts(i))
       select case (lowercase(trim(catt)))
       case ('version')
          version =>  global_atts(i)
       end select
    end do

    if (.NOT.ASSOCIATED(version)) then
        call mpp_error(NOTE,'no version number available for profile file, assuming v0.1a ')
    else
        write(outunit,*) 'Reading profile dataset version = ',trim(catt)
    endif
    
    
    ! 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 (lowercase(trim(axisname)))
       case ('depth_index')
          depth_axis => axes(i)
       case ('station_index')
          station_axis => axes(i)
       end select
    end do

    if (.NOT.ASSOCIATED(depth_axis) .or. .NOT.ASSOCIATED(station_axis)) then
        call mpp_error(FATAL,'depth and/or station axes do not exist in input file')
    endif

    
! get selected field information.
! NOTE: not checking for all variables here.    

    allocate(fields(nvar))
    call mpp_get_fields(unit,fields)
    do i=1,nvar
       call mpp_get_atts(fields(i),name=fldname)
       select case (lowercase(trim(fldname)))
       case ('longitude')
           field_lon => fields(i)
       case ('latitude')
           field_lat => fields(i)
       case ('probe') 
           field_probe => fields(i)
       case ('time')
           field_time => fields(i)
       case ('temp')
           field_temp => fields(i)
       case ('salt')
           field_salt => fields(i)           
       case ('depth')
           field_depth => fields(i)
       case ('link')
           field_link => fields(i)
       case ('error')
           field_error => fields(i)
       end select
   enddo

    call mpp_get_atts(depth_axis,len=nlevs)

    if (nlevs > max_levels) call mpp_error(FATAL,'increase parameter max_levels ')

    if (nlevs < 1) call mpp_error(FATAL)

    outunit = stdout()
    write(outunit,*) 'There are ', nstation, ' records in this dataset'
    write(outunit,*) 'Searching for profiles matching selection criteria ...'

    if (ASSOCIATED(field_temp)) found_temp=.true.
    if (ASSOCIATED(field_salt)) found_salt=.true.

    if (.not. found_temp .and. .not. found_salt) then
        write(outunit,*) 'temp or salt not found in profile file'
        call mpp_error(FATAL)
    endif
    
    call mpp_get_atts(field_time,units=time_units)
    if (found_temp) call mpp_get_atts(field_temp,missing=temp_missing)
    if (found_salt) call mpp_get_atts(field_salt,missing=salt_missing)
    
    call mpp_get_atts(field_depth,missing=depth_missing)        

    if (found_salt) then
        write(outunit,*) 'temperature and salinity where available'
    else
        write(outunit,*) 'temperature only records'
    endif

    
    i=1
    continue=.true.
    
    do while (continue)

       depth(:) = missing_value
       temp(:)  = missing_value
       salt(:)  = missing_value       
! required fields       
       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)
       if (found_temp) call mpp_read(unit,field_temp,temp(1:nlevs),tindex=i)
       if (found_salt) call mpp_read(unit,field_salt,salt(1:nlevs),tindex=i)       
! not required fields
       if (ASSOCIATED(field_error)) then
           call mpp_read(unit,field_error,profile_error,tindex=i)
       endif
       if (ASSOCIATED(field_probe)) then
           call mpp_read(unit, field_probe, rprobe,tindex=i)
       endif
       if (ASSOCIATED(field_link)) then
           call mpp_read(unit,field_link,rlink,tindex=i)
       else
           rlink = 0.0
       endif
       probe=rprobe
       data_is_local = .false.
! NOTE: assuming grid is modulo 360 here. This needs to be generalized.
       
       if (lon .lt. x_grid(isg-1) ) lon = lon + 360.0
       if (lon .gt. x_grid(ieg+1) ) lon = lon - 360.0

! localized data is within region bounded by halo points
! (halo size = 1) adjacent to boundary points of computational domain
       
       if (lon >= x_grid(isc-1) .and. &
           lon <  x_grid(iec+1) .and. &
           lat >= y_grid(jsc-1) .and. &
           lat <  y_grid(jec+1)) data_is_local = .true.

       
       profile_time = get_cal_time(time,time_units,'julian')

       
       if ( data_is_local .OR. .NOT.localize_data) then
           
           num_profiles=num_profiles+1
           if (num_profiles > max_profiles) then
               call mpp_error(FATAL,'maximum number of profiles exceeded.&
                    &Resize parameter max_profiles in ocean_obs_mod')
               
           endif

           call nullify_obs(Profiles(num_profiles))
           
           num_levs_temp = 0
           num_levs_salt = 0           
           do k = 1, nlevs
              
! flag=0 denotes a valid profile level, anything else
! is invalid. See NODC codes.
!================================================================
!0 -     accepted station
!1 -     failed annual standard deviation check
!2 -     two or more density inversions (Levitus, 1982 criteria)
!3 -     flagged cruise
!4 -     failed seasonal standard deviation check
!5 -     failed monthly standard deviation check
!6 -     flag 1 and flag 4
!7 -     bullseye from standard level data or failed annual and monthly
!        standard deviation check
!8 -     failed seasonal and monthly standard deviation check
!9 -     failed annual, seasonal, and monthly standard deviation check
!================================================================
              
              flag_t(k) = 0;flag_s(k) = 0
              
              if (.not.found_salt) then
                  flag_s(k) = 1
              endif
              
              if (depth(k) .eq. depth_missing .or. depth(k) .lt.depth_min&
                   .or. depth(k) .gt. depth_max) then
                  depth(k) = missing_value
                  flag_t(k)=1
                  flag_s(k)=1
              endif
              
              if (found_temp .and. flag_t(k) .eq. 0) then
                  if (temp(k) .eq. temp_missing .or. temp(k) .lt. temp_min&
                       .or. temp(k) .gt. temp_max) then
                      temp(k) = missing_value
                      flag_t(k) = 1
                      flag_s(k) = 1 ! flag salt if no temperature data
                  else 
                      num_levs_temp = num_levs_temp+1
                  endif
              endif
              if (found_salt .and. flag_s(k) .eq. 0) then
                  if (salt(k) .eq. salt_missing .or. salt(k) .lt. salt_min&
                       .or. salt(k) .gt. salt_max) then
                      salt(k) = missing_value
                      flag_s(k) = 1
                  else 
                      num_levs_salt = num_levs_salt+1
                  endif
              endif
              
           enddo

! large profile are stored externally in separate records
! follow the links to get complete profile
           
           ii=i+1
           nlinks = 0
           do while (rlink > 0.0 .and. nlinks .lt. max_links)
              nlinks=nlinks+1
              depth_bfr(nlinks,:) = missing_value
              temp_bfr(nlinks,:) = missing_value
              salt_bfr(nlinks,:) = missing_value              
              call mpp_read(unit,field_depth,depth_bfr(nlinks,1:nlevs),tindex=ii)
              if (found_temp) call mpp_read(unit,field_temp,temp_bfr(nlinks,1:nlevs),tindex=ii)
              if (found_salt) call mpp_read(unit,field_salt,salt_bfr(nlinks,1:nlevs),tindex=ii)
              call mpp_read(unit,field_link,rlink,tindex=ii)
              ii=ii+1
           enddo
           i=ii ! set record counter to start of next profile

           if (nlinks > 0) then
               do nn = 1, nlinks
                  do k=1, nlevs
                     flag_t_bfr(nn,k) = 0
                     flag_s_bfr(nn,k) = 0
                     if (depth_bfr(nn,k) .eq. depth_missing .or.&
                          depth_bfr(nn,k) .lt. depth_min .or. &
                          depth_bfr(nn,k) .gt. depth_max) then
                         depth_bfr(nn,k) = missing_value
                         flag_t_bfr(nn,k)  = 1
                         flag_s_bfr(nn,k)  = 1
                     endif
                     if (found_temp .and. flag_t_bfr(nn,k) .eq. 0) then
                         if (temp_bfr(nn,k) .eq. temp_missing .or.&
                              temp_bfr(nn,k) .lt. temp_min .or.&
                              temp_bfr(nn,k) .gt. temp_max) then
                             temp_bfr(nn,k) = missing_value
                             flag_t_bfr(nn,k) = 1
                             flag_s_bfr(nn,k) = 1                         
                         else 
                             num_levs_temp = num_levs_temp+1
                         endif
                     endif
                     if (found_salt .and. flag_s_bfr(nn,k) .eq. 0) then
                         if (salt_bfr(nn,k) .eq. salt_missing  .or.&
                              salt_bfr(nn,k) .lt. salt_min .or.&
                              salt_bfr(nn,k) .gt. salt_max) then
                             salt_bfr(nn,k) = missing_value
                             flag_t_bfr(nn,k) = 0                         
                             flag_s_bfr(nn,k) = 1
                         else
                             num_levs_salt = num_levs_salt+1
                         endif
                     endif
                  enddo
               enddo
           endif

           num_levs = max(num_levs_temp,num_levs_salt)
           
           if (num_levs == 0) then
               if (i .gt. nstation) continue = .false.
               cycle
           endif

           allocate(profiles(num_profiles)%depth(num_levs))
           profiles(num_profiles)%depth=missing_value
           if (num_levs_temp .gt. 0) then
               allocate(profiles(num_profiles)%data_t(num_levs))
               profiles(num_profiles)%data_t=missing_value
               allocate(profiles(num_profiles)%flag_t(num_levs))
               profiles(num_profiles)%flag_t= 1
               profiles(num_profiles)%nvar=1
           endif

           if (num_levs_salt .gt. 0) then
               allocate(profiles(num_profiles)%data_s(num_levs))
               profiles(num_profiles)%data_s=missing_value
               allocate(profiles(num_profiles)%flag_s(num_levs))
               profiles(num_profiles)%flag_s= 1
               profiles(num_profiles)%nvar=profiles(num_profiles)%nvar + 1
           endif
           

           if (probe < 1 )   probe = 0
           profiles(num_profiles)%probe = probe
           profiles(num_profiles)%levels = num_levs
           profiles(num_profiles)%lat = lat
           profiles(num_profiles)%lon = lon
           allocate(profiles(num_profiles)%ms_t(num_levs))
           profiles(num_profiles)%ms_t(:) = temp_obs_rmse**2.0 ! default error variance for temperature

           if(num_levs_salt .gt. 0) then
               allocate(profiles(num_profiles)%ms_s(num_levs))
               profiles(num_profiles)%ms_s(:) = salt_obs_rmse**2.0  ! default error variance for salinity
           endif
           
           kk= 1
           do k = 1, nlevs
              if (flag_t(k) .eq. 0) then
                  if (kk > profiles(num_profiles)%levels) then
                      call mpp_error(FATAL)
                  endif
                  profiles(num_profiles)%depth(kk) = depth(k)
                  profiles(num_profiles)%data_t(kk) = temp(k)
                  profiles(num_profiles)%flag_t(kk) = 0                  
                  if (found_salt .and. flag_s(k) .eq. 0) then
                      profiles(num_profiles)%data_s(kk) = salt(k)
                      profiles(num_profiles)%flag_s(kk) = 0
                  endif
                  kk=kk+1
              endif
           enddo

           do nn = 1, nlinks
              do k = 1, nlevs
                 if (flag_t_bfr(nn,k) .eq. 0) then
                     if (kk > profiles(num_profiles)%levels) then
                         call mpp_error(FATAL)
                     endif
                     profiles(num_profiles)%depth(kk) = depth_bfr(nn,k)
                     profiles(num_profiles)%data_t(kk) = temp_bfr(nn,k)
                     profiles(num_profiles)%flag_t(kk) = 0
                     if (found_salt .and. flag_s_bfr(nn,k) .eq. 0) then
                         profiles(num_profiles)%data_s(kk) = salt_bfr(nn,k)
                         profiles(num_profiles)%flag_s(kk) = 0
                     endif
                     kk=kk+1
                 endif
              enddo
           enddo
           
           profiles(num_profiles)%time = profile_time
           
! calculate interpolation coefficients (make sure to account for grid offsets here!)
! NOTE: this only works for lat/lon grids. Lower left indices.           
!       
           
           ri0 = frac_index(lon, x_grid(isg-1:ieg+1)) - 1 
           rj0 = frac_index(lat, y_grid(jsg-1:jeg+1)) - 1
           i0 = floor(ri0)
           j0 = floor(rj0)
           Profiles(num_profiles)%i_index = ri0
           Profiles(num_profiles)%j_index = rj0           
           Profiles(num_profiles)%accepted = .true.
           if (i0 < 0 .or. j0 < 0) then
               Profiles(num_profiles)%accepted = .false.
           endif
           if (i0 > ieg .or. j0 > jeg) then  
               call mpp_error(FATAL,'grid lat/lon index is out of bounds ')
           endif
           if ((i0 < isc-1 .or. i0 > iec) .and. localize_data) then
               call mpp_error(FATAL,'grid lat/lon index is out of bounds ')
           endif
           if ((j0 < jsc-1 .or. j0 > jec) .and. localize_data) then
               call mpp_error(FATAL,'grid lat/lon index is out of bounds ')
           endif
!
! flag the profile if it sits on a model land point
!           
           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.
               endif
           endif

           if (Profiles(num_profiles)%accepted) then
               allocate(Profiles(num_profiles)%k_index(Profiles(num_profiles)%levels))
               max_prof_depth=0.0
               do k=1, Profiles(num_profiles)%levels
                  k0=0
                  if (Profiles(num_profiles)%flag_t(k).eq.0) then
                      rk0 = frac_index(Profiles(num_profiles)%depth(k), Grd%z(:))
                      k0 = floor(rk0)
                      if ( k0 == -1) then
                          if (Profiles(num_profiles)%depth(k) .lt. Grd%z(1)) then
                              k0 = 1
                              rk0 = 1.0
                          else if (Profiles(num_profiles)%depth(k) .gt. Grd%z(Grd%nk)) then
                              Profiles(num_profiles)%flag_t(k) = 1
                          endif
                      endif
                  else
                      cycle
                  endif

                  if (k0 .gt. size(Grd%z)-1 ) then
                      write(*,*) 'k0 out of bounds, rk0,k0= ',rk0,k0
                     write(*,*) 'Z_bound= ',Grd%z_bound
                     write(*,*) 'Profile%depth= ',Profiles(num_profiles)%depth
                     
                      call mpp_error(FATAL)
                  endif
                  
                  Profiles(num_profiles)%k_index(k) = rk0

! flag depth level if adjacent model grid points are land

                  if (Profiles(num_profiles)%flag_t(k) .eq. 0) then
                      if (i0 .lt. 0 .or. j0 .lt. 0 .or. k0 .lt. 0) then
                          write(*,*) 'profile index out of bounds'
                          write(*,*) 'i0,j0,k0=',i0,j0,k0
                          write(*,*) 'lon,lat,depth=',Profiles(num_profiles)%lon,&
                               Profiles(num_profiles)%lat,Profiles(num_profiles)%depth(k)
                          call mpp_error(FATAL)
                      endif
                      
                      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_t(k) = 1
                      endif
                      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_t(k) = 1
                      endif
                      if (Profiles(num_profiles)%flag_t(k) .eq. 0) then
                          max_prof_depth = Profiles(num_profiles)%depth(k)
                      endif
                  endif

               enddo ! Prof%levels loop

! Flag profile if it is too shallow.
               
               if (max_prof_depth .lt. min_prof_depth) then
                   Profiles(num_profiles)%accepted = .false.
               endif

               found_neighbor=.false.
               
               do k=2,Profiles(num_profiles)%levels - 1
                  if (Profiles(num_profiles)%flag_t(k) .eq. 0) then
                      kk = k-1
                      found_neighbor = .false.
                      continue_looking = .true.
                      do while (continue_looking .and. kk .ge. 1)
                         if (Profiles(num_profiles)%flag_t(kk) .eq. 0) then
                             zdist = Profiles(num_profiles)%depth(k) - Profiles(num_profiles)%depth(kk)
                             if (zdist .gt. max_prof_spacing) then
                                 Profiles(num_profiles)%accepted = .false.
                                 goto 199
                             else 
                                 continue_looking = .false.
                                 found_neighbor = .true.
                             endif
                         else
                             kk = kk - 1
                         endif
                      end do
                      kk = k+1
                      continue_looking = .true.
                      do while (continue_looking .and. kk .le. Profiles(num_profiles)%levels)
                         if (Profiles(num_profiles)%flag_t(kk).eq. 0) then
                             zdist = Profiles(num_profiles)%depth(kk) - Profiles(num_profiles)%depth(k)
                             if (zdist .gt. max_prof_spacing) then
                                 Profiles(num_profiles)%accepted = .false.
                                 goto 199
                             else
                                 continue_looking = .false.
                                 found_neighbor = .true.
                             endif
                         else
                             kk = kk+1
                         endif
                      enddo
                  endif
               enddo

               if (.not. found_neighbor) Profiles(num_profiles)%accepted = .false.

199            continue
               
           endif ! if Prof%accept
       else ! data is not local
           i = i+1
       endif ! if data_is_local

       
       if (i .gt. nstation) continue = .false.

    enddo

!    a = nprof_in_comp_domain(cpe)


    
!    call mpp_broadcast(nprof_in_comp_domain(cpe),cpe)
    
!    call mpp_sum(a)

!    write(stdout(),*) 'A grand total of ',int(a),' profiles satisify acceptance criteria'

!    do i=0,mpp_npes()
!       write(stdout(),*) 'pe=',i,'profile count=',nprof_in_comp_domain(i)
!    enddo
    
    call mpp_sync_self()
    call mpp_close(unit)

  end subroutine open_profile_dataset

  subroutine get_obs(model_time, Prof, Sfc, nprof, nsfc)


    ! get profiles and sfc
    ! obs relevant to current analysis interval

    type(time_type), intent(in) :: model_time
    type(ocean_profile_type), dimension(:) :: Prof
    type(ocean_surface_type), dimension(:) :: Sfc
    integer, intent(inout) :: nprof, nsfc

    integer :: i,k,yr,mon,day,hr,minu,sec,a,mon_obs,yr_obs, outunit
    type(time_type) :: tdiff
    character(len=1) :: cchar

    nprof=0
    nsfc=0

    outunit = stdout()
    write(outunit,*) 'Gathering profiles for current analysis time'
    call get_date(model_time,yr,mon,day,hr,minu,sec)
    write(outunit,'(a,i4,a,i2,a,i2)') 'Current yyyy/mm/dd= ',yr,'/',mon,'/',day
    write(outunit,*) 'num_profiles=',num_profiles

    
    do i=1,num_profiles

       if (debug .and. i.le.ndebug) then
           call get_date(Profiles(i)%time,yr,mon,day,hr,minu,sec)
           write(*,*) 'in get_obs prof time: yy/mm/dd= ',yr,mon,day
       endif
       
       if (Profiles(i)%time <= model_time) then
           tdiff = model_time - Profiles(i)%time
       else
           tdiff = Profiles(i)%time - model_time
       endif

       if (debug .and. i .le. ndebug) then
           write(*,*) 'Prof%accepted=',Profiles(i)%accepted
       endif
       
       if (tdiff <= time_window(0) .and. &
            Profiles(i)%accepted) then
           nprof=nprof+1
           if (nprof > size(Prof,1)) &
                call mpp_error(FATAL,'increase size of Prof array before call to get_obs')
           call copy_obs(Profiles(i:i),Prof(nprof:nprof))
           Prof(nprof)%tdiff = tdiff
           if (debug .and. nprof .le. ndebug) then
              call get_time(tdiff,sec,day)
              write(*,'(a,i3,a,2f10.5)') 'Accepted profile #',i,' : lon,lat= ',Prof(nprof)%lon,Prof(nprof)%lat
              do k=1,Prof(nprof)%levels
                 if (Prof(nprof)%nvar .eq. 2) then
                     write(*,'(a,i3,a,2f10.5,2i2,2f8.5)') 'Profile #',i,' : temp,salt,flag_t,flag_s,ms_t,ms_s= ',&
                          Prof(nprof)%data_t(k),Prof(nprof)%data_s(k),Prof(nprof)%flag_t(k),Prof(nprof)%flag_s(k),&
                          Prof(nprof)%ms_t(k),Prof(nprof)%ms_s(k)
                 else
                     write(*,'(a,i3,a,2f10.5)') 'Profile #',i,' : temp,flag_t= ',Prof(nprof)%data_t(k),Prof(nprof)%flag_t(k)
                 endif
              enddo
          endif
          
      else
          if (debug .and. i .le. ndebug) then
              call get_time(tdiff,sec,day)
              write(*,'(a,i3,a,2f10.5)') 'Rejected profile #',i,' : lon,lat= ',Prof(i)%lon,Prof(i)%lat
              do k=1,Prof(i)%levels
                 if (Prof(i)%nvar .eq. 2) then
                     write(*,'(a,i3,a,2f10.5,2i2)') 'Profile #',i,' : temp,salt,flag_t,flag_s= ',Prof(i)%data_t(k),Prof(i)%data_s(k),Prof(i)%flag_t(k),Prof(i)%flag_s(k)
                 else
                     write(*,'(a,i3,a,2f10.5)') 'Profile #',i,' : temp,flag_t= ',Prof(i)%data_t(k),Prof(i)%flag_t(k)
                 endif
              enddo
          endif
      endif

    enddo

    a=nprof
    call mpp_sum(a)
    write(outunit,*) 'A total of ',a,'  profiles are being used for the current analysis step'

    return

  end subroutine get_obs

  subroutine oda_core_init(Domain, Grid, localize)

    use fms_mod, only : open_namelist_file, check_nml_error, close_file
    
    type(domain2d), intent(in) :: Domain
    type(grid_type), target, intent(in) :: Grid
    logical, intent(in), optional :: localize

      
    integer :: ioun, ierr, io_status
    
    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)      


!    allocate(nprof_in_comp_domain(0:mpp_npes()-1))
!    nprof_in_comp_domain = 0
    
    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)
    nk = size(Grid%z)

    allocate(sum_wgt(isd:ied,jsd:jed,nk))
    allocate(nobs(isd:ied,jsd:jed,nk))

    if (PRESENT(localize)) localize_data = localize
    
    call init_observations(localize_data)

    call write_ocean_data_init()
    
  end subroutine oda_core_init

  subroutine purge_obs()

    num_profiles=0
    num_sfc_obs=0

  end subroutine purge_obs

  subroutine forward_obs_profile(Model_obs, fg_t, fg_s)

! map first guess to observation locations
! note that forward operator only becomes associated after
! this call

    type(ocean_profile_type), dimension(:), intent(inout) ::Model_obs    
    real, dimension(isd:ied,jsd:jed,nk) , intent(in) :: fg_t ! first guess for temperature
    real, dimension(isd:ied,jsd:jed,nk) , intent(in), optional :: fg_s ! first guess for salinity

    integer :: n, i0, j0, k, num_prof, k0
    real :: a,b,c

    character(len=128) :: mesg
    
    num_prof = size(Model_obs)

    sum_wgt = 0.0

    do n = 1, num_prof
       i0 = floor(Model_obs(n)%i_index)
       j0 = floor(Model_obs(n)%j_index)
       a = Model_obs(n)%i_index - i0 
       b = Model_obs(n)%j_index - j0 

       if (a >= 1.0 .or. a < 0.0 ) call mpp_error(FATAL)
       if (b >= 1.0 .or. b < 0.0 ) call mpp_error(FATAL)


       if (i0 < isc-1 .or. i0 > iec) then
           write(mesg,'(a,i4,2x,i4)') 'out of bounds in forward_obs: i0,j0= ',i0,j0
           call mpp_error(FATAL,trim(mesg))
       endif

       if (j0 < jsc-1 .or. j0 > jec) then
           write(mesg,*) 'out of bounds in forward_obs: i0,j0= ',i0,j0
           call mpp_error(FATAL,trim(mesg))
       endif

       if (ASSOCIATED(Model_obs(n)%data_t) .and. Model_obs(n)%accepted) then

           if (ASSOCIATED(Model_obs(n)%Forward_model%wgt))&
                Model_obs(n)%Forward_model%wgt => NULL()
           allocate(Model_obs(n)%Forward_model%wgt(2,2,2,Model_obs(n)%levels))

           Model_obs(n)%Forward_model%wgt = 0.0
           
           do k = 1, Model_obs(n)%levels
              if (Model_obs(n)%flag_t(k) .eq. 0) then
                  k0 = floor(Model_obs(n)%k_index(k))
                  if (k0 < 1 .or. k0 > Grd%nk-1) then
                      write(mesg,*) 'out of bounds in forward_obs: k0= ',k0
                      call mpp_error(FATAL,trim(mesg))
                  endif          
                  c =	 Model_obs(n)%k_index(k) - k0

                  if (c >= 1.0 .or. c < 0.0 ) call mpp_error(FATAL)

                  Model_obs(n)%Forward_model%wgt(1,1,1,k) = (1.0-a)*(1.0-b)*(1.0-c)
                  Model_obs(n)%Forward_model%wgt(2,1,1,k) = a*(1.0-b)*(1.0-c)
                  Model_obs(n)%Forward_model%wgt(1,2,1,k) = (1.0-a)*b*(1.0-c)
                  Model_obs(n)%Forward_model%wgt(2,2,1,k) = a*b*(1.0-c)
                  Model_obs(n)%Forward_model%wgt(1,1,2,k) = (1.0-a)*(1.0-b)*c
                  Model_obs(n)%Forward_model%wgt(2,1,2,k) = a*(1.0-b)*c
                  Model_obs(n)%Forward_model%wgt(1,2,2,k) = (1.0-a)*b*c
                  Model_obs(n)%Forward_model%wgt(2,2,2,k) = a*b*c

                  sum_wgt(i0,j0,k0) = sum_wgt(i0,j0,k0)+&
                       Model_obs(n)%Forward_model%wgt(1,1,1,k)
                  sum_wgt(i0+1,j0,k0) = sum_wgt(i0+1,j0,k0)+&
                       Model_obs(n)%Forward_model%wgt(2,1,1,k)
                  sum_wgt(i0,j0+1,k0) = sum_wgt(i0,j0+1,k0)+&
                       Model_obs(n)%Forward_model%wgt(1,2,1,k)
                  sum_wgt(i0+1,j0+1,k0) = sum_wgt(i0+1,j0+1,k0)+&
                       Model_obs(n)%Forward_model%wgt(2,2,1,k)
                  sum_wgt(i0,j0,k0+1) = sum_wgt(i0,j0,k0+1)+&
                       Model_obs(n)%Forward_model%wgt(1,1,2,k)
                  sum_wgt(i0+1,j0,k0+1) = sum_wgt(i0+1,j0,k0+1)+&
                       Model_obs(n)%Forward_model%wgt(2,1,2,k)
                  sum_wgt(i0,j0+1,k0+1) = sum_wgt(i0,j0+1,k0+1)+&
                       Model_obs(n)%Forward_model%wgt(1,2,2,k)
                  sum_wgt(i0+1,j0+1,k0+1) = sum_wgt(i0+1,j0+1,k0+1)+&
                       Model_obs(n)%Forward_model%wgt(2,2,2,k)
                  
                      
                  Model_obs(n)%data_t(k) = &
                       fg_t(i0,j0,k0)*Model_obs(n)%Forward_model%wgt(1,1,1,k) &
                     + fg_t(i0+1,j0,k0)*Model_obs(n)%Forward_model%wgt(2,1,1,k) &
                     + fg_t(i0,j0+1,k0)*Model_obs(n)%Forward_model%wgt(1,2,1,k) &
                     + fg_t(i0+1,j0+1,k0)*Model_obs(n)%Forward_model%wgt(2,2,1,k) &
                     + fg_t(i0,j0,k0+1)*Model_obs(n)%Forward_model%wgt(1,1,2,k) &
                     + fg_t(i0+1,j0,k0+1)*Model_obs(n)%Forward_model%wgt(2,1,2,k) &
                     + fg_t(i0,j0+1,k0+1)*Model_obs(n)%Forward_model%wgt(1,2,2,k) &
                     + fg_t(i0+1,j0+1,k0+1)*Model_obs(n)%Forward_model%wgt(2,2,2,k)

                  if (ASSOCIATED(Model_obs(n)%data_s) .and. PRESENT(fg_s)) then
                      Model_obs(n)%data_s(k) = &
                       fg_s(i0,j0,k0)*Model_obs(n)%Forward_model%wgt(1,1,1,k) &
                     + fg_s(i0+1,j0,k0)*Model_obs(n)%Forward_model%wgt(2,1,1,k) &
                     + fg_s(i0,j0+1,k0)*Model_obs(n)%Forward_model%wgt(1,2,1,k) &
                     + fg_s(i0+1,j0+1,k0)*Model_obs(n)%Forward_model%wgt(2,2,1,k) &
                     + fg_s(i0,j0,k0+1)*Model_obs(n)%Forward_model%wgt(1,1,2,k) &
                     + fg_s(i0+1,j0,k0+1)*Model_obs(n)%Forward_model%wgt(2,1,2,k) &
                     + fg_s(i0,j0+1,k0+1)*Model_obs(n)%Forward_model%wgt(1,2,2,k) &
                     + fg_s(i0+1,j0+1,k0+1)*Model_obs(n)%Forward_model%wgt(2,2,2,k) 
                  endif
              else
                  if (ASSOCIATED(Model_obs(n)%data_t)) then
                      Model_obs(n)%data_t(k) = missing_value
                  endif
                  if (ASSOCIATED(Model_obs(n)%data_s)) then
                      Model_obs(n)%data_s(k) = missing_value
                  endif                  
              endif
           enddo
       endif
    enddo

    return

  end subroutine forward_obs_profile





  subroutine forward_obs_sfc(Sfc, Guess, Diff)
    type(ocean_surface_type), intent(in) :: Sfc
    type(ocean_dist_type), intent(in) :: Guess
    type(ocean_surface_type), intent(inout) :: Diff


    return

  end subroutine forward_obs_sfc

  subroutine backward_obs_profile(Obs,model_t,model_s)
!
! map observations back to model locations
!
    type(ocean_profile_type), dimension(:), intent(in) :: Obs
    real, dimension(isd:ied,jsd:jed,nk), intent(inout) :: model_t
    real, dimension(isd:ied,jsd:jed,nk), intent(inout), optional :: model_s    

    real :: tmp
    integer :: i0,j0,k0,k,n

    model_t = 0.0
    if (PRESENT(model_s)) model_s = 0.0
    
    do n=1,size(Obs)
       if (ASSOCIATED(Obs(n)%data_t) .and. Obs(n)%accepted) then
           i0 = floor(Obs(n)%i_index)
           j0 = floor(Obs(n)%j_index) 
! profiles are assumed to lie           
! in domain bounded by first halo points

           if (i0 < isd .or. i0 > ied-1) cycle
           if (j0 < jsd .or. j0 > jed-1) cycle
           
           if (.not.ASSOCIATED(Obs(n)%Forward_model%wgt)) call mpp_error(FATAL,'forward operator not associated with obs')
           
           do k=1, Obs(n)%levels
              if (Obs(n)%flag_t(k) .eq. 0) then
                  k0 = floor(Obs(n)%k_index(k))
                  if (k0 < 1 .or. k0 > Grd%nk-1) call mpp_error(FATAL,'profile k indx out of bnds')
                  
                  tmp = Obs(n)%data_t(k) 
                  model_t(i0,j0,k0)       = model_t(i0,j0,k0) + tmp*Obs(n)%Forward_model%wgt(1,1,1,k)
                  model_t(i0+1,j0,k0)     = model_t(i0+1,j0,k0) + tmp*Obs(n)%Forward_model%wgt(2,1,1,k)
                  model_t(i0,j0+1,k0)     = model_t(i0,j0+1,k0) + tmp*Obs(n)%Forward_model%wgt(1,2,1,k)
                  model_t(i0+1,j0+1,k0)   = model_t(i0+1,j0+1,k0) + tmp*Obs(n)%Forward_model%wgt(2,2,1,k)
                  model_t(i0,j0,k0+1)     = model_t(i0,j0,k0+1) + tmp*Obs(n)%Forward_model%wgt(1,1,2,k)
                  model_t(i0+1,j0,k0+1)   = model_t(i0+1,j0,k0+1) + tmp*Obs(n)%Forward_model%wgt(2,1,2,k)
                  model_t(i0,j0+1,k0+1)   = model_t(i0,j0+1,k0+1) + tmp*Obs(n)%Forward_model%wgt(1,2,2,k)
                  model_t(i0+1,j0+1,k0+1) = model_t(i0+1,j0+1,k0+1) + tmp*Obs(n)%Forward_model%wgt(2,2,2,k)

                  if (PRESENT(model_s)) then

                      tmp = Obs(n)%data_s(k) 
                      model_s(i0,j0,k0)       = model_s(i0,j0,k0) + tmp*Obs(n)%Forward_model%wgt(1,1,1,k)
                      model_s(i0+1,j0,k0)     = model_s(i0+1,j0,k0) + tmp*Obs(n)%Forward_model%wgt(2,1,1,k)
                      model_s(i0,j0+1,k0)     = model_s(i0,j0+1,k0) + tmp*Obs(n)%Forward_model%wgt(1,2,1,k)
                      model_s(i0+1,j0+1,k0)   = model_s(i0+1,j0+1,k0) + tmp*Obs(n)%Forward_model%wgt(2,2,1,k)
                      model_s(i0,j0,k0+1)     = model_s(i0,j0,k0+1) + tmp*Obs(n)%Forward_model%wgt(1,1,2,k)
                      model_s(i0+1,j0,k0+1)   = model_s(i0+1,j0,k0+1) + tmp*Obs(n)%Forward_model%wgt(2,1,2,k)
                      model_s(i0,j0+1,k0+1)   = model_s(i0,j0+1,k0+1) + tmp*Obs(n)%Forward_model%wgt(1,2,2,k)
                      model_s(i0+1,j0+1,k0+1) = model_s(i0+1,j0+1,k0+1) + tmp*Obs(n)%Forward_model%wgt(2,2,2,k)                      
                  endif

              end if

                  
           end do
       end if
    end do


    where(sum_wgt > 0)
        model_t = model_t /sum_wgt
    elsewhere
        model_t = 0.0
    end where

    if (PRESENT(model_s)) then
        where(sum_wgt > 0)
            model_s = model_s /sum_wgt
        elsewhere
            model_s = 0.0
        end where
    endif
    
  end subroutine backward_obs_profile


  subroutine backward_obs_sfc(Obs,model)

    type(ocean_surface_type), dimension(:), intent(in) :: Obs
    real, dimension(isd:ied,jsd:jed,nk), intent(inout) :: model


  end subroutine backward_obs_sfc
  
  subroutine assign_forward_model_profile(Obs1,Obs2)

    type(ocean_profile_type), dimension(:), target, intent(in) :: Obs1
    type(ocean_profile_type), dimension(:), intent(inout) :: Obs2    


    integer :: n

    if (size(Obs1) .ne. size(Obs2)) call mpp_error(FATAL)

    do n=1,size(Obs1)

       Obs2(n)%Forward_model%wgt => Obs1(n)%Forward_model%wgt

    enddo

  end subroutine assign_forward_model_profile


  subroutine assign_forward_model_sfc(Obs1,Obs2)

    type(ocean_surface_type), target, intent(in) :: Obs1
    type(ocean_surface_type), intent(inout) :: Obs2    
    

    return
  end subroutine assign_forward_model_sfc
  
  subroutine diff_obs_profile(prof1, prof2, Diff)

    type(ocean_profile_type), dimension(:), intent(in) :: prof1
    type(ocean_profile_type), dimension(:), intent(in) :: prof2    
    type(ocean_profile_type), dimension(:), intent(inout) :: Diff

    integer :: n,k

    if (size(prof1) .ne. size(prof2) ) call mpp_error(FATAL)
    
    if (size(prof1) .ne. size(Diff) ) call mpp_error(FATAL)    


    do n=1,size(prof1)
       do k=1,prof1(n)%levels
          if (prof1(n)%flag_t(k) .eq. 0) then
              Diff(n)%data_t(k) = prof2(n)%data_t(k)-prof1(n)%data_t(k)
          else
              Diff(n)%data_t(k) = missing_value
          endif
          if (abs(Diff(n)%data_t(k)) .gt. max_misfit) then
              Diff(n)%flag_t(k) = 1
          endif
         if (ASSOCIATED(prof1(n)%data_s)) then
             if (prof1(n)%flag_s(k) .eq. 0) then
                 Diff(n)%data_s(k) = prof2(n)%data_s(k)-prof1(n)%data_s(k)
             else
                 Diff(n)%data_s(k) = missing_value
             endif
             if (abs(Diff(n)%data_s(k)) .gt. max_misfit) then
                 Diff(n)%flag_s(k) = 1
             endif
         endif
       enddo
    enddo


  end subroutine diff_obs_profile

  subroutine diff_obs_sfc(prof1,prof2,Diff)

    type(ocean_surface_type), dimension(:), intent(in) :: prof1, prof2
    type(ocean_surface_type), dimension(:), intent(inout) :: Diff


  end subroutine diff_obs_sfc
  
  subroutine copy_obs_prof(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)) call mpp_error(FATAL,&
         'size mismatch in call to copy_obs_prof')


    do n=1,size(obs_in)
       call nullify_obs(obs_out(n))
       Obs_out(n)%nvar = Obs_in(n)%nvar
       Obs_out(n)%project = Obs_in(n)%project
       Obs_out(n)%probe = Obs_in(n)%probe
       Obs_out(n)%ref_inst = Obs_in(n)%ref_inst
       Obs_out(n)%wod_cast_num = Obs_in(n)%wod_cast_num
       Obs_out(n)%fix_depth = Obs_in(n)%fix_depth
       Obs_out(n)%ocn_vehicle = Obs_in(n)%ocn_vehicle
       Obs_out(n)%database_id = Obs_in(n)%database_id
       Obs_out(n)%levels = Obs_in(n)%levels
       Obs_out(n)%profile_flag = Obs_in(n)%profile_flag
       Obs_out(n)%profile_flag_s = Obs_in(n)%profile_flag_s       
       Obs_out(n)%lon = Obs_in(n)%lon
       Obs_out(n)%lat = Obs_in(n)%lat
       Obs_out(n)%accepted = Obs_in(n)%accepted
       ALLOCATE(Obs_out(n)%depth(Obs_in(n)%levels))
       Obs_out(n)%depth(:) = Obs_in(n)%depth(:)
       ALLOCATE(Obs_out(n)%data_t(Obs_in(n)%levels))
       Obs_out(n)%data_t(:) = Obs_in(n)%data_t(:)
       ALLOCATE(Obs_out(n)%flag_t(Obs_in(n)%levels))
       Obs_out(n)%flag_t(:) = Obs_in(n)%flag_t(:)
       if (ASSOCIATED(Obs_in(n)%data_s)) then
           ALLOCATE(Obs_out(n)%data_s(Obs_in(n)%levels))
           Obs_out(n)%data_s(:) = Obs_in(n)%data_s(:)
           ALLOCATE(Obs_out(n)%flag_s(Obs_in(n)%levels))
           Obs_out(n)%flag_s(:) = Obs_in(n)%flag_s(:)          
       endif
       
       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
       ALLOCATE(Obs_out(n)%k_index(Obs_in(n)%levels))          
       Obs_out(n)%k_index = Obs_in(n)%k_index
       ALLOCATE(Obs_out(n)%ms_t(Obs_in(n)%levels))          
       Obs_out(n)%ms_t = Obs_in(n)%ms_t
       if (ASSOCIATED(Obs_in(n)%ms_s)) then       
           ALLOCATE(Obs_out(n)%ms_s(Obs_in(n)%levels))          
           Obs_out(n)%ms_s = Obs_in(n)%ms_s
       endif
       


       Obs_out(n)%tdiff = Obs_in(n)%tdiff
       Obs_out(n)%nbr_index = Obs_in(n)%nbr_index
       Obs_out(n)%nbr_dist  = Obs_in(n)%nbr_dist
       if (ASSOCIATED(Obs_in(n)%Model_grid)) &
            Obs_out(n)%Model_grid => Obs_in(n)%Model_Grid
    enddo
 
end subroutine copy_obs_prof

  subroutine copy_obs_sfc(Obs_in, Obs_out)
    type(ocean_surface_type), dimension(:), intent(in) :: Obs_in
    type(ocean_surface_type), dimension(:), intent(inout) :: Obs_out


    return

  end subroutine copy_obs_sfc

  subroutine adjust_obs_error_profile(Prof)

    use time_manager_mod, only : get_time
    
    type(ocean_profile_type), dimension(:), intent(inout) :: Prof
    integer :: secs, days, n, k, secs_w, days_w, m
    real :: tfac, Ims
    
    do n=1,size(Prof)
       call get_time(Prof(n)%tdiff,secs, days)
       m=Prof(n)%probe
       call get_time(time_window(m),secs_w,days_w)
       tfac = (days + secs/86400.) / days_w
       tfac = 1. - min(1.,tfac)
       if (tfac > 1.0 ) call mpp_error(FATAL)
       if (tfac < 0.0 ) call mpp_error(FATAL)
       do k=1,Prof(n)%levels
           Prof(n)%ms_t(k) = 1.0/ max(1.e-1,tfac) * Prof(n)%ms_t(k)
           if (ASSOCIATED(Prof(n)%data_s)) then
               Prof(n)%ms_s(k) = 1.0/ max(1.e-1,tfac) * Prof(n)%ms_s(k)
          endif
       end do
    end do

  end subroutine adjust_obs_error_profile

  subroutine adjust_obs_error_sfc(Diff)

    type(ocean_surface_type), intent(inout) :: Diff

    return
    
  end subroutine adjust_obs_error_sfc


  subroutine mult_obs_I_mse_profile(Obs)

    type(ocean_profile_type), dimension(:), intent(inout) :: Obs

    integer :: n,k
    real :: Ims
    
    do n=1,size(Obs)
       do k = 1, Obs(n)%levels
          Ims = 1/Obs(n)%ms_t(k)
          if (Obs(n)%flag_t(k) .eq. 0) Obs(n)%data_t(k) = Ims*Obs(n)%data_t(k)
          if (ASSOCIATED(Obs(n)%data_s)) then
              Ims = 1/Obs(n)%ms_s(k)
              if (Obs(n)%flag_s(k) .eq. 0) Obs(n)%data_s(k) = Ims*Obs(n)%data_s(k)
          endif
       end do
    end do
    
  end subroutine mult_obs_I_mse_profile

  subroutine mult_obs_I_mse_sfc(a, Obs)

    real, dimension(:), intent(in) :: a
    type(ocean_surface_type), intent(inout) :: Obs

  end subroutine mult_obs_I_mse_sfc


       

 ! </SUBROUTINE>
! <FUNCTION NAME="lowercase">
!   <DESCRIPTION>
! Turn a string from uppercase to lowercase, do nothing if the
! string is already in lowercase.
!   </DESCRIPTION>
 function lowercase (cs) 
 character(len=*), intent(in) :: cs
 character(len=len(cs))       :: lowercase 
 character :: ca(len(cs)) 

 integer, parameter :: co=iachar('a')-iachar('A') ! case offset
    
    ca = transfer(cs,"x",len(cs)) 
    where (ca >= "A" .and. ca <= "Z") ca = achar(iachar(ca)+co) 
    lowercase = transfer(ca,cs) 
    
  end function lowercase

  subroutine init_observations(localize)  

    use fms_mod, only : open_namelist_file,close_file,check_nml_error
    use mpp_io_mod, only : mpp_open, MPP_ASCII, MPP_RDONLY, MPP_MULTI, MPP_SINGLE
    use mpp_domains_mod, only : mpp_global_field
    
    logical, intent(in), optional :: localize

    integer :: data_window = 15 ! default data half-window is 15 days

    integer :: i,j
    
    
    type obs_entry_type
       character(len=128) :: filename
       character(len=16)  :: file_type
    end type obs_entry_type

    namelist /ocean_obs_nml/ data_window, max_prof_spacing, min_prof_depth
    
    character(len=128) :: input_files(max_files) = ''
    integer :: nfiles, filetype(max_files), ioun, io_status, ierr,&
                unit, nrecs, n
    character(len=256) :: record
    type(obs_entry_type) :: tbl_entry

    ioun = open_namelist_file()
    read(ioun,nml=ocean_obs_nml,iostat = io_status)
    ierr = check_nml_error(io_status,'ocean_obs_nml')
    call close_file(ioun)    

    time_window(:) = set_time(0,data_window)

    nfiles=0
    
    if (file_exist('ocean_obs_table') ) then
        call mpp_open(unit,'ocean_obs_table',action=MPP_RDONLY)
        nfiles = 0;nrecs=0
        do while (nfiles <= max_files)
           read(unit,'(a)',end=99,err=98) record
           nrecs=nrecs+1
           if (record(1:1) == '#') cycle
           read(record,*,err=98,end=98) tbl_entry
           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 ('idealized')
               filetype(nfiles)    = IDEALIZED_PROFILES
           case default
               call mpp_error(FATAL,'error in obs_table entry format')
           end select
98         continue
        enddo
        call mpp_error(FATAL,' number of obs files exceeds max_files parameter')
99      continue

    endif
    num_profiles=0
    num_sfc_obs=0
    
! get local indices for Model grid
! Since we currently only support regular grids, the
! input 2-d grid array is converted to 1-d
! halo points are added

!    xhalo=isc-isd
!    yhalo=jsc-jsd
!    if (xhalo.ne.ied-iec) call mpp_error(FATAL)
!    if (yhalo.ne.jed-jec) call mpp_error(FATAL)    

    allocate(x_grid(isg-1:ieg+1))
    allocate(y_grid(jsg-1:jeg+1))


    x_grid(isg:ieg) = Grd%x(isg:ieg,jsg)
    y_grid(jsg:jeg) = Grd%y(ieg/4,jsg:jeg)

    allocate(Profiles(max_profiles))

     if (Grd%cyclic) then
         x_grid(isg-1) = x_grid(ieg) - 360. ! assume grid is modulo 360 which is reasonable for data assimilation
         x_grid(ieg+1) = x_grid(isg) + 360.
     else
         x_grid(isg-1) = x_grid(isg) - 1.e-10
         x_grid(ieg+1) = x_grid(ieg) + 1.e-10
     endif

     y_grid(jsg-1) = y_grid(jsg) - 1.e-10
     y_grid(jeg+1) = y_grid(jeg) + 1.e-10
     
    
    do n=1, nfiles
       select case (filetype(n))
       case (PROFILE_FILE)
           call open_profile_dataset(trim(input_files(n)), localize)
       case (IDEALIZED_PROFILES)
           call create_ideal_profiles(localize)
       case default
          call mpp_error(FATAL,'filetype not currently supported')
       end select
    enddo


    return

  end subroutine init_observations
  
   subroutine add_tidal_error(Prof)
! NOT IMPLEMENTED YET !!!
     type(ocean_profile_type), intent(inout) :: Prof

     integer :: k
     real :: dtdz, err, a1, a2
    
     if (.not.ASSOCIATED(prof%ms_t)) then
         allocate(prof%ms_t(prof%levels))
         prof%ms_t(:) = min_obs_err_t
     endif
     
     do k=2,prof%levels - 1
        if (prof%flag_t(k-1) .eq. 0 .and. prof%flag_t(k+1) .eq. 0) then 
            dtdz = (prof%data_t(k+1)-prof%data_t(k-1))/(prof%depth(k+1)-prof%depth(k-1))
            a1 = abs(dtdz) * eta_tide_const
            err = max(a1,min_obs_err_t)
            prof%ms_t(k) = err*err
            if (ASSOCIATED(prof%data_s)) then
                dtdz = (prof%data_s(k+1)-prof%data_s(k-1))/(prof%depth(k+1)-prof%depth(k-1))
                a1 = abs(dtdz) * eta_tide_const
                err = max(a1,min_obs_err_s)
                prof%ms_s(k) = err*err
            endif
        endif
     enddo

   end subroutine add_tidal_error

  subroutine create_ideal_profiles(localize)
!
    use field_manager_mod, only: MODEL_OCEAN, parse, find_field_index, get_field_methods, method_type, get_field_info
    
    logical, intent(in), optional :: localize
    logical :: localize_data = .true.
    integer, parameter :: nlevels = 100 ! number of vertical levels for idealized profiles
    real, parameter :: width_trans = 250.0 ! with over which to transition from sfc value to bottom value
    real, parameter :: bot_depth = 2000.0 ! bottom depth for idealized profiles
    real, allocatable, dimension(:) :: lon,lat, sst, sss, bot_temp, bot_salt, depth
    real, allocatable, dimension(:,:) :: temp, salt, temp_error, salt_error
    integer, allocatable, dimension(:) :: yr, mon, day, hr, mm, ss
    integer :: nstation, unit, n, noobs, i, k
    real :: ri0,rj0,rk0, mid_depth, dtdf, temp_cent, depthC_I_trans, dsdf, salt_cent
    type(time_type) :: profile_time
    logical :: data_is_local
    integer :: model, parse_ok, cpe
    integer :: i0,j0,k0
    real :: dz, a, dx1, dx2, dy1, dy2
    
    real :: temp_missing=missing_value,salt_missing=missing_value,depth_missing=missing_value
    character(len=32) :: fld_type, fld_name
    type(method_type), allocatable, dimension(:) :: ocean_obs_methods
    
    if (PRESENT(localize)) localize_data = localize


    cpe = mpp_pe()

    dx1 = (x_grid(isc)-x_grid(isc-1))/2.0
    dx2 = (x_grid(iec+1)-x_grid(iec))/2.0
    dy1 = (y_grid(jsc)-y_grid(jsc-1))/2.0
    dy2 = (y_grid(jec+1)-y_grid(jec))/2.0
    
    model = model_ocean
    n = find_field_index(model,'ideal_profiles')
    call get_field_info(n,fld_type,fld_name,model,noobs)

    allocate(ocean_obs_methods(noobs))
    allocate(lon(noobs),lat(noobs), yr(noobs), mon(noobs), day(noobs), &
         hr(noobs), mm(noobs), ss(noobs), &
         sst(noobs), sss(noobs), bot_temp(noobs), bot_salt(noobs))
    allocate(temp(noobs,nlevels), salt(noobs,nlevels), temp_error(noobs,nlevels), salt_error(noobs,nlevels))
    allocate(depth(nlevels))
    
    call get_field_methods(n,ocean_obs_methods)
    do i=1,noobs
       parse_ok = parse(ocean_obs_methods(i)%method_control,'lon',lon(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')
       if (lon(i) .lt. x_grid(isg) ) lon(i) = lon(i) + 360.0
       if (lon(i) .gt. x_grid(ieg) ) lon(i) = lon(i) - 360.0       
       parse_ok = parse(ocean_obs_methods(i)%method_control,'lat',lat(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')
       parse_ok = parse(ocean_obs_methods(i)%method_control,'yr',yr(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')              
       parse_ok = parse(ocean_obs_methods(i)%method_control,'mon',mon(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')
       parse_ok = parse(ocean_obs_methods(i)%method_control,'day',day(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')
       parse_ok = parse(ocean_obs_methods(i)%method_control,'hr',hr(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')
       parse_ok = parse(ocean_obs_methods(i)%method_control,'mm',mm(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')
       parse_ok = parse(ocean_obs_methods(i)%method_control,'ss',ss(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')
       parse_ok = parse(ocean_obs_methods(i)%method_control,'sst',sst(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')
       parse_ok = parse(ocean_obs_methods(i)%method_control,'sss',sss(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')
       parse_ok = parse(ocean_obs_methods(i)%method_control,'bot_temp',bot_temp(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')
       parse_ok = parse(ocean_obs_methods(i)%method_control,'bot_salt',bot_salt(i))
       if (parse_ok == 0) call mpp_error(FATAL,'==>Error oda_core_mod: idealized_ocean_profiles table entry error')
    enddo

    if (noobs == 0 ) then
        call mpp_error(FATAL,'==> NOTE from oda_core_mod: no idealized profiles given in field table')
        return
    endif

    dz = bot_depth/(nlevels-1)
    
    do k=1,nlevels
       depth(k) = (k-1)*dz
    enddo
    
    mid_depth = bot_depth/2.0
    depthC_I_trans = mid_depth / width_trans
    do i=1,noobs
       dtdf = (bot_temp(i) - sst(i)) / (2.0*atan(1.0) + atan(depthC_I_trans))
       temp_cent = sst(i) + dtdf * atan(depthC_I_trans)
       temp(i,1) = sst(i)
       do k=2,nlevels-1
          temp(i,k) = temp_cent + dtdf * atan((depth(k)  - mid_depth)/width_trans)
       enddo
       temp(i,nlevels) = bot_temp(i)

       dsdf = (bot_salt(i) - sss(i)) / (2.0*atan(1.0) + atan(depthC_I_trans))
       salt_cent = sss(i) + dsdf * atan(depthC_I_trans)
       salt(i,1) = sss(i)
       do k=2,nlevels-1
          salt(i,k) = salt_cent + dsdf * atan((depth(k)  - mid_depth)/width_trans)
       enddo
       salt(i,nlevels) = bot_salt(i)
    enddo


    num_profiles=0
    do i=1,noobs

       data_is_local = .false.


! localized data is within region bounded by halo points
! (halo size = 1) adjacent to boundary points of computational domain
       
       if (lon(i) >= x_grid(isc-1) .and. &
            lon(i) <  x_grid(iec+1) .and. &
            lat(i) >= y_grid(jsc-1) .and. &
            lat(i) <  y_grid(jec+1)) data_is_local = .true.

      
       profile_time = set_date(yr(i),mon(i),day(i),hr(i),mm(i),ss(i))
       
       if ( data_is_local .OR. .NOT.localize_data) then
           if (lon(i) >= x_grid(isc)-dx1 .and. &
                lon(i) < x_grid(iec)+dx2 .and. &
                lat(i) >= y_grid(jsc)-dy1 .and. &
                lat(i) <  y_grid(jec)+dy2) then
!               nprof_in_comp_domain(cpe) = nprof_in_comp_domain(cpe)+1
           endif
           num_profiles=num_profiles+1
           if (num_profiles > max_profiles) then
               call mpp_error(FATAL,'maximum number of profiles exceeded.  Resize parameter max_profiles in ocean_obs_mod')
           endif


           profiles(num_profiles)%Model_Grid => Grd
           profiles(num_profiles)%nvar = 2
           profiles(num_profiles)%profile_flag = 0
           profiles(num_profiles)%profile_flag_s = 0
           profiles(num_profiles)%accepted = .true.
           allocate(profiles(num_profiles)%depth(nlevels))
           profiles(num_profiles)%depth=depth(1:nlevels)
           allocate(profiles(num_profiles)%data_t(nlevels))
           profiles(num_profiles)%data_t=temp(i,:)
           allocate(profiles(num_profiles)%flag_t(nlevels))
           profiles(num_profiles)%flag_t= 0               

           allocate(profiles(num_profiles)%data_s(nlevels))
           profiles(num_profiles)%data_s=salt(i,:)
           allocate(profiles(num_profiles)%flag_s(nlevels))
           profiles(num_profiles)%flag_s= 0

           profiles(num_profiles)%probe = 0
           profiles(num_profiles)%levels = nlevels
           profiles(num_profiles)%lat = lat(i)
           profiles(num_profiles)%lon = lon(i)
           allocate(profiles(num_profiles)%ms_t(nlevels))
           profiles(num_profiles)%ms_t(:) = min_obs_err_t ! default error variance for temperature
           allocate(profiles(num_profiles)%ms_s(nlevels))
           profiles(num_profiles)%ms_s(:) = min_obs_err_s  ! default error variance for salinity

           profiles(num_profiles)%time = profile_time
           
! calculate interpolation coefficients (make sure to account for grid offsets here!)
! note that this only works for lat/lon grids
           
           ri0 = frac_index(lon(i), x_grid(isg-1:ieg+1)) - 1 
           rj0 = frac_index(lat(i), y_grid(jsg-1:jeg+1)) - 1
           i0 = floor(ri0)
           j0 = floor(rj0)
           Profiles(num_profiles)%i_index = ri0
           Profiles(num_profiles)%j_index = rj0           
           Profiles(num_profiles)%accepted = .true.
           if (i0 < 0 .or. j0 < 0) then
               Profiles(num_profiles)%accepted = .false.
           endif
           if (i0 > ieg+1 .or. j0 > jeg+1) then
               call mpp_error(FATAL,'grid lat/lon index is out of bounds ')
           endif
           if (i0 < isc-1 .or. i0 > iec) then
               call mpp_error(FATAL,'grid lat/lon index is out of bounds ')
           endif
           if (j0 < jsc-1 .or. j0 > jec) then
               call mpp_error(FATAL,'grid lat/lon index is out of bounds ')
           endif
           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.
               endif
           endif


           
           if (Profiles(num_profiles)%accepted) then
               allocate(Profiles(num_profiles)%k_index(Profiles(num_profiles)%levels))
               do k=1, Profiles(num_profiles)%levels
                  rk0 = frac_index(Profiles(num_profiles)%depth(k), Grd%z(:))
                  k0 = floor(rk0)
                  if ( k0 == -1) then
                      if (Profiles(num_profiles)%depth(k) .ne. missing_value .and. &
                           Profiles(num_profiles)%depth(k) .lt. Grd%z(1)) then
                           k0 = 1
                           rk0 = 1.0
                       endif
                   endif


                  if (k0 .gt. size(Grd%z)-1 ) then
                      write(*,*) 'k0 out of bounds, rk0,k0= ',rk0,k0
                     write(*,*) 'Z_bound= ',Grd%z_bound
                      write(*,*) 'Profile%depth= ',Profiles(num_profiles)%depth
                      call mpp_error(FATAL)
                  endif
                  
                  Profiles(num_profiles)%k_index(k) = rk0
                  
                  if (Profiles(num_profiles)%flag_t(k) .eq. 0) 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_t(k) = 1
                      endif
                      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_t(k) = 1
                      endif
                      if (Profiles(num_profiles)%data_t(k) == missing_value &
                         .or. Profiles(num_profiles)%depth(k) == missing_value) then
                          Profiles(num_profiles)%flag_t(k) = 1
                      endif
                  endif
                  
               enddo
           endif          
       endif

    enddo

!    a = nprof_in_comp_domain(cpe)

!    call mpp_broadcast(nprof_in_comp_domain(cpe),cpe)
    
!    call mpp_sum(a)

!    write(stdout(),*) 'A grand total of ',a,' profiles satisify acceptance criteria'

!    do i=0,mpp_npes()-1
!       write(stdout(),*) 'pe=',i,'profile count=',nprof_in_comp_domain(i)
!    enddo
    
  end subroutine create_ideal_profiles


  subroutine nullify_obs_prof(profile)

    type(ocean_profile_type), intent(inout) :: profile


    profile%nvar = 0
    profile%project=0
    profile%probe=0
    profile%ref_inst=0
    profile%wod_cast_num=0
    profile%fix_depth=0
    profile%ocn_vehicle=0
    profile%database_id=0
    profile%levels=0
    profile%profile_flag=-1
    profile%profile_flag_s=-1
    profile%lon=-1.0e10
    profile%lat=-1.0e10
    profile%accepted=.false.
    profile%nlinks=0
    if (ASSOCIATED(profile%next)) profile%next=>NULL()
    if (ASSOCIATED(profile%depth)) profile%depth=>NULL()
    if (ASSOCIATED(profile%data_t)) profile%data_t=>NULL()
    if (ASSOCIATED(profile%data_s)) profile%data_s=>NULL()
    if (ASSOCIATED(profile%flag_t)) profile%flag_t=>NULL()
    if (ASSOCIATED(profile%flag_s)) profile%flag_s=>NULL()
    profile%temp_err=0.0
    profile%salt_err=0.0
    if (ASSOCIATED(profile%ms_t)) profile%ms_t=>NULL()
    if (ASSOCIATED(profile%ms_s)) profile%ms_s=>NULL()    
    profile%time = set_time(0,0)
    profile%yyyy = 0
    profile%mmdd = 0
    if (ASSOCIATED(profile%model_time)) profile%model_time=>NULL()
    if (ASSOCIATED(profile%model_grid)) profile%model_grid=>NULL()
    if (ASSOCIATED(profile%k_index)) profile%k_index=>NULL()
    profile%i_index=-1.0
    profile%j_index=-1.0
    profile%tdiff = set_time(0,0)
    
  end subroutine nullify_obs_prof
  
end module oda_core_mod