function [NLobsCost,TLobsCost]=obs_cost(ObsFile,ModFile);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2008 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% function [iter,Jo,Jb,GradJ]=cost_function(fname,normalize)                %
%                                                                           %
% This function computes nonlinear and tangent linear models cost           %
% functions.                                                                %
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    ObsFile     Observations NetCDF file name (character string).          %
%    ModFile     ROMS 4DVAR NetCDF file name (character string).            %
%                                                                           %
% On Output:                                                                %
%                                                                           %
%    NLobsCost   Nonlinear model cost function.                             %
%    TLobsCost   Tangent linear model cost function.                        %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Read in observation data.

ObsVal=nc_read(ObsFile,'obs_value');
ObsErr=nc_read(ObsFile,'obs_error');
ObsType=nc_read(ObsFile,'obs_type');

%  Read in 4DVAR data.

ObsScale=nc_read(ModFile,'obs_scale');
NLmodVal=nc_read(ModFile,'NLmodel_value');
TLmodVal=nc_read(ModFile,'TLmodel_value');

%  Compute cost functions.

Nobs=length(ObsVal);
NstateVar=length(nc_read(ModFile,'Nobs'));

for ivar=1:NstateVar,
  NLobsCost(ivar)=0;
  TLobsCost(ivar)=0;
end,

for iobs=1:Nobs,

  ivar=ObsType(iobs)+1;

  cff=0.5*ObsScale(iobs)/ObsErr(iobs);
  nlval=cff*(NLmodVal(iobs)-ObsVal(iobs))^2;
  tlval=cff*(NLmodVal(iobs)+TLmodVal(iobs)-ObsVal(iobs))^2;

  NLobsCost(1)=NLobsCost(1)+nlval;
  NLobsCost(ivar)=NLobsCost(ivar)+nlval;

  TLobsCost(1)=TLobsCost(1)+nlval;
  TLobsCost(ivar)=TLobsCost(ivar)+nlval;
end,

%  Report cost function values.

variable={'Total' 'zeta' 'ubar' 'vbar' 'u' 'v' 'temp' 'salt'};

disp('  ');
for ivar=1:NstateVar,
  if (NLobsCost(ivar) > 0 & TLobsCost(ivar) > 0),
    disp(['   NLM Cost = ', sprintf('%15.4f', NLobsCost(ivar)), ...
	  '   TLM Cost = ', sprintf('%15.4f', TLobsCost(ivar)), ...
	  '   ', char(variable(ivar))]);
  end,
end,
disp('  ');

return

