function status=diff_coef(STDname,Kh,Kv);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Hernan G. Arango %%%%
% Copyright (c) 2006 Rutgers University.                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                      %
% function status=diff_coef(STDname,Kh,Kv);                            %
%                                                                      %
% This function writes horizontal and vertical diffusion coefficients, %
% used in data assimilation spatial convolution, into background/model %
% error covariance standard deviation NetCDF file. If need, it defines %
% such NetCDF variables.                                               %
%                                                                      %
% On Input:                                                            %
%                                                                      %
%    STDname     Background/model error covariance standard deviation  %
%                  file name (character string).                       %
%    Kh          Convolution horizontal diffusion coefficients (m/s)   %
%                  at RHO-points. This is a 2D array, 0 < Kh <= 1.     %
%                  If variable not provided, a uniform unity field     %
%                  is written.                                         % 
%    Kv          Convolution vertical diffusion coefficients (m/s)     %
%                  at W-points. This is a 3D array, 0 < Kv <= 1.       %
%                  If variable not provided, a uniform unity field     %
%                  is written.                                         % 
%                                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%-----------------------------------------------------------------------
% Inquire NetCDF file about convolution diffusion coefficients.
%-----------------------------------------------------------------------

got.Kh=0;  define.Kh=1;  Vname.Kh='Kh';
got.Kv=0;  define.Kv=1;  Vname.Kv='Kv';
is3d=0;

[varnam,nvars]=nc_vname(STDname);
for n=1:nvars,
  name=deblank(varnam(n,:));
  switch name
    case {Vname.Kh}
      got.Kh=1;
      define.Kh=0;
    case {Vname.Kv}
      got.Kv=1;
      define.Kv=0;
    case 'temp'
      is3d=1;
  end,
end,

%  Inquire about dimensions.

  Dname.xr='xi_rho';    Dname.yr='eta_rho';
  Dname.sw='s_w';

  [Dnames,Dsizes]=nc_dim(STDname);
  ndims=length(Dsizes);
  for n=1:ndims,
    dimid=n-1;
    name=deblank(Dnames(n,:));
    switch name
      case {Dname.xr}
        Lp=Dsizes(n);
        Dsize.xr=Dsizes(n);
        did.xr=dimid;
      case {Dname.yr}
        Mp=Dsizes(n);
        Dsize.yr=Dsizes(n);
        did.yr=dimid;
      case {Dname.sw}
        Nw=Dsizes(n);  
        Dsize.sw=Dsizes(n);
        did.sw=dimid;
    end,
  end,

%-----------------------------------------------------------------------
%  If appropriate, define convolution diffusion variables.
%-----------------------------------------------------------------------

if (define.Kh | define.Kv),

%  Define NetCDF parameters.

  [ncglobal]=mexcdf('parameter','nc_global');
  [ncdouble]=mexcdf('parameter','nc_double');
  [ncfloat]=mexcdf('parameter','nc_float');
  [ncchar]=mexcdf('parameter','nc_char');

%  Open STD NetCDF file.

  [ncid]=mexcdf('ncopen',STDname,'nc_write');
  if (ncid == -1),
    error(['DIFF_COEF: ncrefdef - unable to open file: ', STDname]);
    return
  end,

%  Put STD NetCDF file into define mode.

  [status]=mexcdf('ncredef',ncid);
  if (status == -1),
    error(['DIFF_COEF: ncrefdef - unable to put into define mode.']);
    return
  end,

  spherical=nc_read(STDname,'spherical');

%  Define convolution horizontal diffusion coefficients.

  if (define.Kh),
    Var.name =Vname.Kh;
    Var.type =ncdouble;
    Var.dimid=[did.yr did.xr];
    Var.long ='convolution horizontal diffusion coefficients';
    Var.units='meter2 second-1';
    Var.time ='ocean_time';
    if (spherical == 'T' | spherical == 't'),
      Var.coord='lon_rho lat_rho';
    else,
      Var.coord='x_rho y_rho';
    end,
    Var.field='Kh, scalar, series';

    [varid,status]=nc_vdef(ncid,Var);
    clear Var
  end,

%  Define convolution vertical diffusion coefficients.

  if (define.Kv),
    Var.name =Vname.Kv;
    Var.type =ncdouble;
    Var.dimid=[did.sw did.yr did.xr];
    Var.long ='convolution vertical diffusion coefficients';
    Var.units='meter2 second-1';
    Var.time ='ocean_time';
    if (spherical == 'T' | spherical == 't'),
      Var.coord='lon_rho lat_rho';
    else,
      Var.coord='x_rho y_rho';
    end,
    Var.field='Kv, scalar, series';

    [varid,status]=nc_vdef(ncid,Var);
    clear Var
  end,

%  Leave definition mode.

  [status]=mexcdf('ncendef',ncid);
  if (status == -1),
    error(['DIFF_COEF: ncendef - unable to leave definition mode.']);
  end,

%  Close STD NetCDF file.

  [status]=mexcdf('ncclose',ncid);
  if (status == -1),
    error(['DIFF_COEF: ncclose - unable to close NetCDF file: ', STDname]);
  end,

end,

%-----------------------------------------------------------------------
%  Write out convolution diffusion coefficients into STD NetCDF file.
%-----------------------------------------------------------------------

if (nargin < 3),
  Kv(1:Lp,1:Mp,1:Nw)=1.0;
end,
if (nargin < 2),
  Kh(1:Lp,1:Mp)=1.0;
end,

[status]=nc_write(STDname,Vname.Kh,Kh);
if (is3d),
  [status]=nc_write(STDname,Vname.Kv,Kv);
end,

return
