function [fldpseudo]=mon2pseudo(fld,startmonth,cyclic)  

if not(exist('cyclic'))
  cyclic=false;
end 

nrec=size(fld,3)
monlen=[31 28 31 30 31 30 31 31 30 31 30 31]; 
for m=1:nrec 
  dm(m)=monlen(mod(m-2+startmonth,12)+1);
end	

% compute monthly averages based on initial field
fldave=nan(size(fld));
for m=1:nrec
  % set month indices and compute weight-coefficients 
  mm1=mod(m-2,12)+1;
  mp1=mod(m,12)+1;
  Cmm1=1/4*dm(m)/(dm(m)+dm(mm1));
  Cmp1=1/4*dm(m)/(dm(m)+dm(mp1));
  Cm=1/4*(dm(mm1)/(dm(m)+dm(mm1))+dm(mp1)/(dm(m)+dm(mp1)))+1/2;
  % compute monthly means assuming input fields are mid-month 
  fldave(:,:,m)=Cmm1*fld(:,:,mm1)+Cm*fld(:,:,m)+Cmp1*fld(:,:,mp1);
end
% replace first and last record with true average 
if not(cyclic)
  fldave(:,:,1)=fld(:,:,1);
  fldave(:,:,end)=fld(:,:,end);
end
% compute global RMSE for initial field
N=sum(~isnan(fldave(:)));
rms_before=sqrt(nansum((fld(:)-fldave(:)).^2)/N);

% compute pseudo, mid-month field 
niter=20; % number of iterations 
fldpseudo=fld;
fldmin=nanmin(fld(:));
fldmax=nanmax(fld(:));
for n=1:niter
  for m=2:nrec-1
    % set month indices and compute weight-coefficients 
    mm1=mod(m-2,12)+1;  
    mp1=mod(m,12)+1;
    Cmm1=1/4*dm(m)/(dm(m)+dm(mm1));
    Cmp1=1/4*dm(m)/(dm(m)+dm(mp1));
    Cm=1/4*(dm(mm1)/(dm(m)+dm(mm1))+dm(mp1)/(dm(m)+dm(mp1)))+1/2;
    % update guess for pseudo field 
    guess=1/Cm*fld(:,:,m)-Cmm1/Cm*fldpseudo(:,:,mm1)-Cmp1/Cm*fldpseudo(:,:,mp1); 
    guess(find(guess<fldmin))=fldmin;
    guess(find(guess>fldmax))=fldmax;
    fldpseudo(:,:,m)=guess;
  end

  % compute monthly averages based on initial and pseudo fields
  fldavepseudo=nan(size(fld));
  for m=2:nrec-1
    % set month indices and compute weight-coefficients 
    mm1=mod(m-2,12)+1;
    mp1=mod(m,12)+1;
    Cmm1=1/4*dm(m)/(dm(m)+dm(mm1));
    Cmp1=1/4*dm(m)/(dm(m)+dm(mp1));
    Cm=1/4*(dm(mm1)/(dm(m)+dm(mm1))+dm(mp1)/(dm(m)+dm(mp1)))+1/2;
    % compute monthly means using pseudo field
    fldavepseudo(:,:,m)=Cmm1*fldpseudo(:,:,mm1)+Cm*fldpseudo(:,:,m)+Cmp1*fldpseudo(:,:,mp1); 
  end
  % replace first and last record with true average 
  if not(cyclic)
    fldavepseudo(:,:,1)=fld(:,:,1);
    fldavepseudo(:,:,end)=fld(:,:,end);
  end 

  % compute global RMSE for initial field
  rms_after=sqrt(nansum((fld(:)-fldavepseudo(:)).^2)/N);
  disp(['Global RMSE computed from pseudo field=' sprintf('%20.15f',rms_after)]) 
  disp(['RMSE fraction=' sprintf('%20.15f',rms_after/rms_before*100) '%']) 

end
