% --- after running this script with Matlab use cdo as below
% --- cdo replace NH3_geos.4x5.nc NH3_geos.4x5.tmp_seabirds.nc NH3_geos.4x5.with_seabirds.nc
% ======= final step: edit HEMCO_Config.rc to replace NH3_geos.4x5.nc with NH3_geos.4x5.with_seabirds.nc
% If using this data please cite
%   https://data.gov.uk/dataset/global-ammonia-emissions-from-seabirds
%   and
%   If you reuse this data, you should cite: Riddick, S.N., Dragosits, U., Blackall, T.D., Daunt, F., Wanless, S., Sutton, M.A. (2012). Global ammonia emissions from seabirds.

% --- read the existing GEOS-Chem input file (natural NH3 sources)
NH3=double(ncread('/data10/ctm/HEMCO/NH3/v2014-07/NH3_geos.4x5.nc','NH3_NATU'));
lat=double(ncread('/data10/ctm/HEMCO/NH3/v2014-07/NH3_geos.4x5.nc','lat'));
lon=double(ncread('/data10/ctm/HEMCO/NH3/v2014-07/NH3_geos.4x5.nc','lon'));

% --- create a new netCDF file
nccreate('NH3_geos.4x5.tmp_seabirds.nc','lon',...
	 'Dimensions',{'lon',72},...
         'Format','classic') 
nccreate('NH3_geos.4x5.tmp_seabirds.nc','lat',...
	 'Dimensions',{'lat',46},...
         'Format','classic') 
nccreate('NH3_geos.4x5.tmp_seabirds.nc','time',...
	 'Dimensions',{'time',inf},...
        'Format','classic') 
nccreate('NH3_geos.4x5.tmp_seabirds.nc','NH3_NATU',...
	 'Dimensions',{'lon',72,'lat',46,'time',12},...
         'Format','classic') 

% -------------------------------------------------------------------
% --- find edges of the GEOS-Chem grid boxes
  lat_edges(1)=-90;
lat_edges(length(lat)+1)=90;
for i=1:length(lat)-1
  lat_edges(i+1)=(lat(i)+lat(i+1))/2;
end

for i=1:length(lon)-1
  lon_edges(i)=(lon(i+1)+lon(i))/2;
end
lon_edges(length(lon))=177.5;

% -------------------------------------------------------------------
% --- read the text file with seabird emissions
load /data1/carlin/seabirds/final_all_bird_emissions.dat

% --- put the seabird emissions on a 4x5 degree grid
% --- use column 5 - the moderate scenario
% --- lon greater than 177.5 or less than -177.5 will go into bin 1

seabird_emiss=double(zeros(72,46));
final_all_bird_emissions=double(final_all_bird_emissions);
count=[];
for i=1:length(final_all_bird_emissions(:,1))
 if final_all_bird_emissions(i,1) >= 50  % --------------- here we select only emissions north of 50N
  for j=1:length(lat_edges)-1
% find the lat index
   if final_all_bird_emissions(i,1) <= lat_edges(j+1) && final_all_bird_emissions(i,1) > lat_edges(j)
       flag=0;
    for k=1:length(lon_edges)-1
% find the lon index
      if final_all_bird_emissions(i,2) <= lon_edges(k+1) && final_all_bird_emissions(i,2) > lon_edges(k)
       seabird_emiss(k,j)=seabird_emiss(k,j)+final_all_bird_emissions(i,5);
       flag=1;
       count=[count i];
      end
    end % end lon
    if flag==0 % a lon box was not found
     if final_all_bird_emissions(i,2) > 177.5 || final_all_bird_emissions(i,2) < -177.5
      seabird_emiss(1,j)=seabird_emiss(1,j)+final_all_bird_emissions(i,5);
      count=[count i];
     end
    end
   end
  end % end lat
 end
end

% --- get a 4x5 GEOS-Chem area file
    GCFile = sprintf('/scratch-space1/carlin/GC_output/ms_std/bpch2hdf_output.20110601.000000.hdf')
    SPEC   = sprintf('DXYP::DXYP')
    AREA=hdfread(GCFile,SPEC );


% --- distribute the emissions over May-September i.e. 5 months  - although May 15 to Sept 15 may be better
% --- divide by area of box and number of seconds in 5 months (may-sep) ---> units change from  Mg NH3/box/yr to kg/m2/s
for i=1:46
  for j=1:72
  seabird_emiss(j,i)=1e3*seabird_emiss(j,i)/(AREA(i,j)*153*24*3600);
end
end

% --- put into a matrix that is the same size as NH3_NATU
  seabird_emiss_monthly=double(zeros(72,46,12));
  for i=5:9
    seabird_emiss_monthly(:,:,i)=seabird_emiss(:,:);
  end
  seabird_emiss_monthly=double(seabird_emiss_monthly);

% ---- add seabird emissions to the existing natural NH3 input data
NH3_NEW=NH3+seabird_emiss_monthly;

time = [0:1:11]
  ncwrite('NH3_geos.4x5.tmp_seabirds.nc','time',time)
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','time','standard_name','time')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','time','long_name','time')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','time','units','months since 1990-01-01 00:00:00')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','time','calendar','standard')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','time','axis','T')

  ncwrite('NH3_geos.4x5.tmp_seabirds.nc','lat',lat)
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','lat','standard_name','Latitude')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','lat','long_name','Latitude')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','lat','units','degrees_north')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','lat','axis','Y')

  ncwrite('NH3_geos.4x5.tmp_seabirds.nc','lon',lon)
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','lon','standard_name','Longitude')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','lon','long_name','Longitude')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','lon','units','degrees_east')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','lon','axis','X')

  ncwrite('NH3_geos.4x5.tmp_seabirds.nc','NH3_NATU',NH3_NEW)
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','NH3_NATU','long_name','natural + seabird NH3 ')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','NH3_NATU','units','kg/m2/s')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','NH3_NATU','grid_type','gaussian')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','NH3_NATU','description','natural + seabird NH3 emissions in kg NH3/m2/s')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','/','Conventions','COARDS')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','/','Title','COARDS/netCDF file')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','/','Delta_Lon','5.f')
  ncwriteatt('NH3_geos.4x5.tmp_seabirds.nc','/','Delta_Lat','4.f')


