% --- before running this script - use the following command to produce NH3_geos.4x5.with_seabirds.nc
%  cdo replace NH3_geos.4x5.nc NH3_geos.4x5.tmp_seabirds.nc NH3_geos.4x5.with_seabirds.nc
% --- after we make NH3_geos.4x5.tmp_seabirds.nc - then check what was added is same as in the original file
                                                                                                                                                                      
NH3new=double(ncread('NH3_geos.4x5.with_seabirds.nc','NH3_NATU'));                                                                                                                       
NH3old=double(ncread('NH3_geos.4x5.nc','NH3_NATU'));                                                                                                                                     
diffNH3=NH3new-NH3old;                                                                                                                                                                   

% find Mg/m2/month
may_seabirdNH3=squeeze(1e-3*diffNH3(:,:,5)*31*24*3600);                                                                                                                   
june_seabirdNH3=squeeze(1e-3*diffNH3(:,:,6)*30*24*3600);                                                                                                                  
july_seabirdNH3=squeeze(1e-3*diffNH3(:,:,7)*31*24*3600);                                                                                                                  
aug_seabirdNH3=squeeze(1e-3*diffNH3(:,:,8)*31*24*3600);                                                                                                                  
sep_seabirdNH3=squeeze(1e-3*diffNH3(:,:,9)*30*24*3600);                                                                                                                   
sum_seabirds=may_seabirdNH3+june_seabirdNH3+july_seabirdNH3+aug_seabirdNH3+sep_seabirdNH3;                                                                                                

AREA2=AREA';                                                                                                                                                                             
sum_seabirds=sum_seabirds.*AREA2;
sumN50_GC=sum(sum(sum_seabirds))

% --- compare with the total emissions of the original file                                                                                                                              
% --- sum_seabirds and sumN50 should be equal if emissions added correctly  (about 35.6 Gg NH3 is added each year north of 50N)                                                                                                                        
sumN50=0;                                                                                                                                                                                
for i=1:length(final_all_bird_emissions(:,1))                                                                                                                                            
if final_all_bird_emissions(i,1) >=50                                                                                                                                                    
sumN50=sumN50+final_all_bird_emissions(i,5);                                                                                                                                             
end                                                                                                                                                                                      
end
sumN50

