diff --git a/tools/rivers/bgc/ARC/Data/ArcticGro/ArcticGRO_Water_Quality_Data.xlsx b/tools/rivers/bgc/ARC/Data/ArcticGro/ArcticGRO_Water_Quality_Data.xlsx new file mode 120000 index 000000000..53601b94a --- /dev/null +++ b/tools/rivers/bgc/ARC/Data/ArcticGro/ArcticGRO_Water_Quality_Data.xlsx @@ -0,0 +1 @@ +../../../NEP/Data/ArcticGro/ArcticGRO_Water_Quality_Data.xlsx \ No newline at end of file diff --git a/tools/rivers/bgc/ARC/Data/ArcticGro/ArcticGro_Process.m b/tools/rivers/bgc/ARC/Data/ArcticGro/ArcticGro_Process.m new file mode 100644 index 000000000..fc9f33be8 --- /dev/null +++ b/tools/rivers/bgc/ARC/Data/ArcticGro/ArcticGro_Process.m @@ -0,0 +1,129 @@ +% Program to process the ArcticGro data + +clear all + +filename = 'ArcticGRO_Water_Quality_Data.xlsx' +% co2 system solver to derive DIC from alkalinity, pH and temperature +addpath /home/cas/matlab/co2sys + +% Longitude and Latitude from GlobalNEWS +river_name{1} = 'Ob'; +lon(1) = 68.5; lat(1) = 68.75; +river_name{2} = 'Yenisey'; +lon(2) = 82.25; lat(2) = 71.25; +river_name{3} = 'Lena'; +lon(3) = 128; lat(3) = 73; +river_name{4} = 'Kolyma'; +lon(4) = 161.25; lat(4) = 69.25; +river_name{5} = 'Yukon'; +lon(5) = -164.75; lat(5) = 62.75; +river_name{6} = 'Mackenzie'; +lon(6) = -134.75; lat(6) = 69.25; + + +for n = 1:6 + Data = readtable(filename,'Sheet',n); + discharge(n) = nanmean(Data{:,5}); + temp(n) = nanmean(Data{:,6}); + % alk in mg CaCO3/L ~ g CaCO3/m3 + alk(n) = nanmean(Data{:,9}); + alk(n) = alk(n)/(40+12+16*3)*2*1e3; % to milliequivalents per m-3 + % pH + pH(n) = nanmean(Data{:,7}); + % tdn in mg N L-1 ~ g N/m3 + tdn(n) = nanmean(Data{:,21}); + tdn(n) = tdn(n)*1e3/14; % mmoles m-3 + % no3 in micrograms N per L ~ mg N m-3 + no3(n) = nanmean(Data{:,22}); + no3(n) = no3(n)/14; % mmoles m-3 + % nh4 in micrograms N per L ~ mg N m-3 + nh4(n) = nanmean(Data{:,23}); + nh4(n) = nh4(n)/14; % mmoles m-3 + % tdp in micrograms P per L ~ mg P m-3 + tdp(n) = nanmean(Data{:,24}); + tdp(n) = tdp(n)/31; % mmoles m-3 + % po4 in micrograms P per L ~ mg P m-3 + po4(n) = nanmean(Data{:,25}); + po4(n) = po4(n)/31; % mmoles m-3 + % sio2 in mg SiO2 per L ~ g P m-3 + sio2(n) = nanmean(Data{:,26}); + sio2(n) = sio2(n)*1e3/(28.06+16*2); + % pon in micrograms N per L ~ mg N m-3 + pon(n) = nanmean(Data{:,47}); + pon(n) = pon(n)/14; %mmoles m-3 + + % calculate DIC from alk, pH and DIC + out = CO2SYS(alk(n),pH(n),1,3,0,temp(n),temp(n), ... + 100,100,0,0,0,0,4,15,1,2,2); + dic(n) = mean(out(:,2)); + + % calculate don from tdn, no3 and nh4 + aa = find(isfinite(Data{:,21}) & isfinite(Data{:,22}) & isfinite(Data{:,23})); + don(n) = nanmean(Data{aa,21}*1e3/14 - Data{aa,22}/14 - Data{aa,23}/14); + + % calculate dop from tdp and po4 + bb = find(isfinite(Data{:,24}/31) & isfinite(Data{:,25}/31)); + dop(n) = nanmean(Data{aa,24}/31 - Data{aa,25}/31); +end + +% Fill in particulate phosphorus from GLOBAL NEWS +pp(1) = 1.29; % Ob +pp(2) = 0.82; % Yenisey +pp(3) = 1.48; % Lena +pp(4) = 1.21; % Kolyma +pp(5) = 1.94; % Yukon +pp(6) = 1.44; % Mackenzie + +% For reference (DIN, DON, PN, DIP, DOP, PP) +% Yukon: 7.30 17.25 29.10 0.07 0.42 1.94 +% Mackenzie: 7.42 19.86 24.34 0.14 0.48 1.44 +% St. Lawrence: 52.80 24.86 3.81 0.98 0.55 0.21 +% Ob: 21.80 24.16 22.35 1.16 0.56 1.29 +% Lena: 7.74 21.34 25.48 0.24 0.52 1.48 +% Yenisey: 8.54 21.65 14.21 0.088 0.52 0.82 +% Kolyma: 10.17 21.92 21.27 0.18 0.53 1.21 + +lon_stations_arcticgro = lon; +lat_stations_arcticgro = lat; +station_names_arcticgro = river_name; + +Q_ann_arcticgro = discharge; +dic_ann_arcticgro = dic; +alk_ann_arcticgro = alk; +no3_ann_arcticgro = no3; +nh4_ann_arcticgro = nh4; +din_ann_arcticgro = no3_ann_arcticgro + nh4_ann_arcticgro; +pn_ann_arcticgro = pon; +don_ann_arcticgro = don; +dip_ann_arcticgro = po4; +dop_ann_arcticgro = dop; +pp_ann_arcticgro = pp; +si_ann_arcticgro = sio2; +o2_ann_arcticgro = ones(size(Q_ann_arcticgro))*NaN; +dfe_ann_arcticgro = ones(size(Q_ann_arcticgro))*NaN; +pfe_ann_arcticgro = ones(size(Q_ann_arcticgro))*NaN; + +% did not try and extract monthly data from arcticgro, so leave as NaNs +dic_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +alk_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +no3_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +nh4_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +din_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +don_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +pn_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +dip_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +dop_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +pp_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +dfe_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +pfe_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +si_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +o2_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; + +save arcticgro_data lon_stations_arcticgro lat_stations_arcticgro station_names_arcticgro ... + Q_ann_arcticgro dic_ann_arcticgro alk_ann_arcticgro no3_ann_arcticgro nh4_ann_arcticgro ... + din_ann_arcticgro pn_ann_arcticgro don_ann_arcticgro dip_ann_arcticgro dop_ann_arcticgro ... + pp_ann_arcticgro si_ann_arcticgro o2_ann_arcticgro dfe_ann_arcticgro pfe_ann_arcticgro ... + dic_monthly_arcticgro alk_monthly_arcticgro no3_monthly_arcticgro nh4_monthly_arcticgro ... + din_monthly_arcticgro don_monthly_arcticgro pn_monthly_arcticgro dip_monthly_arcticgro ... + dop_monthly_arcticgro pp_monthly_arcticgro dfe_monthly_arcticgro pfe_monthly_arcticgro ... + si_monthly_arcticgro o2_monthly_arcticgro diff --git a/tools/rivers/bgc/ARC/Data/GLORICH/Arctic_GLORICH_Process.m b/tools/rivers/bgc/ARC/Data/GLORICH/Arctic_GLORICH_Process.m new file mode 100644 index 000000000..5c50c3019 --- /dev/null +++ b/tools/rivers/bgc/ARC/Data/GLORICH/Arctic_GLORICH_Process.m @@ -0,0 +1,287 @@ +clear all; + +% load in matlab file with GLORICH data +% created with:T = readtable('hydrochemistry.csv','NumHeaderLines',1); +%load GLORICH_Data.mat; +T = readtable('hydrochemistry.csv','NumHeaderLines',1); +stations = T{:,1}; + +%"110009","Skeena River at Usk","BC08EF0001","Canada","BC",54.63,-128.41,"NA1983" +%"110004","Stikine River above Choquette River","BC08CF0002","Canada","BC",56.82,-131.76,"NA1983" +%"102954","COPPER R AT MILLION DOLLAR BRIDGE NR CORDOVA AK","15214000","USA","AK",60.67,-144.74,"NA1983" +%"102983","SUSITNA R AT SUSITNA STATION AK","15294350","USA","AK",61.54,-150.51,"NA1983" +%"102985","NUSHAGAK R AT EKWOK AK","15302500","USA","AK",59.34,-157.47,"NA1983" +%"102986","KUSKOKWIM R AT CROOKED CREEK AK","15304000","USA","AK",61.87,-158.10,"NA1983" +%"102988","YUKON R AT PILOT STATION AK","15565447","USA","AK",61.93,-162.88,"NA1983" +% this one had by far the most readings of the 3 Churchill stations +%"111762","CHURCHILL RIVER ABOVE UPPER MUSKRAT FALLS","NF03OE0001","Canada","New Foundland",53.24,-60.78,"NA1983" +%"101569","ST. LAWRENCE R AT CORNWALL ONT NR MASSENA NY","4264331","USA","NY",45.00,-74.79,"NA1983" + + +% Define rivers. I took the lat/lon from the mouths defined in Global NEWS +% Where stream gages are upstream, this makes the assumption that those +% properties are reasonably indicative of conditions at the river mouth. +num_rivers = 8; +river_name{1} = 'Skeena'; station_num(1) = 110009; +lat(1) = 54.25; lon(1) = 229.75; +river_name{2} = 'Stikine'; station_num(2) = 110004; +lat(2) = 56.75; lon(2) = 227.75; +river_name{3} = 'Copper'; station_num(3) = 102954; +lat(3) = 60.25; lon(3) = 215.25; +river_name{4} = 'Susitna'; station_num(4) = 102983; +lat(4) = 61.75; lon(4) = 209.75; +river_name{5} = 'Nushagak'; station_num(5) = 102985; +lat(5) = 58.75; lon(5) = 201.75; +river_name{6} = 'Kuskokwim'; station_num(6) = 102986; +lat(6) = 60.25; lon(6) = 197.75; +river_name{7} = 'Churchill'; station_num(7) = 111762; +lat(7) = 54.0; lon(7) = -59; +%river_name{8} = 'St. Lawrence'; station_num(8) = 120300; +river_name{8} = 'St. Lawrence'; station_num(8) = 101569; +%river_name{8} = 'St. Lawrence'; station_num(8) = 111573; +lat(8) = 60.25; lon(8) = 197.75; + +% co2 system solver to derive DIC from alkalinity, pH and temperature +addpath /home/cas/matlab/co2sys + +for n = 1:num_rivers + aa = find(stations == station_num(n)); + num_stations(n) = size(aa,1); + discharge(n) = nanmean(T{aa,6}); + alk_vec = T{aa,20}; + alk(n) = nanmean(alk_vec); + % nitrogen species + tn_vec = T{aa,62}; num_tn(n) = size(find(isfinite(tn_vec) == 1),1); tn(n) = nanmean(tn_vec); + dn_vec = T{aa,64}; num_dn(n) = size(find(isfinite(dn_vec) == 1),1); dn(n) = nanmean(dn_vec); + % first method of extracting particulate nitrogen from GLORICH + pn1_vec = tn_vec - dn_vec; pn1_vec(pn1_vec < 0) = 0; + num_pn1(n) = size(find(isfinite(pn1_vec) == 1),1); pn1(n) = nanmean(pn1_vec); + %pn_vec = T{aa,66}; num_pn(n) = size(find(isfinite(pn_vec) == 1),1); pn(n) = nanmean(pn_vec); + tin_vec = T{aa,68}; num_tin(n) = size(find(isfinite(tin_vec) == 1),1); tin(n) = nanmean(tin_vec); + din_vec = T{aa,70}; num_din(n) = size(find(isfinite(din_vec) == 1),1); din(n) = nanmean(din_vec); + ton_vec = T{aa,72}; num_ton(n) = size(find(isfinite(ton_vec) == 1),1); ton(n) = nanmean(ton_vec); + %don_vec = T{aa,74}; num_don(n) = size(find(isfinite(don_vec) == 1),1); don(n) = nanmean(don_vec); + pon_vec = T{aa,76}; num_pon(n) = size(find(isfinite(pon_vec) == 1),1); pon(n) = nanmean(pon_vec); + tkn_vec = T{aa,78}; num_tkn(n) = size(find(isfinite(tkn_vec) == 1),1); tkn(n) = nanmean(tkn_vec); + dkn_vec = T{aa,80}; num_dkn(n) = size(find(isfinite(dkn_vec) == 1),1); dkn(n) = nanmean(dkn_vec); + % second method of extracting particulate nitrogen from GLORICH + pn2_vec = tkn_vec - dkn_vec; pn1_vec(pn2_vec < 0) = 0; + num_pn2(n) = size(find(isfinite(pn2_vec) == 1),1); pn2(n) = nanmean(pn2_vec); + no3_vec = T{aa,82}; num_no3(n) = size(find(isfinite(no3_vec) == 1),1); no3(n) = nanmean(no3_vec); + no2_vec = T{aa,84}; num_no2(n) = size(find(isfinite(no2_vec) == 1),1); no2(n) = nanmean(no2_vec); + no2no3_vec = T{aa,86}; num_no2no3(n) = size(find(isfinite(no2no3_vec) == 1),1); no2no3(n) = nanmean(no2no3_vec); + tnh4_vec = T{aa,88}; num_tnh4(n) = size(find(isfinite(tnh4_vec) == 1),1); tnh4(n) = nanmean(tnh4_vec); + dnh4_vec = T{aa,90}; num_dnh4(n) = size(find(isfinite(dnh4_vec) == 1),1); dnh4(n) = nanmean(dnh4_vec); + don_vec = dn_vec-no2no3_vec-dnh4(n); don_vec(don_vec < 0) = 0; + num_don(n) = size(find(isfinite(don_vec) == 1),1); don(n) = nanmean(don_vec); + % phosphorus + tp_vec = T{aa,92}; num_tp(n) = size(find(isfinite(tp_vec) == 1),1); tp(n) = nanmean(tp_vec); + dp_vec = T{aa,94}; num_dp(n) = size(find(isfinite(dp_vec) == 1),1); dp(n) = nanmean(dp_vec); + pp_vec = tp_vec - dp_vec; pp_vec(pp_vec < 0) = 0; + num_pp(n) = size(find(isfinite(pp_vec) == 1),1); pp(n) = nanmean(pp_vec); + %pp_vec = T{aa,96}; pp(n) = nanmean(pp_vec); + tip_vec = T{aa,98}; num_tip(n) = size(find(isfinite(tip_vec) == 1),1); tip(n) = nanmean(tip_vec); + dip_vec = T{aa,100}; num_dip(n) = size(find(isfinite(dip_vec) == 1),1); dip(n) = nanmean(dip_vec); + dop_vec = dp_vec - dip_vec; dop_vec(dop_vec < 0) = 0; + num_dop(n) = size(find(isfinite(dop_vec) == 1),1); dop(n) = nanmean(dop_vec); + % oxygen + o2_vec = T{aa,12}; num_o2(n) = size(find(isfinite(o2_vec) == 1),1); o2(n) = nanmean(o2_vec); + % silica + si_vec = T{aa,34}; num_si(n) = size(find(isfinite(si_vec) == 1),1); si(n) = nanmean(si_vec); + % dic, fill in with pH if needed + %dic_vec = T{aa,52}; + %bb = find(isfinite(dic_vec)); + pH_vec = T{aa,10}; pH(n) = nanmean(pH_vec); + temp_vec = T{aa,8}; temp(n) = nanmean(temp_vec); + + + % Simple initial DIC calculation based on mean alk and mean pH, room for + % improvement here. + + if isnan(temp(n)); temp(n) = 5; end; + %[RESULT,HEADERS,NICEHEADERS]=CO2SYS(PAR1,PAR2,PAR1TYPE,PAR2TYPE,... + % ...SAL,TEMPIN,TEMPOUT,PRESIN,PRESOUT,SI,PO4,NH4,H2S,pHSCALEIN,... + % ...K1K2CONSTANTS,KSO4CONSTANT,KFCONSTANT,BORON) + out = CO2SYS(alk(n),pH(n),1,3,0,temp(n),temp(n), ... + 100,100,0,0,0,0,4,15,1,2,2); + dic(n) = mean(out(:,2)); +end + +lon_stations_glorich = lon; +lat_stations_glorich = lat; +station_names_glorich = river_name; + +Q_ann_glorich = discharge; +% Set Churchill discharge from GlobalNEWS +Q_ann_glorich(7) = 2281; +dic_ann_glorich = dic; +alk_ann_glorich = alk; +no3_ann_glorich = no2no3; +% use no3 for the Stikine +no3_ann_glorich(2) = no3(1); +nh4_ann_glorich = dnh4; +% assume nh4 in Stikine is equal to nh4 in Skeena +nh4_ann_glorich(2) = nh4_ann_glorich(1); +din_ann_glorich = no3_ann_glorich + nh4_ann_glorich; +% Use weighted average of TN-DN and TKN-DKN +% not that there is no pn data for the Churchill river. Will fill in with +% GlobalNEWS. +pn1(num_pn1 == 0) = 0; pn2(num_pn2 == 0) = 0; +pn_ann_glorich = (pn1.*num_pn1 + pn2.*num_pn2)./(num_pn1 + num_pn2); +% assume particulate nitrogen in Stikine is equal to Skeena +pn_ann_glorich(2) = pn_ann_glorich(1); +don_ann_glorich = don; +% assume dissolved organic nitrogen in Stikine is equal to Skeena +don_ann_glorich(2) = don_ann_glorich(1); + +dip_ann_glorich = dip; +% assume po4 in Stikine is equal to po4 in GLORICH +dip_ann_glorich(2) = dip_ann_glorich(1); + +dop_ann_glorich = dop; +% assume dop in the Stikine and Skeena are equal to that in the Copper +% DOP constraints are generally poor. The Fraser is based on a single +% joint dp, dip measurement. However, the contribution of dop to the total +% phosphorus load is also very low, so this is an issue I'm willing to +% stomach. +dop_ann_glorich(1) = dop_ann_glorich(3); +dop_ann_glorich(2) = dop_ann_glorich(3); + +pp_ann_glorich = pp; +% set particulate load in the Stikine to that in the Skeena +pp_ann_glorich(2) = pp_ann_glorich(1); + +si_ann_glorich = si; +o2_ann_glorich = o2; +dfe_ann_glorich = ones(size(o2_ann_glorich))*NaN; +pfe_ann_glorich = ones(size(o2_ann_glorich))*NaN; + +% did not try and extract monthly data from GLORICH, so leave as NaNs +dic_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +alk_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +no3_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +nh4_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +din_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +don_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +pn_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +dip_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +dop_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +pp_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +dfe_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +pfe_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +si_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +o2_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; + +save Arctic_GLORICH_data lon_stations_glorich lat_stations_glorich station_names_glorich ... + Q_ann_glorich dic_ann_glorich alk_ann_glorich no3_ann_glorich nh4_ann_glorich ... + din_ann_glorich pn_ann_glorich don_ann_glorich dip_ann_glorich dop_ann_glorich ... + pp_ann_glorich si_ann_glorich o2_ann_glorich dfe_ann_glorich pfe_ann_glorich ... + dic_monthly_glorich alk_monthly_glorich no3_monthly_glorich nh4_monthly_glorich ... + din_monthly_glorich don_monthly_glorich pn_monthly_glorich dip_monthly_glorich ... + dop_monthly_glorich pp_monthly_glorich dfe_monthly_glorich pfe_monthly_glorich ... + si_monthly_glorich o2_monthly_glorich + +% Ordering of the GLORICH Data +%1. "STAT_ID", +%2. "RESULT_DATETIME", +%3. "SAMPLE_TIME_DESC", +%4. "SAMPLING_MODE", +%5. "Ref", +%6. "Discharge_inst", +%7. "Discharge_inst_vrc", +%8. "Temp_water", +%9. "Temp_water_vrc", +%10. "pH", +%11. "pH_vrc", +%12. "DO_mgL", +%13. "DO_mgL_vrc", +%14. "DOSAT", +%15. "DOSAT_vrc", +%16. "SpecCond25C", +%17. "SpecCond25C_vrc", +%18. "SPM", +%19. "SPM_vrc", +%20. "Alkalinity", +%21. "Alkalinity_vrc", +%22. "HCO3", +%23. "HCO3_vrc", +%24. "CO3", +%25. "CO3_vrc", +%26. "Ca", +%27. "Ca_vrc", +%28. "Mg", +%29. "Mg_vrc", +%30. "Na", +%31. "Na_vrc", +%32. "K", +%33. "K_vrc", +%34. "SiO2", +%35. "SiO2_vrc", +%36. "Cl", +%37. "Cl_vrc", +%38. "SO4", +%39. "SO4_vrc", +%40. "F", +%41. "F_vrc", +%42. "DSr", +%43. "DSr_vrc", +%44. "TC", +%45. "TC_vrc", +%46. "DC", +%47. "DC_vrc", +%48. "PC", +%49. "PC_vrc", +%50. "TIC", +%51. "TIC_vrc", +%52. "DIC", +%53. "DIC_vrc", +%54. "PIC", +%55. "PIC_vrc", +%56. "TOC", +%57. "TOC_vrc", +%58. "DOC", +%59. "DOC_vrc", +%60. "POC", +%61. "POC_vrc", +%62. "TN", +%63. "TN_vrc", +%64. "DN", +%65. "DN_vrc", +%66. "PN", +%67. "PN_vrc", +%68. "TIN", +%69. "TIN_vrc", +%70. "DIN", +%71. "DIN_vrc", +%72. "TON", +%73. "TON_vrc", +%74. "DON", +%75. "DON_vrc", +%76. "PON", +%77. "PON_vrc", +%78. "TKN", +%79. "TKN_vrc", +%80. "DKN", +%81. "DKN_vrc", +%82. "NO3", +%83. "NO3_vrc", +%84. "NO2", +%85. "NO2_vrc", +%86. "NO2_NO3", +%87. "NO2_NO3_vrc", +%88. "TNH4", +%89. "TNH4_vrc", +%90. "DNH4", +%91. "DNH4_vrc", +%92. "TP", +%93. "TP_vrc", +%94. "DP", +%95. "DP_vrc", +%96. "PP", +%97. "PP_vrc", +%98. "TIP", +%99. "TIP_vrc", +%100. "DIP", +%101. "DIP_vrc", +%102. "PS", +%103. "PS_vrc" diff --git a/tools/rivers/bgc/ARC/README.md b/tools/rivers/bgc/ARC/README.md new file mode 100644 index 000000000..469e88e89 --- /dev/null +++ b/tools/rivers/bgc/ARC/README.md @@ -0,0 +1,112 @@ +## Example Scripts for ARC BGC runoff generation + +This folder contains example scripts for ARC BGC runoff file generation. Users can follow the following instructions to generate BGC runoff file: + +1: Generate a monthly climatology of the river inputs on the model grid +using `make_discharge_climatology.m`. This routine creates a monthly climatology +from the model's freshwater forcing. For example, the default uses GLOFAS/HILL +forcing files created by Kate Hedstrom on 5/11/2023. All you need to do to +regenerate files (or create new ones) is update the file path. The file creates +a matlab file (`*.mat`) with the discharge climatology. This is used in the +assignment of rivers to discharge points. +``` +matlab232 -nodisplay -nosplash -nodesktop -r "run('make_discharge_climatology_arctic.m');exit;" +``` + +2: Use `mapriv_NEWS2.m` to create a river nutrient input file based on the +GlobalNEWS2 estimates (Mayorga et al., 2010). GlobalNEWS contains a +database of global rivers with empirically-derived nutrient inputs. GlobalNEWS +does not, however, have DIC or alkalinity so it cannot be used to provide forcing for +carbon cycling. Also, while globalNEWS is quite skilfull globally, it can +have significant regional biases. Nonetheless, the routine provides a good +way to identify major rivers in the region, and the comprehensive nutrient +estimates that it provides will eventually be used to fill in gaps in river +forcing drawn from direct observations. + +More details about the mapping algorithm can be found below (and in the code). +The first time through, I would recommend setting `inspect_map1` to `y`. +This allows you to step through the mapping of each river and assess its +quality and properties. I have also found that applying a minimum discharge +of 100 m3 sec-1 helps avoid erroneous mapping of very small rivers onto +large discharges. The assignment algorithm was designed to be relatively +insensitive to such instances, but care should still be taken. + +The required inputs are the discharge climatology (from step 1) and a copy +of the globalNEWS data. +``` +matlab232 -nodisplay -nosplash -nodesktop -r "run('mapriv_NEWS2.m');exit;" +``` + +3: Gather/Process direct river measurements: The next step is to get as many +direct measurements as you can. + + - [RC4USCoast](https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.nodc:0260455) +``` +mkdir -p Data/RC4USCoast +cd Data/RC4USCoast +wget https://www.ncei.noaa.gov/archive/archive-management-system/OAS/bin/prd/jquery/download/260455.3.3.tar.gz +tar -xzf 260455.3.3.tar.gz --wildcards --strip-components=4 -C . "*/data/0-data/*.nc" +``` + + - [GLORICH](https://www.geo.uni-hamburg.de/en/geologie/forschung/aquatische-geochemie/glorich.html) +``` +cd Data/GLORICH +wget https://store.pangaea.de/Publications/HartmannJens-etal_2019/Glorich_V01_CSV_plus_Shapefiles_2019_05_24.zip +unzip Glorich_V01_CSV_plus_Shapefiles_2019_05_24.zip +matlab232 -nodisplay -nosplash -nodesktop -r "run('Arctic_GLORICH_Process.m');exit;" +``` + + - [ArcticGro](https://arcticgreatrivers.org/data/) +``` +# A copy of this dataset is provided in the Data directory. +# You may want to download your own to ensure +# that it is up-to-date. +cd Data/ArcticGro +matlab232 -nodisplay -nosplash -nodesktop -r "run('ArcticGro_Process.m');exit;" +``` +4: run `mapriv_combined_Arctic` to create river nutrient and carbon input +estimates based on available observation, while using GlobalNEWS to fill in +some gaps. +``` +cp /archive/ynt//woa_sst_climo.nc ./Data/ +matlab232 -nodisplay -nosplash -nodesktop -r "run('mapriv_combined_Arctic.m');exit;" +``` +As was the case for globalNEWS2, I recommend running these with "inspect_map = 'y'" until +you are satisfied with the results. You can just "click through" each river mapping and +confirm its properties. The routine will then produce numerous final plots for analysis +and quality control, before generating the netcdf for use with MOM6. + +Note: River oxygen levels are set at saturating levels at temperatures provided by the +World Ocean Atlas Climatology. +___________________________________________________________________________________________ + +## RIVER MAPPING ALGORITHM: +Rivers are first filtered to isolate those within the domain and +above a user specified flow threshold. The rationale for this threshold is to minimize the +risk of inadvertantly mapping very small rivers onto large freshwater flows. Small rivers +tend to have more volatile nutrient concentrations, so such mistakes can have large consequences. +To further reduce this risk, rivers are then sorted by size from smallest to largest. Model +outflow points nearest to each river are accumulated until the the value closest to the observed +flow is reached. These points are assigned the concentration for that river. If a larger river +subsequently claims those points, the larger river is given precedence. The concentrations of +any model discharge points left unassigned after all estimates have been cycled through +are assigned using a nearest neighbor algorithm. + +## Note: +One must specify the fraction of particulate phosphorus that is bioavailable (generally between +10-30%, Froelich, Kinetic control of dissolved phosphate in natural rivers and estuaries: A primer +on the phosphate buffer mechanism, Limnology and Oceanography), and the fraction of dissolved organic +inputs that fall into labile, semi-labile and semi-refractory pools. This is set by default to 30%, 35%, +and 35%. This is consistent with the range of bio-availability suggested by Weigner et al. (2006), +Bio-availability of dissolved organic nitrogen and carbon from 9 rivers in the eastern US, Aquatic +Microbial Ecology. + +## VISUALIZATION AND QUALITY CONTROL: +The routines include a y/n toggle option for inspecting the +mapping of each river as it is done. If "inspect_map" is set to 'y', a graphical map of the +model discharge point assigned to each river is presented, along with the outflow and nutrient +characteristics of the river. This can be useful for identifying rivers that may require some +manual editing to ensure the outflow is assigned to the right place. A section for such manual +edits is provided in the code. The visualization pauses until the user presses any key. It +then moves onto the next river. Once all the rivers have been mapped and interpolation completed, +the routine always produces a series of domain-wide plots to evaluate the overall results. diff --git a/tools/rivers/bgc/ARC/make_discharge_climatology_arctic.m b/tools/rivers/bgc/ARC/make_discharge_climatology_arctic.m new file mode 100644 index 000000000..1223376d9 --- /dev/null +++ b/tools/rivers/bgc/ARC/make_discharge_climatology_arctic.m @@ -0,0 +1,61 @@ +clear all; +addpath /home/cas/matlab +nc64startup; + +syear = 1993; +eyear = 2019; +num_years = eyear-syear+1; +readme = '1993-2019 monthly Arctic climatology from GLOFAS/Hill (May 11, 2023 from Kate Hedstrom), kg m-2 sec-1'; + +% get grid information +file_name = ['/archive/cas/Regional_MOM6/Arctic/glofas_hill_05112023/glofas_hill_',num2str(syear,'%4u'),'.nc'] +lon = ncread(file_name,'lon'); +lat = ncread(file_name,'lat'); +nlat = size(lat,2); +nlon = size(lon,1); +area = ncread(file_name,'area'); + +% holds the runoff climatology +runoff_mc = zeros(nlon,nlat,12); + +for yr = syear:eyear + file_name = ['/archive/cas/Regional_MOM6/Arctic/glofas_hill_05112023/glofas_hill_',num2str(yr,'%4u'),'.nc'] + + time = ncread(file_name,'time'); + ntime = size(time,1); + + date = datevec(time+datenum(1950,1,1,0,0,0)); + runoff_days_temp = ncread(file_name,'runoff'); + % files are padded with the first day of the following year, remove + runoff_days = runoff_days_temp(:,:,1:(ntime-1)); + month = date(1:(ntime-1),2); + + for m = 1:12 + aa = find(month == m); + runoff_temp = runoff_days(:,:,aa); + runoff_mc(:,:,m) = runoff_mc(:,:,m) + mean(runoff_temp,3)/num_years; + end + + clear runoff_days runoff_days_temp; + clear runoff_temp aa time date month; + +end + +% modify so that it is time, nlat, nlon +runoff = permute(runoff_mc,[3 2 1]); +area = permute(area,[2 1]); +lon = permute(lon,[2 1]); +lat = permute(lat,[2 1]); + +for m = 1:12; temp = squeeze(runoff(m,:,:)); total_runoff(m) = sum(temp(:).*area(:)); end +figure(1); clf; plot(total_runoff); xlabel('month'); ylabel('runoff, kg sec-1'); + +for m = 1:12; + figure(2); + clf + temp = squeeze(runoff(m,:,:)); + scatter3(lon(:),lat(:),log10(temp(:)),ones(size(temp(:)))*10,log10(temp(:))); view(2); caxis([-3 -1]); colorbar; + pause +end + +save glofas_hill_runoff_monthlyclim_arctic12k_05112023 runoff area lat lon readme; diff --git a/tools/rivers/bgc/ARC/mapriv_NEWS2.m b/tools/rivers/bgc/ARC/mapriv_NEWS2.m new file mode 100644 index 000000000..53542f4c8 --- /dev/null +++ b/tools/rivers/bgc/ARC/mapriv_NEWS2.m @@ -0,0 +1,859 @@ +% Routine to map Global NEWS nutrient data onto the MOM6 Arctic grid @ + +clear all; +addpath /home/cas/matlab +nc64startup; + +% name of netcdf file to be created +%nc_file_name = 'RiverNutrients_GlobalNEWS2_plusFe_Q100_NEP10k.nc'; +% Arctic +nc_file_name = 'RiverNutrients_GlobalNEWS2_plusFe_Q100_Arctic12k.nc'; + +% Parameters for the assignment algorithm. +Q_min = 100; % minimum flow in m3 sec +plot_width = 15; % width of window (in degrees) for inspecting locations + % of rivers and outflow points that have been assigned to + % them. +min_dist = 2.5; % minimum distance (degrees) of the closest outflow point + % for the river to be considered in the domain (useful + % for preventing the algorithm from trying to map rivers + % flowing to different ocean basins. +max_dist = 3; % maximum distance to search for a point +inspect_map = 'n'; % flag enabling you to pause and inspect each river + % mapping as it is being done. + +% set the bio-availability of phosphorus and the fractionation of dissolved +% organic; The code assumes that 30% of particulate phosphorus (PP) is +% bioavailable (e.g., Frolich et al., 1988). One also needs to divide the +% dissolved organic nitrogen and phosphorus components into fractions with +% different lability. The default values are based losely on Wiegner et +% al., (2006). If you have better information, feel free to modify. +frac_PP = 0.3; +frac_ldon = 0.3; +frac_sldon = 0.35; +frac_srdon = 0.35; +frac_ldop = 0.3; +frac_sldop = 0.35; +frac_srdop = 0.35; +% 40 nM dissolved iron concentration from De Baar and De Jong + 30nM +% Colloidal and nanoparticle flux as reported in Canfield and Raiswell +const_fed = 70.0e-6; + +% GlobalNEWS2 data obtained from Emilio Mayorga +filename = '/archive/cas/COBALT_EVAL/River_Data/GlobalNEWS2/GlobalNEWS2_RH2000Dataset-version1.0.xls' +basin = readtable(filename,'Sheet',2); +hydrology = readtable(filename,'Sheet',3); +load = readtable(filename,'Sheet',4); + +% find all the river basins that empty into "land", e.g., lakes +ocean = basin.ocean; +land_index = zeros(size(ocean)); +for n = 1:size(ocean,1); + test = strcmp('Land',ocean(n)); + land_index(n) = single(test); +end + +river_names_all = basin.basinname; + +% basin area in +area = basin.A; +lon_news_all = basin.mouth_lon; +lat_news_all = basin.mouth_lat; +% Depending on the model output, may need to adjust longitudes +lon_news_all(lon_news_all < 0) = lon_news_all(lon_news_all < 0) + 360; + +% Loads in Mg yr-1 converted to moles per sec +DIN_load_all = load.Ld_DIN*1e6/14/86400/365; +DIP_load_all = load.Ld_DIP*1e6/31/86400/365; +DON_load_all = load.Ld_DON*1e6/14/86400/365; +DOP_load_all = load.Ld_DOP*1e6/31/86400/365; +Si_load_all = load.Ld_DSi*1e6/28.1/86400/365; +PN_load_all = (load.Ld_PN*1e6/14/86400/365); +PP_load_all = (load.Ld_PP*1e6/31/86400/365)*frac_PP; + +% actual and natural discharge (convert from km3/yr to m3/sec) +% Used the actual hydrology to calculate concentrations +Qact_all = hydrology.Qact*1e9/(86400*365); +Qnat_all = hydrology.Qnat*1e9/(86400*365); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Load in monthly climatology of river forcing from the regional grid. % +% File contains: % +% runoff: monthly average runoff in kg m-2 sec-1 % +% area: area of grid cell in m-2 % +% lon: longitude (0-360 degrees) % +% lat: latitude % +% % +% The file was calculated from daily output using the routine % +% "make_climatology.m". This routine and all associated files can be % +% found in the same directory as the data file % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% NEP +%load glofas_hill_runoff_monthlyclim_NEP10k_05122023.mat; +% Arctic +load glofas_hill_runoff_monthlyclim_arctic12k_05112023.mat; +lon_mod = lon; +lat_mod = lat; +area_mod = area; +% convert runoff from kg m-2 sec-1 to m3 sec-1 +Q_mod_monthly = zeros(size(runoff)); +for m = 1:12 + Q_mod_monthly(m,:,:) = squeeze(runoff(m,:,:)).*area_mod./1000; +end +Q_mod_ann = squeeze(mean(Q_mod_monthly,1)); +clear lon lat runoff area; + +%grid_file = '/archive/cas/Regional_MOM6/NWA12/nwa12_ocean_static.nc'; +%temp = ncread(grid_file,'deptho'); +%depth = permute(temp,[2,1]); clear temp; +%depth(isnan(depth)) = -1; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Filter for rivers in the region, set thresholds for minimum river size, % +% set parameters for plotting routines. % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% use grid to filter rivers outside domain +lat_mod_max = max(lat_mod(:)); +lat_mod_min = min(lat_mod(:)); +lon_mod_max = max(lon_mod(:)); +lon_mod_min = min(lon_mod(:)); + +in_region = find( (lon_news_all <= lon_mod_max) & (lon_news_all >= lon_mod_min) & ... + (lat_news_all <= lat_mod_max) & (lat_news_all >= lat_mod_min) & ... + (isfinite(Qact_all)) & (Qact_all > Q_min) ); + +% If you are using a high threshold, grab one smaller river to constrain +% Carribean Islands +%if Q_min > 100 +% for n = 1:size(lon_news_all,1); +% test = strcmp('GHAASBasin1808',river_names_all{n}); +% if test == 1; +% cuba_ind = n; +% end +% end +% in_region = [in_region; cuba_ind]; +%end + +num_rivers = size(in_region,1); + +% Establish vectors of flow and nutrient loads for the NWA +Qact = Qact_all(in_region); +lon_news = lon_news_all(in_region); +lat_news = lat_news_all(in_region); +DIN_load = DIN_load_all(in_region); +DON_load = DON_load_all(in_region); +PN_load = PN_load_all(in_region); +DIP_load = DIP_load_all(in_region); +DOP_load = DOP_load_all(in_region); +PP_load = PP_load_all(in_region); +Si_load = Si_load_all(in_region); +river_names = river_names_all(in_region); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Following inspection of initial mapping, add any manual edits here to % +% prevent anomalous extrapolations, etc. % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %Add in 2 "dummy" rivers to handle Cuba, giving it properties like + %Jamaica or Haiti rather than Florida. + + %GHAASBasin1808 is in Haiti. Fluxes are characterized by particularly + %high particulate phosphorus inputs. +%for n = 1:num_rivers; +% n +% test = strcmp('GHAASBasin1808',river_names{n}) +% if test == 1; +% cuba_ind = n; +% end + + % Move the Susquehanna a bit south so that it catches the Chesapeake + % and not the Delaware. +% if strcmp('Susquehanna',river_names{n}) +% lat_news(n) = 38.5; +% lon_news(n) = -76.67; +% end +%end + +% Two "rivers" with properties like Haiti used to specify Cuba +%Qact(num_rivers+1) = Qact(cuba_ind)/2; +%lon_news(num_rivers+1) = -81; +%lat_news(num_rivers+1) = 22.6; +%DIN_load(num_rivers+1) = DIN_load(cuba_ind)/2; +%DON_load(num_rivers+1) = DON_load(cuba_ind)/2; +%PN_load(num_rivers+1) = PN_load(cuba_ind)/2; +%DIP_load(num_rivers+1) = DIP_load(cuba_ind)/2; +%DOP_load(num_rivers+1) = DOP_load(cuba_ind)/2; +%PP_load(num_rivers+1) = PP_load(cuba_ind)/2; +%Si_load(num_rivers+1) = Si_load(cuba_ind)/2; +%river_names(num_rivers+1) = river_names(cuba_ind); + +%Qact(num_rivers+2) = Qact(cuba_ind)/2; +%lon_news(num_rivers+2) = -83.25; +%lat_news(num_rivers+2) = 22.6; +%DIN_load(num_rivers+2) = DIN_load(cuba_ind)/2; +%DON_load(num_rivers+2) = DON_load(cuba_ind)/2; +%PN_load(num_rivers+2) = PN_load(cuba_ind)/2; +%DIP_load(num_rivers+2) = DIP_load(cuba_ind)/2; +%DOP_load(num_rivers+2) = DOP_load(cuba_ind)/2; +%PP_load(num_rivers+2) = PP_load(cuba_ind)/2; +%Si_load(num_rivers+2) = Si_load(cuba_ind)/2; +%river_names(num_rivers+2) = river_names(cuba_ind); + +%num_rivers = num_rivers + 1; + +% Adjust location of cfilename = '/archive/cas/COBALT_EVAL/River_Data/GlobalNEWS2/GlobalNEWS2_RH2000Dataset-version1.0.xls' +basin = readtable(filename,'Sheet',2); +hydrology = readtable(filename,'Sheet',3); +load = readtable(filename,'Sheet',4); +% Cuba; This has little effect on patterns in Florida. +%for n = 1:num_rivers; +% n +% test = strcmp('GHAASBasin448',river_names{n}) +% if test == 1; +% fla_ind = n; +% end +%end + +%lon_news(fla_ind) = -80.5; +%lat_news(fla_ind) = 26.6; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% END MANUAL EDITS % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Assigning outflow points to rivers. % +% 1. Assignment starts with the rivers with the smallest flow and works % +% to the largest, w/larger river characteristics taking precedence to % +% ensure the most significant rivers are well represented. % +% 2. The algorithm keeps choosing the closest points to each river mouth % +% until the assigned flow is as close as possible to that observed % +% 3. Once the outflow points are assigned using the mean flow values, % +% monthly concentrations are assigned to those points. % +% 4. A simple "nearest neighbor" algorithm is used to fill in the gaps % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Sort rivers by discharge +[Qact_sort,sort_ind] = sort(Qact,'ascend'); +lon_news_sort = lon_news(sort_ind); +lat_news_sort = lat_news(sort_ind); +DIN_load_sort = DIN_load(sort_ind); +DON_load_sort = DON_load(sort_ind); +PN_load_sort = PN_load(sort_ind); +DIP_load_sort = DIP_load(sort_ind); +DOP_load_sort = DOP_load(sort_ind); +PP_load_sort = PP_load(sort_ind); +Si_load_sort = Si_load(sort_ind); +river_names_sort = river_names(sort_ind); + +% Total N and P load diagnostics +N_load_sort = DIN_load_sort + DON_load_sort + PN_load_sort; +P_load_sort = DIP_load_sort + DOP_load_sort + PP_load_sort; + +% Calculate concentrations +% Loads are in moles N sec-1, Q in m3 s-1; conc in moles N m-3 +DIN_conc_sort = DIN_load_sort./Qact_sort; +DON_conc_sort = DON_load_sort./Qact_sort; +DIP_conc_sort = DIP_load_sort./Qact_sort; +DOP_conc_sort = DOP_load_sort./Qact_sort; +PN_conc_sort = PN_load_sort./Qact_sort; +PP_conc_sort = PP_load_sort./Qact_sort; +Si_conc_sort = Si_load_sort./Qact_sort; + +% initialize vectors to hold nutrient concentrations at eac runoff +% point. +aa = find(Q_mod_ann > 0); +Q_mod_vec = Q_mod_ann(aa); +din_conc_vec = zeros(size(Q_mod_vec)); +don_conc_vec = zeros(size(Q_mod_vec)); +pn_conc_vec = zeros(size(Q_mod_vec)); +dip_conc_vec = zeros(size(Q_mod_vec)); +dop_conc_vec = zeros(size(Q_mod_vec)); +pp_conc_vec = zeros(size(Q_mod_vec)); +si_conc_vec = zeros(size(Q_mod_vec)); + +lon_mod_runoff_vec = double(lon_mod(aa)); +lat_mod_runoff_vec = double(lat_mod(aa)); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Loop identifies points assigned to each river % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +for k=1:num_rivers + k + dist = pdist2([lon_news_sort(k) lat_news_sort(k)], ... + [lon_mod_runoff_vec lat_mod_runoff_vec]); + [dist_sort, dist_sort_ind] = sort(dist,'ascend'); + + if dist_sort(1) < min_dist; % filters out rivers lying outside the domain + Q_sum1 = 0; + Q_sum2 = 0; + n = 0; + while ((Q_sum2 < Qact_sort(k)) & (dist_sort(n+1) < max_dist)) + Q_sum1 = Q_sum2; + n = n+1; + Q_sum2 = Q_sum1 + Q_mod_vec(dist_sort_ind(n)); + end + % I generally find that the search algorithm works best when you keep + % grabbing points until the total flow captured exceeds the flow in + % the river. If you uncomment the the first part of the "if" + % statement here will pick the last below if it is closer than the + % first above. However, I've found that this option sometimes fails + % to map important rivers. + %if abs(Q_sum1 - Qact_sort(k)) < abs(Q_sum2 - Qact_sort(k)) + % nrp = n-1; % number of runoff points + % [Q_sum1 Qact_sort(k)] % a quick check for comparable flow + % din_conc_vec(dist_sort_ind(1:nrp)) = DIN_conc_sort(k); + % don_conc_vec(dist_sort_ind(1:nrp)) = DON_conc_sort(k); + % dip_conc_vec(dist_sort_ind(1:nrp)) = DIP_conc_sort(k); + % dop_conc_vec(dist_sort_ind(1:nrp)) = DOP_conc_sort(k); + % pn_conc_vec(dist_sort_ind(1:nrp)) = PN_conc_sort(k); + % pp_conc_vec(dist_sort_ind(1:nrp)) = PP_conc_sort(k); + % si_conc_vec(dist_sort_ind(1:nrp)) = Si_conc_sort(k); + %else + nrp = n; % number of runoff points + [Q_sum2 Qact_sort(k)] % a quick check for comparable flow + din_conc_vec(dist_sort_ind(1:nrp)) = DIN_conc_sort(k); + don_conc_vec(dist_sort_ind(1:nrp)) = DON_conc_sort(k); + dip_conc_vec(dist_sort_ind(1:nrp)) = DIP_conc_sort(k); + dop_conc_vec(dist_sort_ind(1:nrp)) = DOP_conc_sort(k); + pn_conc_vec(dist_sort_ind(1:nrp)) = PN_conc_sort(k); + pp_conc_vec(dist_sort_ind(1:nrp)) = PP_conc_sort(k); + si_conc_vec(dist_sort_ind(1:nrp)) = Si_conc_sort(k); + %end + + if inspect_map == 'y' + figure(1) + clf + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(Q_mod_vec),3,log10(Q_mod_vec)); + hold on + scatter3(lon_mod_runoff_vec(dist_sort_ind(1:nrp)),lat_mod_runoff_vec(dist_sort_ind(1:nrp)), ... + log10(Q_mod_vec(dist_sort_ind(1:nrp))),40, ... + log10(Q_mod_vec(dist_sort_ind(1:nrp))),'filled'); + view(2); + plot3(lon_news_sort(k),lat_news_sort(k),1e5,'k.','MarkerSize',20); + %contour(lon_mod,lat_mod,depth,[0 0],'k-'); + axis([lon_news_sort(k)-plot_width lon_news_sort(k)+plot_width ... + lat_news_sort(k)-plot_width lat_news_sort(k)+plot_width]); + caxis([-4 3]); + titl = ['river number: ',num2str(k),' name: ',river_names_sort{k}]; + title(titl); + colorbar; + + % provide a few diagnostics to ensure the calculation was done + % correctly (remove semicolon to inspect as they are mapped in + % the matlab output line. Feel free to add more here. + N_conc = DIN_conc_sort(k) + DON_conc_sort(k) + PN_conc_sort(k); + P_conc = DIP_conc_sort(k) + DOP_conc_sort(k) + PP_conc_sort(k); + Si_conc = Si_conc_sort(k); + river_names_sort(k) + [lon_news_sort(k) lat_news_sort(k)] + ind = dist_sort_ind(1:nrp); + 'lon lat' + [lon_news_sort(k) lat_news_sort(k)] + 'total flow in m3 sec' + [Qact_sort(k) sum(Q_mod_vec(ind))] + 'Nitrogen and phosphorus in Gg per year'; + [N_load_sort(k) sum(Q_mod_vec(ind))*N_conc]*14*86400*365/1e9; + [P_load_sort(k) sum(Q_mod_vec(ind))*P_conc]*31*86400*365/1e9; + 'N, P conc (mmoles m-3), DI, DO, P' + [DIN_conc_sort(k) DON_conc_sort(k) PN_conc_sort(k)]*1e3 + [DIP_conc_sort(k) DOP_conc_sort(k) PP_conc_sort(k)]*1e3 + 'Total N, Total P, Total N: Total P'; + [N_conc*1e3 P_conc*1e3 N_conc/P_conc] + 'DO:DI and P:DI ratios'; + [DON_conc_sort(k)/DIN_conc_sort(k) PN_conc_sort(k)/DIN_conc_sort(k)]; + [DOP_conc_sort(k)/DIP_conc_sort(k) PP_conc_sort(k)/DIP_conc_sort(k)]; + 'silica concentration (mmoles m-3)'; + [Si_conc_sort(k)]*1e3; + pause + end + + else + % This is for rivers that were captured by the coarse initial filter + % to determine if they were in the domain; but are actually outside + % the domain. Calibration suggested that a threshold of 0.75 degrees + % from the specified river mouth reliably identified these cases. + + if inspect_map == 'y' + figure(1) + clf + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(Q_mod_vec),3,log10(Q_mod_vec)); + hold on + view(2); + plot3(lon_news_sort(k),lat_news_sort(k),1e5,'k.','MarkerSize',20); + axis([lon_news_sort(k)-15 lon_news_sort(k)+15 ... + lat_news_sort(k)-15 lat_news_sort(k)+15]); + caxis([-4 3]); + titl = ['OUTSIDE: river number: ',num2str(k),' name: ',river_names_sort{k}]; + title(titl); + colorbar; + pause + end + + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% nearest neighbor search to fill in any 0 values left for each input field +% after the river mapping is done. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +aa = find(din_conc_vec == 0); +bb = find(din_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),din_conc_vec(bb), ... + 'nearest','nearest'); +din_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(don_conc_vec == 0); +bb = find(don_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),don_conc_vec(bb), ... + 'nearest','nearest'); +don_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(pn_conc_vec == 0); +bb = find(pn_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),pn_conc_vec(bb), ... + 'nearest','nearest'); +pn_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(dip_conc_vec == 0); +bb = find(dip_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),dip_conc_vec(bb), ... + 'nearest','nearest'); +dip_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(dop_conc_vec == 0); +bb = find(dop_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),dop_conc_vec(bb), ... + 'nearest','nearest'); +dop_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(pp_conc_vec == 0); +bb = find(pp_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),pp_conc_vec(bb), ... + 'nearest','nearest'); +pp_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(si_conc_vec == 0); +bb = find(si_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),si_conc_vec(bb), ... + 'nearest','nearest'); +si_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +totn_conc_vec = din_conc_vec + don_conc_vec + pn_conc_vec; +totp_conc_vec = dip_conc_vec + dop_conc_vec + pp_conc_vec; + +% calculate ratios of dissolved and particulate to inorganic +din_flux_vec = din_conc_vec.*Q_mod_vec; +dip_flux_vec = dip_conc_vec.*Q_mod_vec; +don_flux_vec = don_conc_vec.*Q_mod_vec; +dop_flux_vec = dop_conc_vec.*Q_mod_vec; +pn_flux_vec = pn_conc_vec.*Q_mod_vec; +pp_flux_vec = pp_conc_vec.*Q_mod_vec; + +don_ratio = sum(don_flux_vec)/sum(din_flux_vec) +dop_ratio = sum(dop_flux_vec)/sum(dip_flux_vec) +pn_ratio = sum(pn_flux_vec)/sum(din_flux_vec) +pp_ratio = sum(pp_flux_vec)/sum(dip_flux_vec) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Produce plots to evaluate the mapping % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Use total fluxes to scale dots, m3 sec-1 * moles m-3 = moles sec-1 +totn_flux_vec = totn_conc_vec.*Q_mod_vec; +totp_flux_vec = totp_conc_vec.*Q_mod_vec; +totsi_flux_vec = si_conc_vec.*Q_mod_vec; + +% scale marker size with the total nitrogen flux +ms_vec = zeros(size(Q_mod_vec)); +ms_vec(log10(Q_mod_vec) < 0) = 1; +ms_vec((log10(Q_mod_vec) > 0) & (log10(Q_mod_vec) < 1)) = 1; +ms_vec((log10(Q_mod_vec) > 1) & (log10(Q_mod_vec) < 2)) = 1; +ms_vec((log10(Q_mod_vec) > 2) & (log10(Q_mod_vec) < 3)) = 1; +ms_vec(log10(Q_mod_vec) > 3) = 100; + +figure(1) +clf + +subplot(1,3,1); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,totn_conc_vec*1e3,ms_vec, ... + totn_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([30 150]); +colorbar +title('total nitrogen concentration, mmoles m-3'); + +subplot(1,3,2); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,totp_conc_vec*1e3,ms_vec, ... + totp_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 10]); +colorbar +title('total phosphorus concentration, mmoles m-3'); + +subplot(1,3,3) +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,totn_conc_vec./totp_conc_vec,ms_vec, ... + totn_conc_vec./totp_conc_vec,'filled'); +hold on +view(2); +caxis([10 100]); +colorbar +title('N:P ratio'); + +figure(2) +clf + +subplot(1,3,1); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,din_conc_vec*1e3,ms_vec, ... + din_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 100]); +colorbar +title('din concentration, mmoles m-3'); + +subplot(1,3,2); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,don_conc_vec*1e3,ms_vec, ... + don_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 100]); +colorbar +title('don concentration, mmoles m-3'); + +subplot(1,3,3) +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pn_conc_vec*1e3,ms_vec, ... + pn_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 100]); +colorbar +title('pn concentration, mmoles m-3'); + +figure(3) +clf + +subplot(1,3,1); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dip_conc_vec*1e3,ms_vec, ... + dip_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 5]); +colorbar +title('dip concentration, mmoles m-3'); + +subplot(1,3,2); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dop_conc_vec*1e3,ms_vec, ... + dop_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 3]); +colorbar +title('dop concentration, mmoles m-3'); + +subplot(1,3,3) +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pp_conc_vec*1e3,ms_vec, ... + pp_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 3]); +colorbar +title('pp concentration, mmoles m-3'); + +figure(10) +clf + +subplot(2,2,1); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,don_conc_vec./din_conc_vec,ms_vec, ... + don_conc_vec./din_conc_vec,'filled'); +hold on +view(2); +caxis([0 2]); +colorbar +title('don:din'); + +subplot(2,2,2); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pn_conc_vec./din_conc_vec,ms_vec, ... + pn_conc_vec./din_conc_vec,'filled'); +hold on +view(2); +caxis([0 2]); +colorbar +title('pn:din'); + +subplot(2,2,3) +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dop_conc_vec./dip_conc_vec,ms_vec, ... + dop_conc_vec./dip_conc_vec,'filled'); +hold on +view(2); +caxis([0 2]); +colorbar +title('dop:dip'); + +subplot(2,2,4) +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pp_conc_vec./dip_conc_vec,ms_vec, ... + pp_conc_vec./dip_conc_vec,'filled'); +hold on +view(2); +caxis([0 2]); +colorbar +title('pp:dip'); + +figure(5) +clf +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,si_conc_vec*1e3,ms_vec, ... + si_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([30 300]); +colorbar +title('si concentration, mmoles m-3'); + +figure(6) +clf + +subplot(1,2,1); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(totn_flux_vec),ms_vec, ... + log10(totn_flux_vec),'filled'); +hold on +view(2); +caxis([-1 2]); +colorbar +title('N flux, log10(moles sec-1)'); + +subplot(1,2,2); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(totp_flux_vec),ms_vec, ... + log10(totp_flux_vec),'filled'); +hold on +view(2); +caxis([-3 1]); +colorbar +title('P flux, log10(moles sec-1)'); + +% output total fluxes in gigagrams of N +total_annual_n = sum(totn_flux_vec)*86400*365*14/1e12 +total_annual_p = sum(totp_flux_vec)*86400*365*31/1e12 +total_annual_si = sum(totsi_flux_vec)*86400*365*28.1/1e12 + +'press any key to continue' +pause + +% Initialize 2D concentration arrays; these are the ones read into MOM6 to +% specify the nutrient concentrations of river inputs. +din_conc = zeros(size(lon_mod)); +don_conc = zeros(size(lon_mod)); +dip_conc = zeros(size(lon_mod)); +dop_conc = zeros(size(lon_mod)); +pn_conc = zeros(size(lon_mod)); +pp_conc = zeros(size(lon_mod)); +si_conc = zeros(size(lon_mod)); + +% Map concentration vectors onto 2D arrays. +aa = find(Q_mod_ann > 0); +din_conc(aa) = din_conc_vec; +don_conc(aa) = don_conc_vec; +pn_conc(aa) = pn_conc_vec; +dip_conc(aa) = dip_conc_vec; +dop_conc(aa) = dop_conc_vec; +pp_conc(aa) = pp_conc_vec; +si_conc(aa) = si_conc_vec; + +NO3_CONC = din_conc; +LDON_CONC = frac_ldon*don_conc; +SLDON_CONC = frac_sldon*don_conc; +SRDON_CONC = frac_srdon*don_conc; +PO4_CONC = dip_conc; +LDOP_CONC = frac_ldop*dop_conc; +SLDOP_CONC = frac_sldop*dop_conc; +SRDOP_CONC = frac_srdop*dop_conc; +NDET_CONC = pn_conc; +PDET_CONC = pp_conc; % The bioavailability of particulate P has already + % been accounted for. + +SI_CONC = si_conc; + +% Add iron concentrations - initialize with nitrate and then overwrite +FED_CONC = NO3_CONC; +FEDET_CONC = NO3_CONC; +% 40 nM dissolved iron concentration from De Baar and De Jong + 30nM +% Colloidal and nanoparticle flux as reported in Canfield and Raiswell +FED_CONC(FED_CONC > 0) = const_fed; +FEDET_CONC(FEDET_CONC > 0) = 0.0; + +% series of quick figures to check the 2D nutrient input files. +ms = 8; +figure(7); +clf +subplot(3,2,1); +title('log10(NO3 CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(NO3_CONC(:)),ms,log10(NO3_CONC(:)),'filled'); +caxis([-4 -1]); colorbar; + +subplot(3,2,2); +title('log10(LDON CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(LDON_CONC(:)),ms,log10(LDON_CONC(:)),'filled'); +caxis([-4 -1]); colorbar; + +subplot(3,2,3); +title('log10(SLDON CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(SLDON_CONC(:)),ms,log10(SLDON_CONC(:)),'filled'); +caxis([-4 -1]); colorbar; + +subplot(3,2,4); +title('log10(SRDON CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(SRDON_CONC(:)),ms,log10(SRDON_CONC(:)),'filled'); +caxis([-4 -1]); colorbar; + +subplot(3,2,5); +title('log10(NDET CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(NDET_CONC(:)),ms,log10(NDET_CONC(:)),'filled'); +caxis([-4 -1]); colorbar; + +figure(8); +clf +subplot(3,2,1); +title('log10(PO4 CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(PO4_CONC(:)),ms,log10(PO4_CONC(:)),'filled'); +caxis([-4 -2]); colorbar; + +subplot(3,2,2); +title('log10(LDOP CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(LDOP_CONC(:)),ms,log10(LDOP_CONC(:)),'filled'); +caxis([-4 -2]); colorbar; + +subplot(3,2,3); +title('log10(SLDOP CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(SLDOP_CONC(:)),ms,log10(SLDOP_CONC(:)),'filled'); +caxis([-4 -2]); colorbar; + +subplot(3,2,4); +title('log10(SRDOP CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(SRDOP_CONC(:)),ms,log10(SRDOP_CONC(:)),'filled'); +caxis([-4 -2]); colorbar; + +subplot(3,2,5); +title('log10(PDET CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(PDET_CONC(:)),ms,log10(PDET_CONC(:)),'filled'); +caxis([-4 -2]); colorbar; + +figure(9) +clf +clf +subplot(3,2,1); +title('log10(FED CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(FED_CONC(:)),ms,log10(FED_CONC(:)),'filled'); +caxis([-5 -3]); colorbar; + +subplot(3,2,2); +title('log10(FEDET CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(FEDET_CONC(:)),ms,log10(FEDET_CONC(:)),'filled'); +caxis([-5 -3]); colorbar; + +subplot(3,2,3); +title('log10(SI CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(SI_CONC(:)),ms,log10(SI_CONC(:)),'filled'); +caxis([-3 -0]); colorbar; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Save Files % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% option to save matlab file +% save River_DIC_ALK_USGS_NWA ALK_CONC DIC_CONC + +% Construct netcdf file following format used by nutrient input files to +% MOM6 +time = 0; +nlat = size(lat_mod,1); +nlon = size(lat_mod,2); + +ncid = netcdf.create(nc_file_name,'CLOBBER'); +dimid0 = netcdf.defDim(ncid,'time',netcdf.getConstant('NC_UNLIMITED')); + +dimid1 = netcdf.defDim(ncid,'y',nlat); +dimid2 = netcdf.defDim(ncid,'x',nlon); + +%attributes inherited from the old river nutrient forcing file. The +%calendar needs to be specified even though the nutrient values are +%constant. Not sure how much of the rest is needed. +varid0 = netcdf.defVar(ncid,'time','double',dimid0); +netcdf.putAtt(ncid,varid0,'calendar','NOLEAP'); +netcdf.putAtt(ncid,varid0,'calendar_type','NOLEAP'); +netcdf.putAtt(ncid,varid0,'modulo','T'); +netcdf.putAtt(ncid,varid0,'units','days since 1900-1-1 0:00:00'); +netcdf.putAtt(ncid,varid0,'time_origin','01-JAN-1990 00:00:00'); + +varid1 = netcdf.defVar(ncid,'y','int',dimid1); +netcdf.putAtt(ncid,varid1,'cartesian_axis','Y'); +varid2 = netcdf.defVar(ncid,'x','int',dimid2); +netcdf.putAtt(ncid,varid2,'cartesian_axis','X'); +varid3 = netcdf.defVar(ncid,'lat','double',[dimid2 dimid1]); +netcdf.putAtt(ncid,varid3,'units','degrees north'); +varid4 = netcdf.defVar(ncid,'lon','double',[dimid2 dimid1]); +netcdf.putAtt(ncid,varid4,'units','degrees east'); +varid5 = netcdf.defVar(ncid,'NO3_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid5,'units','mol m-3'); +netcdf.putAtt(ncid,varid5,'long_name','DIN_CONC'); +varid6 = netcdf.defVar(ncid,'LDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid6,'units','mol m-3'); +netcdf.putAtt(ncid,varid6,'long_name','0.3*DON_CONC'); +varid7 = netcdf.defVar(ncid,'SLDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid7,'units','mol m-3'); +netcdf.putAtt(ncid,varid7,'long_name','0.35*DON_CONC'); +varid8 = netcdf.defVar(ncid,'SRDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid8,'units','mol m-3'); +netcdf.putAtt(ncid,varid8,'long_name','0.35*DON_CONC'); +varid9 = netcdf.defVar(ncid,'NDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid9,'units','mol m-3'); +netcdf.putAtt(ncid,varid9,'long_name','1.0*PN_CONC'); +varid10 = netcdf.defVar(ncid,'PO4_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid10,'units','mol m-3'); +netcdf.putAtt(ncid,varid10,'long_name','PO4_CONC'); +varid11 = netcdf.defVar(ncid,'LDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid11,'units','mol m-3'); +netcdf.putAtt(ncid,varid11,'long_name','0.3*DOP_CONC'); +varid12 = netcdf.defVar(ncid,'SLDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid12,'units','mol m-3'); +netcdf.putAtt(ncid,varid12,'long_name','0.35*DOP_CONC'); +varid13 = netcdf.defVar(ncid,'SRDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid13,'units','mol m-3'); +netcdf.putAtt(ncid,varid13,'long_name','0.35*DOP_CONC'); +varid14 = netcdf.defVar(ncid,'PDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid14,'units','mol m-3'); +netcdf.putAtt(ncid,varid14,'long_name','0.3*PP_CONC'); +varid15 = netcdf.defVar(ncid,'FED_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid15,'units','mol m-3'); +netcdf.putAtt(ncid,varid15,'long_name','FED_CONC'); +varid16 = netcdf.defVar(ncid,'FEDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid16,'units','mol m-3'); +netcdf.putAtt(ncid,varid16,'long_name','FEDET_CONC'); +varid17 = netcdf.defVar(ncid,'SI_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid17,'units','mol m-3'); +netcdf.putAtt(ncid,varid17,'long_name','SI_CONC'); +netcdf.close(ncid) + +ncid = netcdf.open(nc_file_name,'NC_WRITE'); +netcdf.putVar(ncid,varid0,0,1,time); +% nutrient input files appear seem to need dummy axes to be read in +% properly, but eventually do a grid by grid mapping that doesn't require +% these. +netcdf.putVar(ncid,varid1,1:nlat); +netcdf.putVar(ncid,varid2,1:nlon); +netcdf.putVar(ncid,varid3,permute(lon_mod,[2,1])); +netcdf.putVar(ncid,varid4,permute(lat_mod,[2,1])); +netcdf.putVar(ncid,varid5,permute(NO3_CONC,[3,2,1])); +netcdf.putVar(ncid,varid6,permute(LDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid7,permute(SLDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid8,permute(SRDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid9,permute(NDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid10,permute(PO4_CONC,[3,2,1])); +netcdf.putVar(ncid,varid11,permute(LDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid12,permute(SLDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid13,permute(SRDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid14,permute(PDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid15,permute(FED_CONC,[3,2,1])); +netcdf.putVar(ncid,varid16,permute(FEDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid17,permute(SI_CONC,[3,2,1])); +netcdf.close(ncid) diff --git a/tools/rivers/bgc/ARC/mapriv_combined_Arctic.m b/tools/rivers/bgc/ARC/mapriv_combined_Arctic.m new file mode 100644 index 000000000..dfcb9cafe --- /dev/null +++ b/tools/rivers/bgc/ARC/mapriv_combined_Arctic.m @@ -0,0 +1,1343 @@ +% Routine to map USGS nutrient data onto the MOM6 Northwest Atlantic (NWA) +% grid. Run on matlab97 or above. + +clear all; +addpath /home/cas/matlab +nc64startup + +% name of netcdf file to be created +nc_file_name = 'RiverNutrients_Combined_Q100_Arctic12k.nc'; + +% GLOBAL NEWS based map for filling in gaps +NEWS_file = 'RiverNutrients_GlobalNEWS2_plusFe_Q100_Arctic12k.nc'; + +% load in monthly world ocean T, S climatology for saturated oxygen calculation +temp = ncread('Data/woa_sst_climo.nc','t_an'); +woa_temp = permute(temp,[3 2 1]); + +% Parameters for the assignment algorithm. +Q_min = 100; % minimum flow in m3 sec +plot_width = 30; % width of window (in degrees) for inspecting locations + % of rivers and outflow points that have been assigned to + % them. +min_dist = 1.5; % minimum distance (degrees) of the closest outflow point + % for the river to be considered in the domain (useful + % for preventing the algorithm from trying to map rivers + % flowing to different ocean basins. +max_dist = 2.0; % maximum distance (degrees) away that the algorithm + % looks for points for rivers that are in the domain +nutrient_option = 2; % option for deriving dissolved organic nutrients in RC4US +inspect_map = 'y'; % flag enabling you to pause and inspect each river + % mapping as it is being done. + +min_lon_ref = -180;% set to either 0 if model grid contains no negative + % values; set to -180 if model is on a -180-180 grid. + + +% set the bio-availability of phosphorus and the fractionation of dissolved +% organic; PP is set to 30% based on Froelich; Partitioning of detritus +% between +frac_PP = 0.3; +frac_ldon = 0.3; +frac_sldon = 0.35; +frac_srdon = 0.35; +frac_ldop = 0.3; +frac_sldop = 0.35; +frac_srdop = 0.35; +% 40 nM dissolved iron concentration from De Baar and De Jong + 30nM +% Colloidal and nanoparticle flux as reported in Canfield and Raiswell +const_fed = 70.0e-6; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% USGS data compiled by Fabian Gomez % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +filename_chem = 'Data/RC4USCoast/mclim_19902022_chem.nc'; +alk_monthly_RC4US = ncread(filename_chem,'alk'); alk_monthly_RC4US = permute(alk_monthly_RC4US,[2 1]); +dic_monthly_RC4US = ncread(filename_chem,'dic'); dic_monthly_RC4US = permute(dic_monthly_RC4US,[2 1]); +no3_monthly_RC4US = ncread(filename_chem,'no3'); no3_monthly_RC4US = permute(no3_monthly_RC4US,[2 1]); +nh4_monthly_RC4US = ncread(filename_chem,'nh4'); nh4_monthly_RC4US = permute(nh4_monthly_RC4US,[2 1]); +din_monthly_RC4US = no3_monthly_RC4US + nh4_monthly_RC4US; +dip_monthly_RC4US = ncread(filename_chem,'po4'); dip_monthly_RC4US = permute(dip_monthly_RC4US,[2 1]); +si_monthly_RC4US = ncread(filename_chem,'sio2'); si_monthly_RC4US = permute(si_monthly_RC4US,[2 1]); +% The RC4US database seems to be in mmoles O m-3 rather than mmoles O2 m-3, +% divide by 2.0 for consistency with other O2 data sources +o2_monthly_RC4US = ncread(filename_chem,'do'); o2_monthly_RC4US = permute(o2_monthly_RC4US,[2 1])/2.0; +don_monthly_RC4US = ncread(filename_chem,'don'); don_monthly_RC4US = permute(don_monthly_RC4US,[2 1]); +temp_monthly_RC4US = ncread(filename_chem,'temp'); temp_monthly_RC4US = permute(temp_monthly_RC4US,[2 1]); +if nutrient_option == 1 + % This option will eventually set all river values for pn, dop and pp using GlobalNEWS + pn_monthly_RC4US = no3_monthly_RC4US*NaN; + dop_monthly_RC4US = no3_monthly_RC4US*NaN; + pp_monthly_RC4US = no3_monthly_RC4US*NaN; +elseif nutrient_option == 2 + % This option will use differences between total filtered and unfiltered + % and other properties to derive pn, dop and pp. This unfortunately + % generates negative values in some cases. + tnf_monthly_RC4US = ncread(filename_chem,'tnf'); tnf_monthly_RC4US = permute(tnf_monthly_RC4US,[2 1]); + don2_monthly_RC4US = tnf_monthly_RC4US - din_monthly_RC4US; + tnu_monthly_RC4US = ncread(filename_chem,'tnu'); tnu_monthly_RC4US = permute(tnu_monthly_RC4US,[2 1]); + pn_monthly_RC4US = tnu_monthly_RC4US - tnf_monthly_RC4US; + pn_monthly_RC4US(pn_monthly_RC4US < 0) = NaN; + tpf_monthly_RC4US = ncread(filename_chem,'tpf'); tpf_monthly_RC4US = permute(tpf_monthly_RC4US,[2 1]); + dop_monthly_RC4US = tpf_monthly_RC4US - dip_monthly_RC4US; + dop_monthly_RC4US(dop_monthly_RC4US < 0) = NaN; + tpu_monthly_RC4US = ncread(filename_chem,'tpu'); tpu_monthly_RC4US = permute(tpu_monthly_RC4US,[2 1]); + pp_monthly_RC4US = (tpu_monthly_RC4US - tpf_monthly_RC4US)*frac_PP; + pp_monthly_RC4US(pp_monthly_RC4US < 0) = NaN; +end +dfe_monthly_RC4US = no3_monthly_RC4US*NaN; +pfe_monthly_RC4US = no3_monthly_RC4US*NaN; + +filename_discharge = 'Data/RC4USCoast/mclim_19902022_disc.nc'; +Q_monthly_RC4US = ncread(filename_discharge,'disc'); Q_monthly_RC4US = permute(Q_monthly_RC4US,[2 1]); % m-3 sec-1 +station_names_RC4US = h5read(filename_discharge,'/river_name'); +lon_stations_RC4US = ncread(filename_discharge,'mouth_lon'); +lat_stations_RC4US = ncread(filename_discharge,'mouth_lat'); + +Q_ann_RC4US = mean(Q_monthly_RC4US,1,'native','omitnan')'; +dic_ann_RC4US = mean(dic_monthly_RC4US,1,'native','omitnan')'; +alk_ann_RC4US = mean(alk_monthly_RC4US,1,'native','omitnan')'; +no3_ann_RC4US = mean(no3_monthly_RC4US,1,'native','omitnan')'; +nh4_ann_RC4US = mean(nh4_monthly_RC4US,1,'native','omitnan')'; +o2_ann_RC4US = mean(o2_monthly_RC4US,1,'native','omitnan')'; +dip_ann_RC4US = mean(dip_monthly_RC4US,1,'native','omitnan')'; +si_ann_RC4US = mean(si_monthly_RC4US,1,'native','omitnan')'; +din_ann_RC4US = no3_ann_RC4US + nh4_ann_RC4US; +don_ann_RC4US = mean(don_monthly_RC4US,1,'native','omitnan')'; +if nutrient_option == 1 + don_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + pn_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + dop_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + pp_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; +elseif nutrient_option == 2 + don2_ann_RC4US = mean(don2_monthly_RC4US,1,'native','omitnan')'; + pn_ann_RC4US = mean(pn_monthly_RC4US,1,'native','omitnan')'; + dop_ann_RC4US = mean(dop_monthly_RC4US,1,'native','omitnan')'; + pp_ann_RC4US = mean(pp_monthly_RC4US,1,'native','omitnan')'; +end +dfe_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; +pfe_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + +for n = 1:size(lon_stations_RC4US,1) + + % Make any adjustments to the river locations here, e.g: + % + % Move the Susquehanna a bit south so that it catches the Chesapeake + % and not the Delaware. + %if strcmp('Susquehanna',station_names_RC4US{n}) + % lat_stations_RC4US(n) = 38.5; + % lon_stations_RC4US(n) = -77.5; + % %pause + %end + +end + +Q_ann_RC4US = mean(Q_monthly_RC4US,1,'native','omitnan')'; +dic_ann_RC4US = mean(dic_monthly_RC4US,1,'native','omitnan')'; +alk_ann_RC4US = mean(alk_monthly_RC4US,1,'native','omitnan')'; +no3_ann_RC4US = mean(no3_monthly_RC4US,1,'native','omitnan')'; +nh4_ann_RC4US = mean(nh4_monthly_RC4US,1,'native','omitnan')'; +o2_ann_RC4US = mean(o2_monthly_RC4US,1,'native','omitnan')'; +dip_ann_RC4US = mean(dip_monthly_RC4US,1,'native','omitnan')'; +si_ann_RC4US = mean(si_monthly_RC4US,1,'native','omitnan')'; +din_ann_RC4US = no3_ann_RC4US + nh4_ann_RC4US; +don_ann_RC4US = mean(don_monthly_RC4US,1,'native','omitnan')'; +if nutrient_option == 1 + don_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + pn_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + dop_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + pp_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; +elseif nutrient_option == 2 + don2_ann_RC4US = mean(don2_monthly_RC4US,1,'native','omitnan')'; + pn_ann_RC4US = mean(pn_monthly_RC4US,1,'native','omitnan')'; + dop_ann_RC4US = mean(dop_monthly_RC4US,1,'native','omitnan')'; + pp_ann_RC4US = mean(pp_monthly_RC4US,1,'native','omitnan')'; +end +dfe_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; +pfe_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + +% Created by Process_GLORICH_NEP.m, includes following variables: +%save Arctic_GLORICH_data lon_stations_glorich lat_stations_glorich station_names_glorich ... +% Q_ann_glorich dic_ann_glorich alk_ann_glorich no3_ann_glorich nh4_ann_glorich ... +% din_ann_glorich pn_ann_glorich don_ann_glorich dip_ann_glorich dop_ann_glorich ... +% pp_ann_glorich si_ann_glorich o2_ann_glorich dfe_ann_glorich pfe_ann_glorich ... +% dic_monthly_glorich alk_monthly_glorich no3_monthly_glorich nh4_monthly_glorich ... +% din_monthly_glorich don_monthly_glorich pn_monthly_glorich dip_monthly_glorich ... +% dop_monthly_glorich pp_monthly_glorich dfe_monthly_glorich pfe_monthly_glorich ... +% si_monthly_glorich o2_monthly_glorich +load Data/GLORICH/Arctic_GLORICH_data.mat; + +% Created by Process_ARCTICGRO.m, includes following variables: +%save arcticgro_data lon_stations_arcticgro lat_stations_arcticgro station_names_arcticgro ... +% Q_ann_arcticgro dic_ann_arcticgro alk_ann_arcticgro no3_ann_arcticgro nh4_ann_arcticgro ... +% din_ann_arcticgro pn_ann_arcticgro don_ann_arcticgro dip_ann_arcticgro dop_ann_arcticgro ... +% pp_ann_arcticgro si_ann_arcticgro o2_ann_arcticgro dfe_ann_arcticgro pfe_ann_arcticgro ... +% dic_monthly_arcticgro alk_monthly_arcticgro no3_monthly_arcticgro nh4_monthly_arcticgro ... +% din_monthly_arcticgro don_monthly_arcticgro pn_monthly_arcticgro dip_monthly_arcticgro ... +% dop_monthly_arcticgro pp_monthly_arcticgro dfe_monthly_arcticgro pfe_monthly_arcticgro ... +% si_monthly_arcticgro o2_monthly_arcticgro +load Data/ArcticGro/arcticgro_data.mat; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Combine all annual and monthly station data % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +station_names_all = [station_names_RC4US; string(station_names_glorich)'; string(station_names_arcticgro)']; +lon_stations_all = [lon_stations_RC4US; lon_stations_glorich'; lon_stations_arcticgro']; +if min_lon_ref == 0 + lon_stations_all(lon_stations_all < 0) = lon_stations_all(lon_stations_all < 0) + 360; +elseif min_lon_ref == -180 + lon_stations_all(lon_stations_all > 180) = lon_stations_all(lon_stations_all > 180) - 360; +end +lat_stations_all = [lat_stations_RC4US; lat_stations_glorich'; lat_stations_arcticgro']; +Q_ann_all = [Q_ann_RC4US; Q_ann_glorich'; Q_ann_arcticgro']; + +dic_ann_all = [dic_ann_RC4US; dic_ann_glorich'; dic_ann_arcticgro']; +alk_ann_all = [alk_ann_RC4US; alk_ann_glorich'; alk_ann_arcticgro']; +no3_ann_all = [no3_ann_RC4US; no3_ann_glorich'; no3_ann_arcticgro']; +nh4_ann_all = [nh4_ann_RC4US; nh4_ann_glorich'; nh4_ann_arcticgro']; +din_ann_all = [din_ann_RC4US; din_ann_glorich'; din_ann_arcticgro']; +don_ann_all = [don_ann_RC4US; don_ann_glorich'; don_ann_arcticgro']; +pn_ann_all = [pn_ann_RC4US; pn_ann_glorich'; pn_ann_arcticgro']; +dip_ann_all = [dip_ann_RC4US; dip_ann_glorich'; dip_ann_arcticgro']; +dop_ann_all = [dop_ann_RC4US; dop_ann_glorich'; dop_ann_arcticgro']; +pp_ann_all = [pp_ann_RC4US; pp_ann_glorich'; pp_ann_arcticgro']; +dfe_ann_all = [dfe_ann_RC4US; dfe_ann_glorich'; dfe_ann_arcticgro']; +pfe_ann_all = [pfe_ann_RC4US; pfe_ann_glorich'; pfe_ann_arcticgro']; +si_ann_all = [si_ann_RC4US; si_ann_glorich'; si_ann_arcticgro']; +o2_ann_all = [o2_ann_RC4US; o2_ann_glorich'; o2_ann_arcticgro']; + +dic_monthly_all = [dic_monthly_RC4US dic_monthly_glorich dic_monthly_arcticgro]; +alk_monthly_all = [alk_monthly_RC4US alk_monthly_glorich alk_monthly_arcticgro]; +no3_monthly_all = [no3_monthly_RC4US no3_monthly_glorich no3_monthly_arcticgro]; +nh4_monthly_all = [nh4_monthly_RC4US nh4_monthly_glorich nh4_monthly_arcticgro]; +din_monthly_all = [din_monthly_RC4US din_monthly_glorich din_monthly_arcticgro]; +don_monthly_all = [don_monthly_RC4US don_monthly_glorich don_monthly_arcticgro]; +pn_monthly_all = [pn_monthly_RC4US pn_monthly_glorich pn_monthly_arcticgro]; +dip_monthly_all = [dip_monthly_RC4US dip_monthly_glorich dip_monthly_arcticgro]; +dop_monthly_all = [dop_monthly_RC4US dop_monthly_glorich dop_monthly_arcticgro]; +pp_monthly_all = [pp_monthly_RC4US pp_monthly_glorich pp_monthly_arcticgro]; +dfe_monthly_all = [dfe_monthly_RC4US dfe_monthly_glorich dfe_monthly_arcticgro]; +pfe_monthly_all = [pfe_monthly_RC4US pfe_monthly_glorich pfe_monthly_arcticgro]; +si_monthly_all = [si_monthly_RC4US si_monthly_glorich si_monthly_arcticgro]; +o2_monthly_all = [o2_monthly_RC4US o2_monthly_glorich o2_monthly_arcticgro]; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Load in monthly climatology of river forcing from the regional grid. % +% File contains: % +% runoff: monthly average runoff in kg m-2 sec-1 % +% area: area of grid cell in m-2 % +% lon: longitude (0-360 degrees) % +% lat: latitude % % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +load glofas_hill_runoff_monthlyclim_arctic12k_05112023.mat; + +lon_mod = lon; +lat_mod = lat; +area_mod = area; +% convert runoff from kg m-2 sec-1 to m3 sec-1 +Q_mod_monthly = zeros(size(runoff)); +for m = 1:12 + Q_mod_monthly(m,:,:) = squeeze(runoff(m,:,:)).*area_mod./1000; +end +Q_mod_ann = squeeze(mean(Q_mod_monthly,1)); +clear lon lat runoff area; + +%grid_file = '/archive/cas/Regional_MOM6/NWA12/nwa12_ocean_static.nc'; +%temp = ncread(grid_file,'deptho'); +%depth = permute(temp,[2,1]); clear temp; +%depth(isnan(depth)) = -1; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Filter for rivers in the region, set thresholds for minimum river size, % +% set parameters for plotting routines. % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% use grid to filter rivers outside domain +lat_mod_max = max(lat_mod(:)); +lat_mod_min = min(lat_mod(:)); +lon_mod_max = max(lon_mod(:)); +lon_mod_min = min(lon_mod(:)); + +in_region = find( (lon_stations_all <= lon_mod_max) & (lon_stations_all >= lon_mod_min) & ... + (lat_stations_all <= lat_mod_max) & (lat_stations_all >= lat_mod_min) & ... + (isfinite(Q_ann_all)) & (Q_ann_all > Q_min) ); +num_rivers = size(in_region,1); + +station_names_reg = station_names_all(in_region); +lon_stations_reg = lon_stations_all(in_region); +lat_stations_reg = lat_stations_all(in_region); +Q_ann_reg = Q_ann_all(in_region); +dic_ann_reg = dic_ann_all(in_region); +alk_ann_reg = alk_ann_all(in_region); +no3_ann_reg = no3_ann_all(in_region); +nh4_ann_reg = nh4_ann_all(in_region); +din_ann_reg = din_ann_all(in_region); +don_ann_reg = don_ann_all(in_region); +pn_ann_reg = pn_ann_all(in_region); +dip_ann_reg = dip_ann_all(in_region); +dop_ann_reg = dop_ann_all(in_region); +pp_ann_reg = pp_ann_all(in_region); +dfe_ann_reg = dfe_ann_all(in_region); +pfe_ann_reg = pfe_ann_all(in_region); +si_ann_reg = si_ann_all(in_region); +o2_ann_reg = o2_ann_all(in_region); + +for m = 1:12 + dic_monthly_reg(m,:) = dic_monthly_all(m,in_region); + alk_monthly_reg(m,:) = alk_monthly_all(m,in_region); + no3_monthly_reg(m,:) = no3_monthly_all(m,in_region); + nh4_monthly_reg(m,:) = nh4_monthly_all(m,in_region); + din_monthly_reg(m,:) = din_monthly_all(m,in_region); + don_monthly_reg(m,:) = don_monthly_all(m,in_region); + pn_monthly_reg(m,:) = pn_monthly_all(m,in_region); + dip_monthly_reg(m,:) = dip_monthly_all(m,in_region); + dop_monthly_reg(m,:) = dop_monthly_all(m,in_region); + pp_monthly_reg(m,:) = pp_monthly_all(m,in_region); + dfe_monthly_reg(m,:) = dfe_monthly_all(m,in_region); + pfe_monthly_reg(m,:) = pfe_monthly_all(m,in_region); + si_monthly_reg(m,:) = si_monthly_all(m,in_region); + o2_monthly_reg(m,:) = o2_monthly_all(m,in_region); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Assigning outflow points to rivers. % +% 1. Assignment starts with the rivers with the smallest flow and works % +% to the largest, w/larger river characteristics taking precedence to % +% ensure the most significant rivers are well represented. % +% 2. The algorithm keeps choosing the closest points to each river mouth % +% until the assigned flow is as close as possible to that observed % +% 3. Once the outflow points are assigned using the mean flow values, % +% monthly concentrations are assigned to those points. % +% 4. A simple "nearest neighbor" algorithm is used to fill in the gaps % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Sort rivers by discharge +[Q_ann_sort,sort_ind] = sort(Q_ann_reg,'ascend'); + +station_names_sort = station_names_reg(sort_ind); +lon_stations_sort = lon_stations_reg(sort_ind); +lat_stations_sort = lat_stations_reg(sort_ind); +Q_ann_sort = Q_ann_reg(sort_ind); +dic_ann_sort = dic_ann_reg(sort_ind); +alk_ann_sort = alk_ann_reg(sort_ind); +no3_ann_sort = no3_ann_reg(sort_ind); +nh4_ann_sort = nh4_ann_reg(sort_ind); +din_ann_sort = din_ann_reg(sort_ind); +don_ann_sort = don_ann_reg(sort_ind); +pn_ann_sort = pn_ann_reg(sort_ind); +dip_ann_sort = dip_ann_reg(sort_ind); +dop_ann_sort = dop_ann_reg(sort_ind); +pp_ann_sort = pp_ann_reg(sort_ind); +dfe_ann_sort = dfe_ann_reg(sort_ind); +pfe_ann_sort = pfe_ann_reg(sort_ind); +si_ann_sort = si_ann_reg(sort_ind); +o2_ann_sort = o2_ann_reg(sort_ind); + +for m = 1:12 + dic_monthly_sort(m,:) = dic_monthly_reg(m,sort_ind); + alk_monthly_sort(m,:) = alk_monthly_reg(m,sort_ind); + no3_monthly_sort(m,:) = no3_monthly_reg(m,sort_ind); + nh4_monthly_sort(m,:) = nh4_monthly_reg(m,sort_ind); + din_monthly_sort(m,:) = din_monthly_reg(m,sort_ind); + don_monthly_sort(m,:) = don_monthly_reg(m,sort_ind); + pn_monthly_sort(m,:) = pn_monthly_reg(m,sort_ind); + dip_monthly_sort(m,:) = dip_monthly_reg(m,sort_ind); + dop_monthly_sort(m,:) = dop_monthly_reg(m,sort_ind); + pp_monthly_sort(m,:) = pp_monthly_reg(m,sort_ind); + dfe_monthly_sort(m,:) = dfe_monthly_reg(m,sort_ind); + pfe_monthly_sort(m,:) = pfe_monthly_reg(m,sort_ind); + si_monthly_sort(m,:) = si_monthly_reg(m,sort_ind); + o2_monthly_sort(m,:) = o2_monthly_reg(m,sort_ind); +end + +% Create vectors of values at the runoff points from the model grid. These +% are used to accelerate the mapping relative to wrangling the full grid +% with all the zeros included. "ind_ro" are the grid indexes with runoff +ind_ro = find(Q_mod_ann > 0); +Q_mod_vec = Q_mod_ann(ind_ro); +lon_mod_runoff_vec = lon_mod(ind_ro); +lat_mod_runoff_vec = lat_mod(ind_ro); +Q_mod_monthly_vecs = zeros(12,size(lon_mod_runoff_vec,1)); +for m = 1:12 + temp = squeeze(Q_mod_monthly(m,:,:)); + Q_mod_monthly_vecs(m,:) = temp(ind_ro); +end + +% Create a grid of saturated oxygen values using the world ocean atlas data +temp_woa_monthly_vecs = zeros(12,size(lon_mod_runoff_vec,1)); +o2sat_woa_monthly_vecs = zeros(12,size(lon_mod_runoff_vec,1)); + +% Constants for o2 saturation calculation (taken from COBALT) +a_0 = 2.00907; +a_1 = 3.22014; +a_2 = 4.05010; +a_3 = 4.94457; +a_4 = -2.56847e-1; +a_5 = 3.88767; +sal = 0; +b_0 = -6.24523e-3; +b_1 = -7.37614e-3; +b_2 = -1.03410e-2; +b_3 = -8.17083e-3; +c_0 = -4.88682e-7; + +for m = 1:12 + temp = squeeze(woa_temp(m,:,:)); + % limit for validity of o2sat calculation + temp(temp > 40) = 40; temp(temp < 0) = 0; + temp_woa_monthly_vecs(m,:) = temp(ind_ro); + % calculate the oxygen saturation at a given temperature and salinity = 0 + % code taken from COBALT with limits applied above + tt = 298.15 - temp_woa_monthly_vecs(m,:); + tkb = 273.15 + temp_woa_monthly_vecs(m,:); + ts = log(tt / tkb); + ts2 = ts * ts; + ts3 = ts2 * ts; + ts4 = ts3 * ts; + ts5 = ts4 * ts; + + o2sat_woa_monthly_vecs(m,:) = (1000.0/22391.6) * 1000 * ... %convert from ml/l to mmol m-3 + exp(a_0 + a_1*ts + a_2*ts2 + a_3*ts3 + a_4*ts4 + a_5*ts5 + ... + (b_0 + b_1*ts + b_2*ts2 + b_3*ts3 + c_0*sal)*sal); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Load in fields generated from global NEWS. Where necessary, the ratio % +% of constituents relative to DIN will be used to fill forcing gaps % +% The m-file used to generate the NEWS forcing file is included in this % +% directory and uses an analogous mapping algorithm to this one % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +din_ann_NEWS = nc_varget(NEWS_file,'NO3_CONC'); +aa = find(din_ann_NEWS > 0); +temp1 = nc_varget(NEWS_file,'LDON_CONC'); +temp2 = nc_varget(NEWS_file,'SLDON_CONC'); +temp3 = nc_varget(NEWS_file,'SRDON_CONC'); +don_ann_NEWS = temp1 + temp2 + temp3; +don_ratio_NEWS_vec = don_ann_NEWS(aa)./din_ann_NEWS(aa); +clear temp1 temp2 temp3; +pn_ann_NEWS = nc_varget(NEWS_file,'NDET_CONC'); +pn_ratio_NEWS_vec = pn_ann_NEWS(aa)./din_ann_NEWS(aa); + +dip_ann_NEWS = nc_varget(NEWS_file,'PO4_CONC'); +dip_ratio_NEWS_vec = dip_ann_NEWS(aa)./din_ann_NEWS(aa); +temp1 = nc_varget(NEWS_file,'LDOP_CONC'); +temp2 = nc_varget(NEWS_file,'SLDOP_CONC'); +temp3 = nc_varget(NEWS_file,'SRDOP_CONC'); +dop_ann_NEWS = temp1 + temp2 + temp3; +dop_ratio_NEWS_vec = dop_ann_NEWS(aa)./din_ann_NEWS(aa); +clear temp1 temp2 temp3; +pp_ann_NEWS = nc_varget(NEWS_file,'PDET_CONC'); +pp_ratio_NEWS_vec = pp_ann_NEWS(aa)./din_ann_NEWS(aa); +si_ann_NEWS = nc_varget(NEWS_file,'SI_CONC'); +si_ratio_NEWS_vec = si_ann_NEWS(aa)./din_ann_NEWS(aa); + +% Vectors to hold monthly values mapped onto model runoff points +dic_mod_monthly_vecs = zeros(12,size(lon_mod_runoff_vec,1)); +alk_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +no3_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +nh4_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +din_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +don_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +pn_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +dip_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +dop_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +pp_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +dfe_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +pfe_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +si_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +o2_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Loop identifies points assigned to each river % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +for k=1:num_rivers + dist = pdist2([lon_stations_sort(k) lat_stations_sort(k)], ... + [lon_mod_runoff_vec lat_mod_runoff_vec]); + [dist_sort, dist_sort_ind] = sort(dist,'ascend'); + + if dist_sort(1) < min_dist + Q_sum1 = 0; + Q_sum2 = 0; + n = 0; + while (Q_sum2 < Q_ann_sort(k) && (dist_sort(n+1) < max_dist)) + Q_sum1 = Q_sum2; + n = n+1; + Q_sum2 = Q_sum1 + Q_mod_vec(dist_sort_ind(n)); + end + %if abs(Q_sum1 - Q_ann_sort(k)) < abs(Q_sum2 - Q_ann_sort(k)) + % nrp = n-1; % number of runoff points + % [Q_sum1 Q_ann_sort(k)] % a quick check for comparable flow + %else + nrp = n; + [Q_sum2 Q_ann_sort(k)] + %end + + % enter monthly concentration values into an array of monthly values + % if no monthly value is available, use annuals. + for m = 1:12 + + if isnan(dic_monthly_sort(m,k)) + if isfinite(dic_ann_sort(k)) + dic_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dic_ann_sort(k); + else + dic_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = 0; + end + else + dic_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dic_monthly_sort(m,k); + end + + if isnan(alk_monthly_sort(m,k)) + if isfinite(dic_ann_sort(k)) + alk_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = alk_ann_sort(k); + else + alk_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = 0; + end + else + alk_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = alk_monthly_sort(m,k); + end + + % mapping assumes that DIN is defined for nutrient calculations since + % ratios relative to DIN are used to fill in other components. If + % DIN is not defined, values are left at 0 and eventually filled with + % a nearest neighbor filling (next section) + + if (isfinite(din_monthly_sort(m,k)) || isfinite(din_ann_sort(k)) ) + if isfinite(din_monthly_sort(m,k)) + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = din_monthly_sort(m,k); + elseif isfinite(din_ann_sort(k)) + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = din_ann_sort(k); + end + + if isfinite(no3_monthly_sort(m,k)) + no3_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = no3_monthly_sort(m,k); + elseif isfinite(no3_ann_sort(k)) + no3_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = no3_ann_sort(k); + end + + if isfinite(nh4_monthly_sort(m,k)) + nh4_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = nh4_monthly_sort(m,k); + elseif isfinite(nh4_ann_sort(k)) + nh4_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = nh4_ann_sort(k); + end + + if (isnan(don_monthly_sort(m,k)) && isnan(don_ann_sort(k))) + don_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + don_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(don_monthly_sort(m,k)) && ~isnan(don_ann_sort(k))) + don_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = don_ann_sort(k); + else + don_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = don_monthly_sort(m,k); + end + + if (isnan(pn_monthly_sort(m,k)) && isnan(pn_ann_sort(k))) + pn_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + pn_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(pn_monthly_sort(m,k)) && ~isnan(pn_ann_sort(k))) + pn_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = pn_ann_sort(k); + else + pn_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = pn_monthly_sort(m,k); + end + + if (isnan(dip_monthly_sort(m,k)) && isnan(dip_ann_sort(k))) + dip_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + dip_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(dip_monthly_sort(m,k)) && ~isnan(dip_ann_sort(k))) + dip_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dip_ann_sort(k); + else + dip_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dip_monthly_sort(m,k); + end + + if (isnan(dop_monthly_sort(m,k)) && isnan(dop_ann_sort(k))) + dop_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + dop_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(dop_monthly_sort(m,k)) && ~isnan(dop_ann_sort(k))) + dop_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dop_ann_sort(k); + else + dop_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dop_monthly_sort(m,k); + end + + if (isnan(pp_monthly_sort(m,k)) && isnan(pp_ann_sort(k))) + pp_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + pp_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(pp_monthly_sort(m,k)) && ~isnan(pp_ann_sort(k))) + pp_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = pp_ann_sort(k); + else + pp_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = pp_monthly_sort(m,k); + end + + if (isnan(si_monthly_sort(m,k)) && isnan(si_ann_sort(k))) + si_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + si_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(si_monthly_sort(m,k)) && ~isnan(si_ann_sort(k))) + si_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = si_ann_sort(k); + else + si_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = si_monthly_sort(m,k); + end + + if isnan(o2_monthly_sort(m,k)) + o2_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + o2sat_woa_monthly_vecs(m,dist_sort_ind(1:nrp)); + else + o2_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = o2_monthly_sort(m,k); + end + + end + end + + % plot to check location if inspect_map == 'y'. The plot puts + % open circles at each runoff location in the model grid and fills + % those that are assigned to each river. Note that some of the smaller + % rivers may be replaced with larger ones as the fitting process + % continues. + + if inspect_map == 'y' + figure(1) + clf + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(Q_mod_vec),3,log10(Q_mod_vec)); + hold on + scatter3(lon_mod_runoff_vec(dist_sort_ind(1:nrp)),lat_mod_runoff_vec(dist_sort_ind(1:nrp)), ... + log10(Q_mod_vec(dist_sort_ind(1:nrp))),25, ... + log10(Q_mod_vec(dist_sort_ind(1:nrp))),'filled'); + view(2); + plot3(lon_stations_sort(k),lat_stations_sort(k),1e5,'k.','MarkerSize',20); + %contour(lon_mod,lat_mod,depth,[0 0],'k-'); + axis([lon_stations_sort(k)-plot_width/2 lon_stations_sort(k)+plot_width/2 ... + lat_stations_sort(k)-plot_width/2 lat_stations_sort(k)+plot_width/2]); + caxis([-4 3]); + titl = ['river number: ',num2str(k),' name: ',station_names_sort{k}]; + title(titl); + colorbar; + + % check the values of each mapping + N_check = mean(din_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all') + ... + mean(don_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all') + ... + mean(pn_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + P_check = mean(dip_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all') + ... + mean(dop_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all') + ... + mean(pp_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + SI_check = mean(si_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + DIN_check = mean(din_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + DON_check = mean(don_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + PN_check = mean(pn_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + DIP_check = mean(dip_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + DOP_check = mean(dop_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + PP_check = mean(pp_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + SI_check = mean(si_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + + station_names_sort(k) + ind = dist_sort_ind(1:nrp); + 'total flow in m3 sec' + [Q_ann_sort(k) sum(Q_mod_vec(dist_sort_ind(1:nrp)))] + 'N, P conc (mmoles m-3), DI, DO, P' + [DIN_check DON_check PN_check] + [DIP_check DOP_check PP_check] + 'Total N, Total P, Total N: Total P' + [N_check P_check N_check/P_check] + 'DO:DI and P:DI ratios'; + [DON_check/DIN_check PN_check/DIN_check]; + [DOP_check/DIP_check PP_check/DIP_check]; + 'silica concentration (mmoles m-3)'; + SI_check; + + pause + end + + % If river is outside the domain, skip all of the calculations above and + % just plot for inspection/evaluation + else + % This is for rivers that were outside of the domain + if inspect_map == 'y' + figure(1) + clf + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(Q_mod_vec),3,log10(Q_mod_vec)); + hold on + view(2); + plot3(lon_stations_sort(k),lat_stations_sort(k),1e5,'k.','MarkerSize',20); + axis([lon_stations_sort(k)-10 lon_stations_sort(k)+10 ... + lat_stations_sort(k)-10 lat_stations_sort(k)+10]); + caxis([-4 3]); + titl = ['OUTSIDE: river number: ',num2str(k),' name: ',station_names_sort{k}]; + title(titl); + colorbar; + + pause + end + + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% nearest neighbor search to fill in any runoff points that were not % +% assigned after the runoff mapping step % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +lon_mod_runoff_vec = double(lon_mod_runoff_vec); +lat_mod_runoff_vec = double(lat_mod_runoff_vec); + +for m = 1:12 + aa = find(dic_mod_monthly_vecs(m,:) == 0); + bb = find(dic_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + dic_mod_monthly_vecs(m,bb)','nearest','nearest'); + dic_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(alk_mod_monthly_vecs(m,:) == 0); + bb = find(alk_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + alk_mod_monthly_vecs(m,bb)','nearest','nearest'); + alk_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(no3_mod_monthly_vecs(m,:) == 0); + bb = find(no3_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + no3_mod_monthly_vecs(m,bb)','nearest','nearest'); + no3_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(nh4_mod_monthly_vecs(m,:) == 0); + bb = find(nh4_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + nh4_mod_monthly_vecs(m,bb)','nearest','nearest'); + nh4_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(din_mod_monthly_vecs(m,:) == 0); + bb = find(din_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + din_mod_monthly_vecs(m,bb)','nearest','nearest'); + din_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(don_mod_monthly_vecs(m,:) == 0); + bb = find(don_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + don_mod_monthly_vecs(m,bb)','nearest','nearest'); + don_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(pn_mod_monthly_vecs(m,:) == 0); + bb = find(pn_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + pn_mod_monthly_vecs(m,bb)','nearest','nearest'); + pn_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(dip_mod_monthly_vecs(m,:) == 0); + bb = find(dip_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + dip_mod_monthly_vecs(m,bb)','nearest','nearest'); + dip_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(dop_mod_monthly_vecs(m,:) == 0); + bb = find(dop_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + dop_mod_monthly_vecs(m,bb)','nearest','nearest'); + dop_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(pp_mod_monthly_vecs(m,:) == 0); + bb = find(pp_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + pp_mod_monthly_vecs(m,bb)','nearest','nearest'); + pp_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(si_mod_monthly_vecs(m,:) == 0); + bb = find(si_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + si_mod_monthly_vecs(m,bb)','nearest','nearest'); + si_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); +end + +% For o2sat, fill in any 0 values with saturated o2 at the world ocean +% atlas climatology +for m = 1:12 + aa = find(o2_mod_monthly_vecs(m,:) == 0); + o2_mod_monthly_vecs(m,aa) = o2sat_woa_monthly_vecs(m,aa); +end + +totn_mod_monthly_vecs = din_mod_monthly_vecs + don_mod_monthly_vecs + pn_mod_monthly_vecs; +totp_mod_monthly_vecs = dip_mod_monthly_vecs + dop_mod_monthly_vecs + pp_mod_monthly_vecs; + +dicflux_mod_monthly_vecs = dic_mod_monthly_vecs.*Q_mod_monthly_vecs; +alkflux_mod_monthly_vecs = alk_mod_monthly_vecs.*Q_mod_monthly_vecs; +dinflux_mod_monthly_vecs = din_mod_monthly_vecs.*Q_mod_monthly_vecs; +no3flux_mod_monthly_vecs = no3_mod_monthly_vecs.*Q_mod_monthly_vecs; +nh4flux_mod_monthly_vecs = nh4_mod_monthly_vecs.*Q_mod_monthly_vecs; +dipflux_mod_monthly_vecs = dip_mod_monthly_vecs.*Q_mod_monthly_vecs; +donflux_mod_monthly_vecs = don_mod_monthly_vecs.*Q_mod_monthly_vecs; +dopflux_mod_monthly_vecs = dop_mod_monthly_vecs.*Q_mod_monthly_vecs; +pnflux_mod_monthly_vecs = pn_mod_monthly_vecs.*Q_mod_monthly_vecs; +ppflux_mod_monthly_vecs = pp_mod_monthly_vecs.*Q_mod_monthly_vecs; +siflux_mod_monthly_vecs = si_mod_monthly_vecs.*Q_mod_monthly_vecs; +totnflux_mod_monthly_vecs = totn_mod_monthly_vecs.*Q_mod_monthly_vecs; +totpflux_mod_monthly_vecs = totp_mod_monthly_vecs.*Q_mod_monthly_vecs; +o2flux_mod_monthly_vecs = o2_mod_monthly_vecs.*Q_mod_monthly_vecs; + +dicflux_mod_ann_vec = mean(dicflux_mod_monthly_vecs,1); +alkflux_mod_ann_vec = mean(alkflux_mod_monthly_vecs,1); +dinflux_mod_ann_vec = mean(dinflux_mod_monthly_vecs,1); +no3flux_mod_ann_vec = mean(no3flux_mod_monthly_vecs,1); +nh4flux_mod_ann_vec = mean(nh4flux_mod_monthly_vecs,1); +dipflux_mod_ann_vec = mean(dipflux_mod_monthly_vecs,1); +donflux_mod_ann_vec = mean(donflux_mod_monthly_vecs,1); +dopflux_mod_ann_vec = mean(dopflux_mod_monthly_vecs,1); +pnflux_mod_ann_vec = mean(pnflux_mod_monthly_vecs,1); +ppflux_mod_ann_vec = mean(ppflux_mod_monthly_vecs,1); +siflux_mod_ann_vec = mean(siflux_mod_monthly_vecs,1); +totnflux_mod_ann_vec = mean(totnflux_mod_monthly_vecs,1); +totpflux_mod_ann_vec = mean(totpflux_mod_monthly_vecs,1); +o2flux_mod_ann_vec = mean(o2flux_mod_monthly_vecs,1); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Produce plots to evaluate the mapping % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% scale marker size with the freshwater flux +ms_vec = zeros(size(Q_mod_vec)); +ms_vec(log10(Q_mod_vec) < 0) = 1; +ms_vec((log10(Q_mod_vec) > 0) & (log10(Q_mod_vec) < 1)) = 1; +ms_vec((log10(Q_mod_vec) > 1) & (log10(Q_mod_vec) < 2)) = 1; +ms_vec((log10(Q_mod_vec) > 2) & (log10(Q_mod_vec) < 3)) = 1; +ms_vec(log10(Q_mod_vec) > 3) = 100; + +% DIC, Alk concentrations and DIC:Alk +for m = 1:12 + + figure(1) + clf + + subplot(1,3,1); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dic_mod_monthly_vecs(m,:),ms_vec, ... + dic_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 2500]); + colorbar + titlestr = ['DIC, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(1,3,2); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,alk_mod_monthly_vecs(m,:),ms_vec, ... + alk_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 2500]); + colorbar + titlestr = ['Alk, meq m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(1,3,3) + dic_alk_ratio = dic_mod_monthly_vecs(m,:)./alk_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dic_alk_ratio, ... + ms_vec,dic_alk_ratio,'filled'); + hold on + view(2); + caxis([0.8 1.2]); + colorbar + title('DIC:Alk ratio'); + + pause +end + +% Nitrogen Concentrations +for m = 1:12 + + figure(1) + clf + + subplot(2,2,1); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,din_mod_monthly_vecs(m,:),ms_vec, ... + din_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 100]); + colorbar + titlestr = ['DIN, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,2,2); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,no3_mod_monthly_vecs(m,:),ms_vec, ... + din_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 100]); + colorbar + titlestr = ['no3, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,2,3); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,nh4_mod_monthly_vecs(m,:),ms_vec, ... + nh4_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 20]); + colorbar + titlestr = ['nh4, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,2,4); + no3_din_ratio = no3_mod_monthly_vecs(m,:)./din_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,no3_din_ratio, ... + ms_vec,no3_din_ratio,'filled'); + hold on + view(2); + caxis([0 1.0]); + colorbar + title('NO3:DIN ratio'); + + pause +end + +for m = 1:12 + + subplot(2,3,1); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,din_mod_monthly_vecs(m,:),ms_vec, ... + din_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 100]); + colorbar + titlestr = ['DIN, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,2); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,don_mod_monthly_vecs(m,:),ms_vec, ... + don_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 100]); + colorbar + titlestr = ['DON, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,3); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pn_mod_monthly_vecs(m,:),ms_vec, ... + pn_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 100]); + colorbar + titlestr = ['PN, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,5) + don_din_ratio = don_mod_monthly_vecs(m,:)./din_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,don_din_ratio, ... + ms_vec,don_din_ratio,'filled'); + hold on + view(2); + caxis([0 2]); + colorbar + title('DON:DIN ratio'); + + subplot(2,3,6) + pn_din_ratio = pn_mod_monthly_vecs(m,:)./din_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pn_din_ratio, ... + ms_vec,pn_din_ratio,'filled'); + hold on + view(2); + caxis([0 2]); + colorbar + title('PN:DIN ratio'); + + pause +end + +% Phosphorus Concentrations +for m = 1:12 + + figure(1) + clf + + subplot(2,3,1); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dip_mod_monthly_vecs(m,:),ms_vec, ... + dip_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 3]); + colorbar + titlestr = ['DIP, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,2); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dop_mod_monthly_vecs(m,:),ms_vec, ... + dop_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 3]); + colorbar + titlestr = ['DOP, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,3); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pp_mod_monthly_vecs(m,:),ms_vec, ... + pp_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 3]); + colorbar + titlestr = ['PP, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,5) + dop_dip_ratio = dop_mod_monthly_vecs(m,:)./dip_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dop_dip_ratio, ... + ms_vec,dop_dip_ratio,'filled'); + hold on + view(2); + caxis([0 3]); + colorbar + title('DOP:DIP ratio'); + + subplot(2,3,6) + pp_dip_ratio = pp_mod_monthly_vecs(m,:)./dip_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pp_dip_ratio, ... + ms_vec,pp_dip_ratio,'filled'); + hold on + view(2); + caxis([0 2]); + colorbar + title('PP:DIP ratio'); + + pause +end + +% silica and oxygen concentrations +for m = 1:12 + figure(1) + clf + + subplot(2,1,1); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,si_mod_monthly_vecs(m,:),ms_vec, ... + si_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 200]); + colorbar + titlestr = ['Si, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,1,2); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,o2_mod_monthly_vecs(m,:),ms_vec, ... + o2_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 350]); + colorbar + titlestr = ['o2, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + pause +end + +% Initialize 2D concentration arrays; these are the ones read into MOM6 to +% specify the nutrient concentrations of river inputs. +DIC_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +ALK_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +NO3_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +NH4_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +LDON_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +SLDON_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +SRDON_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +PO4_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +LDOP_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +SLDOP_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +SRDOP_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +NDET_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +PDET_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +SI_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +O2_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); + +% Map concentration vectors onto 2D arrays +temp = zeros(size(lon_mod)); +for m = 1:12 + temp(ind_ro) = dic_mod_monthly_vecs(m,:); + DIC_CONC(m,:,:) = temp; + temp(ind_ro) = alk_mod_monthly_vecs(m,:); + ALK_CONC(m,:,:) = temp; + + temp(ind_ro) = no3_mod_monthly_vecs(m,:); % contains all NEWS DIN + %temp(ind_ro) = din_mod_monthly_vecs(m,:); + NO3_CONC(m,:,:) = temp; + temp(ind_ro) = nh4_mod_monthly_vecs(m,:); + NH4_CONC(m,:,:) = temp; + temp(ind_ro) = don_mod_monthly_vecs(m,:); + LDON_CONC(m,:,:) = frac_ldon*temp; + SLDON_CONC(m,:,:) = frac_sldon*temp; + SRDON_CONC(m,:,:) = frac_srdon*temp; + temp(ind_ro) = pn_mod_monthly_vecs(m,:); + NDET_CONC(m,:,:) = temp; + + temp(ind_ro) = dip_mod_monthly_vecs(m,:); + PO4_CONC(m,:,:) = temp; + temp(ind_ro) = dop_mod_monthly_vecs(m,:); + LDOP_CONC(m,:,:) = frac_ldop*temp; + SLDOP_CONC(m,:,:) = frac_sldop*temp; + SRDOP_CONC(m,:,:) = frac_srdop*temp; + temp(ind_ro) = pp_mod_monthly_vecs(m,:); + PDET_CONC(m,:,:) = temp; + + temp(ind_ro) = si_mod_monthly_vecs(m,:); + SI_CONC(m,:,:) = temp; + temp(ind_ro) = o2_mod_monthly_vecs(m,:); + O2_CONC(m,:,:) = temp; +end + +% MOM6 is taking river values in moles m-3 for other constituents. Change +% for consistency across river constituents +DIC_CONC = DIC_CONC./1e3; +ALK_CONC = ALK_CONC./1e3; +NO3_CONC = NO3_CONC./1e3; +NH4_CONC = NH4_CONC./1e3; +LDON_CONC = LDON_CONC./1e3; +SLDON_CONC = SLDON_CONC./1e3; +SRDON_CONC = SRDON_CONC./1e3; +NDET_CONC = NDET_CONC./1e3; +PO4_CONC = PO4_CONC./1e3; +LDOP_CONC = LDON_CONC./1e3; +SLDOP_CONC = SLDON_CONC./1e3; +SRDOP_CONC = SRDON_CONC./1e3; +PDET_CONC = PDET_CONC./1e3; +SI_CONC = SI_CONC./1e3; +O2_CONC = O2_CONC./1e3; + +% Add iron concentrations - initialize with nitrate and then overwrite +FED_CONC = NO3_CONC; +FEDET_CONC = NO3_CONC; +% 40 nM dissolved iron concentration from De Baar and De Jong + 30nM +% Colloidal and nanoparticle flux as reported in Canfield and Raiswell +FED_CONC(FED_CONC > 0) = const_fed; +FEDET_CONC(FEDET_CONC > 0) = 0.0; + +ms = 8; + +% quick set of plots to ensure the mapping to the model grid was successful +for m = 1:12 + % DIC and alkalinity + figure(1) + clf + subplot(2,1,1); + title('log10(DIC CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(DIC_CONC(m,:)),ms,log10(DIC_CONC(m,:)),'filled'); + caxis([-1 1]); colorbar; + + subplot(2,1,2); + title('log10(ALK CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(ALK_CONC(m,:)),ms,log10(ALK_CONC(m,:)),'filled'); + caxis([-1 1]); colorbar; + + %pause +end + +% Nitrogen +for m = 1:12 + figure(1); + clf + subplot(3,2,1); + title('log10(NO3 CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(NO3_CONC(m,:)),ms,log10(NO3_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + subplot(3,2,2); + title('log10(NH4 CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(NH4_CONC(m,:)),ms,log10(NH4_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + subplot(3,2,3); + title('log10(LDON CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(LDON_CONC(m,:)),ms,log10(LDON_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + subplot(3,2,4); + title('log10(SLDON CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(SLDON_CONC(m,:)),ms,log10(SLDON_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + subplot(3,2,5); + title('log10(SRDON CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(SRDON_CONC(m,:)),ms,log10(SRDON_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + subplot(3,2,6); + title('log10(NDET CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(NDET_CONC(m,:)),ms,log10(NDET_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + %pause +end + +% Phosphorus +for m = 1:12 + figure(1); + clf + subplot(3,2,1); + title('log10(PO4 CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(PO4_CONC(m,:)),ms,log10(PO4_CONC(m,:)),'filled'); + caxis([-4 -2]); colorbar; + + subplot(3,2,2); + title('log10(LDOP CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(LDOP_CONC(m,:)),ms,log10(LDOP_CONC(m,:)),'filled'); + caxis([-4 -2]); colorbar; + + subplot(3,2,3); + title('log10(SLDOP CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(SLDOP_CONC(m,:)),ms,log10(SLDOP_CONC(m,:)),'filled'); + caxis([-4 -2]); colorbar; + + subplot(3,2,4); + title('log10(SRDOP CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(SRDOP_CONC(m,:)),ms,log10(SRDOP_CONC(m,:)),'filled'); + caxis([-4 -2]); colorbar; + + subplot(3,2,5); + title('log10(PDET CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(PDET_CONC(m,:)),ms,log10(PDET_CONC(m,:)),'filled'); + caxis([-4 -2]); colorbar; + + %pause; +end + +% Iron, Silica, Oxygen +for m = 1:12 + figure(1) + clf + subplot(3,2,1); + title('log10(FED CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(FED_CONC(m,:)),ms,log10(FED_CONC(m,:)),'filled'); + caxis([-5 -3]); colorbar; + + subplot(3,2,2); + title('log10(FEDET CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(FEDET_CONC(m,:)),ms,log10(FEDET_CONC(m,:)),'filled'); + caxis([-5 -3]); colorbar; + + subplot(3,2,3); + title('log10(SI CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(SI_CONC(m,:)),ms,log10(SI_CONC(m,:)),'filled'); + caxis([-3 0]); colorbar; + + subplot(3,2,4); + title('log10(O2 CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(O2_CONC(m,:)),ms,log10(O2_CONC(m,:)),'filled'); + caxis([-3 0]); colorbar; + + %pause; +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Save Files % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% option to save matlab file +% save River_DIC_ALK_RC4US_NWA ALK_CONC DIC_CONC + +% Construct netcdf file following format used by nutrient input files to +% MOM6 + +% Make reference dates for standard non-leap year +dates = [1990 1 16 12 0 0; 1990 2 15 0 0 0; 1990 3 16 12 0 0; ... + 1990 4 16 0 0 0; 1990 5 16 12 0 0; 1990 6 16 0 0 0; ... + 1990 7 16 12 0 0; 1990 8 16 12 0 0; 1990 9 16 0 0 0; ... + 1990 10 16 12 0 0; 1990 11 16 0 0 0; 1990 12 16 12 0 0]; +time = datenum(dates) - datenum([1990 1 1 0 0 0]); + +nlat = size(lat_mod,1); +nlon = size(lat_mod,2); + +ncid = netcdf.create(nc_file_name,'CLOBBER'); +dimid0 = netcdf.defDim(ncid,'time',netcdf.getConstant('NC_UNLIMITED')); +dimid1 = netcdf.defDim(ncid,'y',nlat); +dimid2 = netcdf.defDim(ncid,'x',nlon); + +varid0 = netcdf.defVar(ncid,'time','double',dimid0); +netcdf.putAtt(ncid,varid0,'calendar','NOLEAP'); +netcdf.putAtt(ncid,varid0,'calendar_type','NOLEAP'); +netcdf.putAtt(ncid,varid0,'modulo','T'); +netcdf.putAtt(ncid,varid0,'units','days since 1990-1-1 0:00:00'); +netcdf.putAtt(ncid,varid0,'time_origin','01-JAN-1990 00:00:00'); +varid1 = netcdf.defVar(ncid,'y','int',dimid1); +netcdf.putAtt(ncid,varid1,'cartesian_axis','Y'); +varid2 = netcdf.defVar(ncid,'x','int',dimid2); +netcdf.putAtt(ncid,varid2,'cartesian_axis','X'); +varid3 = netcdf.defVar(ncid,'lat','double',[dimid2 dimid1]); +netcdf.putAtt(ncid,varid3,'units','degrees north'); +varid4 = netcdf.defVar(ncid,'lon','double',[dimid2 dimid1]); +netcdf.putAtt(ncid,varid4,'units','degrees east'); +varid5 = netcdf.defVar(ncid,'DIC_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid5,'units','mol m-3'); +netcdf.putAtt(ncid,varid5,'long_name','DIC_CONC'); +varid6 = netcdf.defVar(ncid,'ALK_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid6,'units','mole Eq. m-3'); +netcdf.putAtt(ncid,varid6,'long_name','ALK_CONC'); +varid7 = netcdf.defVar(ncid,'NO3_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid7,'units','mol m-3'); +netcdf.putAtt(ncid,varid7,'long_name','NO3_CONC'); +varid8 = netcdf.defVar(ncid,'NH4_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid8,'units','mol m-3'); +netcdf.putAtt(ncid,varid8,'long_name','NH4_CONC'); +varid9 = netcdf.defVar(ncid,'LDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid9,'units','mol m-3'); +netcdf.putAtt(ncid,varid9,'long_name','0.3*DON_CONC'); +varid10 = netcdf.defVar(ncid,'SLDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid10,'units','mol m-3'); +netcdf.putAtt(ncid,varid10,'long_name','0.35*DON_CONC'); +varid11 = netcdf.defVar(ncid,'SRDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid11,'units','mol m-3'); +netcdf.putAtt(ncid,varid11,'long_name','0.35*DON_CONC'); +varid12 = netcdf.defVar(ncid,'NDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid12,'units','mol m-3'); +netcdf.putAtt(ncid,varid12,'long_name','1.0*PN_CONC'); +varid13 = netcdf.defVar(ncid,'PO4_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid13,'units','mol m-3'); +netcdf.putAtt(ncid,varid13,'long_name','PO4_CONC+0.3*PP_CONC'); +varid14 = netcdf.defVar(ncid,'LDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid14,'units','mol m-3'); +netcdf.putAtt(ncid,varid14,'long_name','0.3*DOP_CONC'); +varid15 = netcdf.defVar(ncid,'SLDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid15,'units','mol m-3'); +netcdf.putAtt(ncid,varid15,'long_name','0.35*DOP_CONC'); +varid16 = netcdf.defVar(ncid,'SRDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid16,'units','mol m-3'); +netcdf.putAtt(ncid,varid16,'long_name','0.35*DOP_CONC'); +varid17 = netcdf.defVar(ncid,'PDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid17,'units','mol m-3'); +netcdf.putAtt(ncid,varid17,'long_name','0*PP_CONC'); +varid18 = netcdf.defVar(ncid,'FED_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid18,'units','mol m-3'); +netcdf.putAtt(ncid,varid18,'long_name','FED_CONC'); +varid19 = netcdf.defVar(ncid,'FEDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid19,'units','mol m-3'); +netcdf.putAtt(ncid,varid19,'long_name','FEDET_CONC'); +varid20 = netcdf.defVar(ncid,'O2_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid20,'units','mol m-3'); +netcdf.putAtt(ncid,varid20,'long_name','O2_CONC'); +netcdf.close(ncid) + +ncid = netcdf.open(nc_file_name,'NC_WRITE'); +netcdf.putVar(ncid,varid0,0,12,time); +% nutrient input files appear seem to need dummy axes to be read in +% properly, but eventually do a grid by grid mapping that doesn't require +% these. +netcdf.putVar(ncid,varid1,1:nlat); +netcdf.putVar(ncid,varid2,1:nlon); +netcdf.putVar(ncid,varid3,permute(lat_mod,[2,1])); +netcdf.putVar(ncid,varid4,permute(lon_mod,[2,1])); +netcdf.putVar(ncid,varid5,permute(DIC_CONC,[3,2,1])); +netcdf.putVar(ncid,varid6,permute(ALK_CONC,[3,2,1])); +netcdf.putVar(ncid,varid7,permute(NO3_CONC,[3,2,1])); +netcdf.putVar(ncid,varid8,permute(NH4_CONC,[3,2,1])); +netcdf.putVar(ncid,varid9,permute(LDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid10,permute(SLDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid11,permute(SRDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid12,permute(NDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid13,permute(PO4_CONC,[3,2,1])); +netcdf.putVar(ncid,varid14,permute(LDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid15,permute(SLDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid16,permute(SRDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid17,permute(PDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid18,permute(FED_CONC,[3,2,1])); +netcdf.putVar(ncid,varid19,permute(FEDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid20,permute(O2_CONC,[3,2,1])); +netcdf.close(ncid) diff --git a/tools/rivers/bgc/NEP/Data/ArcticGro/ArcticGRO_Water_Quality_Data.xlsx b/tools/rivers/bgc/NEP/Data/ArcticGro/ArcticGRO_Water_Quality_Data.xlsx new file mode 100644 index 000000000..9a60d7aec Binary files /dev/null and b/tools/rivers/bgc/NEP/Data/ArcticGro/ArcticGRO_Water_Quality_Data.xlsx differ diff --git a/tools/rivers/bgc/NEP/Data/ArcticGro/ArcticGro_Process.m b/tools/rivers/bgc/NEP/Data/ArcticGro/ArcticGro_Process.m new file mode 100644 index 000000000..fc9f33be8 --- /dev/null +++ b/tools/rivers/bgc/NEP/Data/ArcticGro/ArcticGro_Process.m @@ -0,0 +1,129 @@ +% Program to process the ArcticGro data + +clear all + +filename = 'ArcticGRO_Water_Quality_Data.xlsx' +% co2 system solver to derive DIC from alkalinity, pH and temperature +addpath /home/cas/matlab/co2sys + +% Longitude and Latitude from GlobalNEWS +river_name{1} = 'Ob'; +lon(1) = 68.5; lat(1) = 68.75; +river_name{2} = 'Yenisey'; +lon(2) = 82.25; lat(2) = 71.25; +river_name{3} = 'Lena'; +lon(3) = 128; lat(3) = 73; +river_name{4} = 'Kolyma'; +lon(4) = 161.25; lat(4) = 69.25; +river_name{5} = 'Yukon'; +lon(5) = -164.75; lat(5) = 62.75; +river_name{6} = 'Mackenzie'; +lon(6) = -134.75; lat(6) = 69.25; + + +for n = 1:6 + Data = readtable(filename,'Sheet',n); + discharge(n) = nanmean(Data{:,5}); + temp(n) = nanmean(Data{:,6}); + % alk in mg CaCO3/L ~ g CaCO3/m3 + alk(n) = nanmean(Data{:,9}); + alk(n) = alk(n)/(40+12+16*3)*2*1e3; % to milliequivalents per m-3 + % pH + pH(n) = nanmean(Data{:,7}); + % tdn in mg N L-1 ~ g N/m3 + tdn(n) = nanmean(Data{:,21}); + tdn(n) = tdn(n)*1e3/14; % mmoles m-3 + % no3 in micrograms N per L ~ mg N m-3 + no3(n) = nanmean(Data{:,22}); + no3(n) = no3(n)/14; % mmoles m-3 + % nh4 in micrograms N per L ~ mg N m-3 + nh4(n) = nanmean(Data{:,23}); + nh4(n) = nh4(n)/14; % mmoles m-3 + % tdp in micrograms P per L ~ mg P m-3 + tdp(n) = nanmean(Data{:,24}); + tdp(n) = tdp(n)/31; % mmoles m-3 + % po4 in micrograms P per L ~ mg P m-3 + po4(n) = nanmean(Data{:,25}); + po4(n) = po4(n)/31; % mmoles m-3 + % sio2 in mg SiO2 per L ~ g P m-3 + sio2(n) = nanmean(Data{:,26}); + sio2(n) = sio2(n)*1e3/(28.06+16*2); + % pon in micrograms N per L ~ mg N m-3 + pon(n) = nanmean(Data{:,47}); + pon(n) = pon(n)/14; %mmoles m-3 + + % calculate DIC from alk, pH and DIC + out = CO2SYS(alk(n),pH(n),1,3,0,temp(n),temp(n), ... + 100,100,0,0,0,0,4,15,1,2,2); + dic(n) = mean(out(:,2)); + + % calculate don from tdn, no3 and nh4 + aa = find(isfinite(Data{:,21}) & isfinite(Data{:,22}) & isfinite(Data{:,23})); + don(n) = nanmean(Data{aa,21}*1e3/14 - Data{aa,22}/14 - Data{aa,23}/14); + + % calculate dop from tdp and po4 + bb = find(isfinite(Data{:,24}/31) & isfinite(Data{:,25}/31)); + dop(n) = nanmean(Data{aa,24}/31 - Data{aa,25}/31); +end + +% Fill in particulate phosphorus from GLOBAL NEWS +pp(1) = 1.29; % Ob +pp(2) = 0.82; % Yenisey +pp(3) = 1.48; % Lena +pp(4) = 1.21; % Kolyma +pp(5) = 1.94; % Yukon +pp(6) = 1.44; % Mackenzie + +% For reference (DIN, DON, PN, DIP, DOP, PP) +% Yukon: 7.30 17.25 29.10 0.07 0.42 1.94 +% Mackenzie: 7.42 19.86 24.34 0.14 0.48 1.44 +% St. Lawrence: 52.80 24.86 3.81 0.98 0.55 0.21 +% Ob: 21.80 24.16 22.35 1.16 0.56 1.29 +% Lena: 7.74 21.34 25.48 0.24 0.52 1.48 +% Yenisey: 8.54 21.65 14.21 0.088 0.52 0.82 +% Kolyma: 10.17 21.92 21.27 0.18 0.53 1.21 + +lon_stations_arcticgro = lon; +lat_stations_arcticgro = lat; +station_names_arcticgro = river_name; + +Q_ann_arcticgro = discharge; +dic_ann_arcticgro = dic; +alk_ann_arcticgro = alk; +no3_ann_arcticgro = no3; +nh4_ann_arcticgro = nh4; +din_ann_arcticgro = no3_ann_arcticgro + nh4_ann_arcticgro; +pn_ann_arcticgro = pon; +don_ann_arcticgro = don; +dip_ann_arcticgro = po4; +dop_ann_arcticgro = dop; +pp_ann_arcticgro = pp; +si_ann_arcticgro = sio2; +o2_ann_arcticgro = ones(size(Q_ann_arcticgro))*NaN; +dfe_ann_arcticgro = ones(size(Q_ann_arcticgro))*NaN; +pfe_ann_arcticgro = ones(size(Q_ann_arcticgro))*NaN; + +% did not try and extract monthly data from arcticgro, so leave as NaNs +dic_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +alk_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +no3_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +nh4_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +din_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +don_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +pn_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +dip_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +dop_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +pp_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +dfe_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +pfe_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +si_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; +o2_monthly_arcticgro = ones(12,size(lon_stations_arcticgro,2))*NaN; + +save arcticgro_data lon_stations_arcticgro lat_stations_arcticgro station_names_arcticgro ... + Q_ann_arcticgro dic_ann_arcticgro alk_ann_arcticgro no3_ann_arcticgro nh4_ann_arcticgro ... + din_ann_arcticgro pn_ann_arcticgro don_ann_arcticgro dip_ann_arcticgro dop_ann_arcticgro ... + pp_ann_arcticgro si_ann_arcticgro o2_ann_arcticgro dfe_ann_arcticgro pfe_ann_arcticgro ... + dic_monthly_arcticgro alk_monthly_arcticgro no3_monthly_arcticgro nh4_monthly_arcticgro ... + din_monthly_arcticgro don_monthly_arcticgro pn_monthly_arcticgro dip_monthly_arcticgro ... + dop_monthly_arcticgro pp_monthly_arcticgro dfe_monthly_arcticgro pfe_monthly_arcticgro ... + si_monthly_arcticgro o2_monthly_arcticgro diff --git a/tools/rivers/bgc/NEP/Data/GLORICH/NEP_GLORICH_Process.m b/tools/rivers/bgc/NEP/Data/GLORICH/NEP_GLORICH_Process.m new file mode 100644 index 000000000..76c4186e1 --- /dev/null +++ b/tools/rivers/bgc/NEP/Data/GLORICH/NEP_GLORICH_Process.m @@ -0,0 +1,288 @@ +clear all; + +% load in matlab file with GLORICH data +% created with:T = readtable('hydrochemistry.csv','NumHeaderLines',1); +%load GLORICH_Data.mat; +T = readtable('hydrochemistry.csv','NumHeaderLines',1); +stations = T{:,1}; + +%"102481","COLORADO R AT NIB AB MORELOS DAM NR ANDRADE, CA.","9522000","USA","AZ",32.71,-114.71,"NA1983" +%"110042","Fraser River at Hope","BC08MF0001","Canada","BC",49.38,-121.45,"NA1983 +%"110009","Skeena River at Usk","BC08EF0001","Canada","BC",54.63,-128.41,"NA1983" +%"110004","Stikine River above Choquette River","BC08CF0002","Canada","BC",56.82,-131.76,"NA1983" +%"102954","COPPER R AT MILLION DOLLAR BRIDGE NR CORDOVA AK","15214000","USA","AK",60.67,-144.74,"NA1983" +%"102983","SUSITNA R AT SUSITNA STATION AK","15294350","USA","AK",61.54,-150.51,"NA1983" +%"102985","NUSHAGAK R AT EKWOK AK","15302500","USA","AK",59.34,-157.47,"NA1983" +%"102986","KUSKOKWIM R AT CROOKED CREEK AK","15304000","USA","AK",61.87,-158.10,"NA1983" + +% Define rivers. I took the lat/lon from the mouths defined in Global NEWS +% Where stream gages are upstream, this makes the assumption that those +% properties are reasonably indicative of conditions at the river mouth. +num_rivers = 8; +river_name{1} = 'Colorado'; station_num(1) = 102481; +lat(1) = 31.75; lon(1) = 245.25; +river_name{2} = 'Fraser'; station_num(2) = 110042; +lat(2) = 49.25; lon(2) = 237.25; +river_name{3} = 'Skeena'; station_num(3) = 110009; +lat(3) = 54.25; lon(3) = 229.75; +river_name{4} = 'Stikine'; station_num(4) = 110004; +lat(4) = 56.75; lon(4) = 227.75; +river_name{5} = 'Copper'; station_num(5) = 102954; +lat(5) = 60.25; lon(5) = 215.25; +river_name{6} = 'Susitna'; station_num(6) = 102983; +lat(6) = 61.75; lon(6) = 209.75; +river_name{7} = 'Nushagak'; station_num(7) = 102985; +lat(7) = 58.75; lon(7) = 201.75; +river_name{8} = 'Kuskokwim'; station_num(8) = 102986; +lat(8) = 60.25; lon(8) = 197.75; +% Rely on ArcticGRO for the Arctic rivers +% Yukon at Pilot Station +%river_name{9} = 'Yukon'; station_num(9) = 102988; +%lat(9) = 60.25; lon(9) = 197.75; +% Mackenzie at Tsiigehtchic is the same as ArcticGRO +%river_name{10} = 'Mackenzie'; station_num(10) = 120036; +%lat(10) = 60.25; lon(10) = 197.75; +%river_name{9} = 'Anadyr'; station_num(9) = 510090; +%lat(9) = 64.75; lon(9) = 177.75; + +% co2 system solver to derive DIC from alkalinity, pH and temperature +addpath /home/cas/matlab/co2sys + +for n = 1:num_rivers + aa = find(stations == station_num(n)); + num_stations(n) = size(aa,1); + discharge(n) = nanmean(T{aa,6}); + alk_vec = T{aa,20}; + alk(n) = nanmean(alk_vec); + % nitrogen species + tn_vec = T{aa,62}; num_tn(n) = size(find(isfinite(tn_vec) == 1),1); tn(n) = nanmean(tn_vec); + dn_vec = T{aa,64}; num_dn(n) = size(find(isfinite(dn_vec) == 1),1); dn(n) = nanmean(dn_vec); + % first method of extracting particulate nitrogen from GLORICH + pn1_vec = tn_vec - dn_vec; pn1_vec(pn1_vec < 0) = 0; + num_pn1(n) = size(find(isfinite(pn1_vec) == 1),1); pn1(n) = nanmean(pn1_vec); + %pn_vec = T{aa,66}; num_pn(n) = size(find(isfinite(pn_vec) == 1),1); pn(n) = nanmean(pn_vec); + tin_vec = T{aa,68}; num_tin(n) = size(find(isfinite(tin_vec) == 1),1); tin(n) = nanmean(tin_vec); + din_vec = T{aa,70}; num_din(n) = size(find(isfinite(din_vec) == 1),1); din(n) = nanmean(din_vec); + ton_vec = T{aa,72}; num_ton(n) = size(find(isfinite(ton_vec) == 1),1); ton(n) = nanmean(ton_vec); + %don_vec = T{aa,74}; num_don(n) = size(find(isfinite(don_vec) == 1),1); don(n) = nanmean(don_vec); + pon_vec = T{aa,76}; num_pon(n) = size(find(isfinite(pon_vec) == 1),1); pon(n) = nanmean(pon_vec); + tkn_vec = T{aa,78}; num_tkn(n) = size(find(isfinite(tkn_vec) == 1),1); tkn(n) = nanmean(tkn_vec); + dkn_vec = T{aa,80}; num_dkn(n) = size(find(isfinite(dkn_vec) == 1),1); dkn(n) = nanmean(dkn_vec); + % second method of extracting particulate nitrogen from GLORICH + pn2_vec = tkn_vec - dkn_vec; pn1_vec(pn2_vec < 0) = 0; + num_pn2(n) = size(find(isfinite(pn2_vec) == 1),1); pn2(n) = nanmean(pn2_vec); + no3_vec = T{aa,82}; num_no3(n) = size(find(isfinite(no3_vec) == 1),1); no3(n) = nanmean(no3_vec); + no2_vec = T{aa,84}; num_no2(n) = size(find(isfinite(no2_vec) == 1),1); no2(n) = nanmean(no2_vec); + no2no3_vec = T{aa,86}; num_no2no3(n) = size(find(isfinite(no2no3_vec) == 1),1); no2no3(n) = nanmean(no2no3_vec); + tnh4_vec = T{aa,88}; num_tnh4(n) = size(find(isfinite(tnh4_vec) == 1),1); tnh4(n) = nanmean(tnh4_vec); + dnh4_vec = T{aa,90}; num_dnh4(n) = size(find(isfinite(dnh4_vec) == 1),1); dnh4(n) = nanmean(dnh4_vec); + don_vec = dn_vec-no2no3_vec-dnh4(n); don_vec(don_vec < 0) = 0; + num_don(n) = size(find(isfinite(don_vec) == 1),1); don(n) = nanmean(don_vec); + % phosphorus + tp_vec = T{aa,92}; num_tp(n) = size(find(isfinite(tp_vec) == 1),1); tp(n) = nanmean(tp_vec); + dp_vec = T{aa,94}; num_dp(n) = size(find(isfinite(dp_vec) == 1),1); dp(n) = nanmean(dp_vec); + pp_vec = tp_vec - dp_vec; pp_vec(pp_vec < 0) = 0; + num_pp(n) = size(find(isfinite(pp_vec) == 1),1); pp(n) = nanmean(pp_vec); + %pp_vec = T{aa,96}; pp(n) = nanmean(pp_vec); + tip_vec = T{aa,98}; num_tip(n) = size(find(isfinite(tip_vec) == 1),1); tip(n) = nanmean(tip_vec); + dip_vec = T{aa,100}; num_dip(n) = size(find(isfinite(dip_vec) == 1),1); dip(n) = nanmean(dip_vec); + dop_vec = dp_vec - dip_vec; dop_vec(dop_vec < 0) = 0; + num_dop(n) = size(find(isfinite(dop_vec) == 1),1); dop(n) = nanmean(dop_vec); + % oxygen + o2_vec = T{aa,12}; num_o2(n) = size(find(isfinite(o2_vec) == 1),1); o2(n) = nanmean(o2_vec); + % silica + si_vec = T{aa,34}; num_si(n) = size(find(isfinite(si_vec) == 1),1); si(n) = nanmean(si_vec); + % dic, fill in with pH if needed + %dic_vec = T{aa,52}; + %bb = find(isfinite(dic_vec)); + pH_vec = T{aa,10}; pH(n) = nanmean(pH_vec); + temp_vec = T{aa,8}; temp(n) = nanmean(temp_vec); + + + % Simple initial DIC calculation based on mean alk and mean pH, room for + % improvement here. + + if isnan(temp(n)); temp(n) = 5; end; + %[RESULT,HEADERS,NICEHEADERS]=CO2SYS(PAR1,PAR2,PAR1TYPE,PAR2TYPE,... + % ...SAL,TEMPIN,TEMPOUT,PRESIN,PRESOUT,SI,PO4,NH4,H2S,pHSCALEIN,... + % ...K1K2CONSTANTS,KSO4CONSTANT,KFCONSTANT,BORON) + out = CO2SYS(alk(n),pH(n),1,3,0,temp(n),temp(n), ... + 100,100,0,0,0,0,4,15,1,2,2); + dic(n) = mean(out(:,2)); +end + +lon_stations_glorich = lon; +lat_stations_glorich = lat; +station_names_glorich = river_name; + +Q_ann_glorich = discharge; +dic_ann_glorich = dic; +alk_ann_glorich = alk; +no3_ann_glorich = no2no3; +% use no3 for the Stikine +no3_ann_glorich(4) = no3(4); +nh4_ann_glorich = dnh4; +% assume nh4 in Stikine is equal to nh4 in Skeena +nh4_ann_glorich(4) = nh4_ann_glorich(3); +din_ann_glorich = no3_ann_glorich + nh4_ann_glorich; +% Use weighted average of TN-DN and TKN-DKN +pn1(num_pn1 == 0) = 0; pn2(num_pn2 == 0) = 0; +pn_ann_glorich = (pn1.*num_pn1 + pn2.*num_pn2)./(num_pn1 + num_pn2); +% assume particulate nitrogen in Stikine is equal to Skeena +pn_ann_glorich(4) = pn_ann_glorich(3); +don_ann_glorich = don; +% assume dissolved organic nitrogen in Stikine is equal to Skeena +don_ann_glorich(4) = don_ann_glorich(3); + +dip_ann_glorich = dip; +% assume po4 in Stikine is equal to po4 in GLORICH +dip_ann_glorich(4) = dip_ann_glorich(3); + +dop_ann_glorich = dop; +% assume dop in the Stikine and Skeena are equal to that in the Copper +% DOP constraints are generally poor. The Fraser is based on a single +% joint dp, dip measurement. However, the contribution of dop to the total +% phosphorus load is also very low, so this is an issue I'm willing to +% stomach. +dop_ann_glorich(3) = dop_ann_glorich(5); +dop_ann_glorich(4) = dop_ann_glorich(5); + +pp_ann_glorich = pp; +% set particulate load in the Stikine to that in the Skeena +pp_ann_glorich(4) = pp_ann_glorich(3); + +si_ann_glorich = si; +o2_ann_glorich = o2; +dfe_ann_glorich = ones(size(o2_ann_glorich))*NaN; +pfe_ann_glorich = ones(size(o2_ann_glorich))*NaN; + +% did not try and extract monthly data from GLORICH, so leave as NaNs +dic_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +alk_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +no3_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +nh4_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +din_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +don_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +pn_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +dip_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +dop_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +pp_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +dfe_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +pfe_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +si_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; +o2_monthly_glorich = ones(12,size(lon_stations_glorich,2))*NaN; + +save NEP_GLORICH_data lon_stations_glorich lat_stations_glorich station_names_glorich ... + Q_ann_glorich dic_ann_glorich alk_ann_glorich no3_ann_glorich nh4_ann_glorich ... + din_ann_glorich pn_ann_glorich don_ann_glorich dip_ann_glorich dop_ann_glorich ... + pp_ann_glorich si_ann_glorich o2_ann_glorich dfe_ann_glorich pfe_ann_glorich ... + dic_monthly_glorich alk_monthly_glorich no3_monthly_glorich nh4_monthly_glorich ... + din_monthly_glorich don_monthly_glorich pn_monthly_glorich dip_monthly_glorich ... + dop_monthly_glorich pp_monthly_glorich dfe_monthly_glorich pfe_monthly_glorich ... + si_monthly_glorich o2_monthly_glorich + + +% Ordering of the GLORICH Data +%1. "STAT_ID", +%2. "RESULT_DATETIME", +%3. "SAMPLE_TIME_DESC", +%4. "SAMPLING_MODE", +%5. "Ref", +%6. "Discharge_inst", +%7. "Discharge_inst_vrc", +%8. "Temp_water", +%9. "Temp_water_vrc", +%10. "pH", +%11. "pH_vrc", +%12. "DO_mgL", +%13. "DO_mgL_vrc", +%14. "DOSAT", +%15. "DOSAT_vrc", +%16. "SpecCond25C", +%17. "SpecCond25C_vrc", +%18. "SPM", +%19. "SPM_vrc", +%20. "Alkalinity", +%21. "Alkalinity_vrc", +%22. "HCO3", +%23. "HCO3_vrc", +%24. "CO3", +%25. "CO3_vrc", +%26. "Ca", +%27. "Ca_vrc", +%28. "Mg", +%29. "Mg_vrc", +%30. "Na", +%31. "Na_vrc", +%32. "K", +%33. "K_vrc", +%34. "SiO2", +%35. "SiO2_vrc", +%36. "Cl", +%37. "Cl_vrc", +%38. "SO4", +%39. "SO4_vrc", +%40. "F", +%41. "F_vrc", +%42. "DSr", +%43. "DSr_vrc", +%44. "TC", +%45. "TC_vrc", +%46. "DC", +%47. "DC_vrc", +%48. "PC", +%49. "PC_vrc", +%50. "TIC", +%51. "TIC_vrc", +%52. "DIC", +%53. "DIC_vrc", +%54. "PIC", +%55. "PIC_vrc", +%56. "TOC", +%57. "TOC_vrc", +%58. "DOC", +%59. "DOC_vrc", +%60. "POC", +%61. "POC_vrc", +%62. "TN", +%63. "TN_vrc", +%64. "DN", +%65. "DN_vrc", +%66. "PN", +%67. "PN_vrc", +%68. "TIN", +%69. "TIN_vrc", +%70. "DIN", +%71. "DIN_vrc", +%72. "TON", +%73. "TON_vrc", +%74. "DON", +%75. "DON_vrc", +%76. "PON", +%77. "PON_vrc", +%78. "TKN", +%79. "TKN_vrc", +%80. "DKN", +%81. "DKN_vrc", +%82. "NO3", +%83. "NO3_vrc", +%84. "NO2", +%85. "NO2_vrc", +%86. "NO2_NO3", +%87. "NO2_NO3_vrc", +%88. "TNH4", +%89. "TNH4_vrc", +%90. "DNH4", +%91. "DNH4_vrc", +%92. "TP", +%93. "TP_vrc", +%94. "DP", +%95. "DP_vrc", +%96. "PP", +%97. "PP_vrc", +%98. "TIP", +%99. "TIP_vrc", +%100. "DIP", +%101. "DIP_vrc", +%102. "PS", +%103. "PS_vrc" diff --git a/tools/rivers/bgc/NEP/README.md b/tools/rivers/bgc/NEP/README.md new file mode 100644 index 000000000..98f02bb6f --- /dev/null +++ b/tools/rivers/bgc/NEP/README.md @@ -0,0 +1,112 @@ +## Example Scripts for NEP BGC runoff generation + +This folder contains example scripts for NEP BGC runoff file generation. Users can follow the following instructions to generate BGC runoff file: + +1: Generate a monthly climatology of the river inputs on the model grid +using `make_discharge_climatology.m`. This routine creates a monthly climatology +from the model's freshwater forcing. For example, the default uses GLOFAS/HILL +forcing files created by Kate Hedstrom on 5/11/2023. All you need to do to +regenerate files (or create new ones) is update the file path. The file creates +a matlab file (`*.mat`) with the discharge climatology. This is used in the +assignment of rivers to discharge points. +``` +matlab232 -nodisplay -nosplash -nodesktop -r "run('make_discharge_climatology_nep.m');exit;" +``` + +2: Use `mapriv_NEWS2.m` to create a river nutrient input file based on the +GlobalNEWS2 estimates (Mayorga et al., 2010). GlobalNEWS contains a +database of global rivers with empirically-derived nutrient inputs. GlobalNEWS +does not, however, have DIC or alkalinity so it cannot be used to provide forcing for +carbon cycling. Also, while globalNEWS is quite skilfull globally, it can +have significant regional biases. Nonetheless, the routine provides a good +way to identify major rivers in the region, and the comprehensive nutrient +estimates that it provides will eventually be used to fill in gaps in river +forcing drawn from direct observations. + +More details about the mapping algorithm can be found below (and in the code). +The first time through, I would recommend setting `inspect_map1` to `y`. +This allows you to step through the mapping of each river and assess its +quality and properties. I have also found that applying a minimum discharge +of 100 m3 sec-1 helps avoid erroneous mapping of very small rivers onto +large discharges. The assignment algorithm was designed to be relatively +insensitive to such instances, but care should still be taken. + +The required inputs are the discharge climatology (from step 1) and a copy +of the globalNEWS data. +``` +matlab232 -nodisplay -nosplash -nodesktop -r "run('mapriv_NEWS2.m');exit;" +``` + +3: Gather/Process direct river measurements: The next step is to get as many +direct measurements as you can. + + - [RC4USCoast](https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.nodc:0260455) +``` +mkdir -p Data/RC4USCoast +cd Data/RC4USCoast +wget https://www.ncei.noaa.gov/archive/archive-management-system/OAS/bin/prd/jquery/download/260455.3.3.tar.gz +tar -xzf 260455.3.3.tar.gz --wildcards --strip-components=4 -C . "*/data/0-data/*.nc" +``` + + - [GLORICH](https://www.geo.uni-hamburg.de/en/geologie/forschung/aquatische-geochemie/glorich.html) +``` +cd Data/GLORICH +wget https://store.pangaea.de/Publications/HartmannJens-etal_2019/Glorich_V01_CSV_plus_Shapefiles_2019_05_24.zip +unzip Glorich_V01_CSV_plus_Shapefiles_2019_05_24.zip +matlab232 -nodisplay -nosplash -nodesktop -r "run('NEP_GLORICH_Process.m');exit;" +``` + + - [ArcticGro](https://arcticgreatrivers.org/data/) +``` +# A copy of this dataset is provided in the Data directory. +# You may want to download your own to ensure +# that it is up-to-date. +cd Data/ArcticGro +matlab232 -nodisplay -nosplash -nodesktop -r "run('ArcticGro_Process.m');exit;" +``` +4: run `mapriv_combined_NEP10k` to create river nutrient and carbon input +estimates based on available observation, while using GlobalNEWS to fill in +some gaps. +``` +cp /archive/ynt//woa_sst_climo.nc ./Data/ +matlab232 -nodisplay -nosplash -nodesktop -r "run('mapriv_combined_NEP10k.m');exit;" +``` +As was the case for globalNEWS2, I recommend running these with "inspect_map = 'y'" until +you are satisfied with the results. You can just "click through" each river mapping and +confirm its properties. The routine will then produce numerous final plots for analysis +and quality control, before generating the netcdf for use with MOM6. + +Note: River oxygen levels are set at saturating levels at temperatures provided by the +World Ocean Atlas Climatology. +___________________________________________________________________________________________ + +## RIVER MAPPING ALGORITHM: +Rivers are first filtered to isolate those within the domain and +above a user specified flow threshold. The rationale for this threshold is to minimize the +risk of inadvertantly mapping very small rivers onto large freshwater flows. Small rivers +tend to have more volatile nutrient concentrations, so such mistakes can have large consequences. +To further reduce this risk, rivers are then sorted by size from smallest to largest. Model +outflow points nearest to each river are accumulated until the the value closest to the observed +flow is reached. These points are assigned the concentration for that river. If a larger river +subsequently claims those points, the larger river is given precedence. The concentrations of +any model discharge points left unassigned after all estimates have been cycled through +are assigned using a nearest neighbor algorithm. + +## Note: +One must specify the fraction of particulate phosphorus that is bioavailable (generally between +10-30%, Froelich, Kinetic control of dissolved phosphate in natural rivers and estuaries: A primer +on the phosphate buffer mechanism, Limnology and Oceanography), and the fraction of dissolved organic +inputs that fall into labile, semi-labile and semi-refractory pools. This is set by default to 30%, 35%, +and 35%. This is consistent with the range of bio-availability suggested by Weigner et al. (2006), +Bio-availability of dissolved organic nitrogen and carbon from 9 rivers in the eastern US, Aquatic +Microbial Ecology. + +## VISUALIZATION AND QUALITY CONTROL: +The routines include a y/n toggle option for inspecting the +mapping of each river as it is done. If "inspect_map" is set to 'y', a graphical map of the +model discharge point assigned to each river is presented, along with the outflow and nutrient +characteristics of the river. This can be useful for identifying rivers that may require some +manual editing to ensure the outflow is assigned to the right place. A section for such manual +edits is provided in the code. The visualization pauses until the user presses any key. It +then moves onto the next river. Once all the rivers have been mapped and interpolation completed, +the routine always produces a series of domain-wide plots to evaluate the overall results. diff --git a/tools/rivers/bgc/NEP/make_discharge_climatology_nep.m b/tools/rivers/bgc/NEP/make_discharge_climatology_nep.m new file mode 100644 index 000000000..9537a15d3 --- /dev/null +++ b/tools/rivers/bgc/NEP/make_discharge_climatology_nep.m @@ -0,0 +1,61 @@ +clear all; +addpath /home/cas/matlab +nc64startup; + +syear = 1993; +eyear = 2019; +num_years = eyear-syear+1; +readme = '1993-2019 monthly NEP10k climatology from GLOFAS/Hill (May 12, 2023 from Liz Drenkard), kg m-2 sec-1'; + +% get grid information +file_name = ['/archive/cas/Regional_MOM6/NEP10k/glofas_hill_05122023/glofas_hill_dis_runoff_mac_',num2str(syear,'%4u'),'.nc'] +lon = ncread(file_name,'lon'); +lat = ncread(file_name,'lat'); +nlat = size(lat,2); +nlon = size(lon,1); +area = ncread(file_name,'area'); + +% holds the runoff climatology +runoff_mc = zeros(nlon,nlat,12); + +for yr = syear:eyear + file_name = ['/archive/cas/Regional_MOM6/NEP10k/glofas_hill_05122023/glofas_hill_dis_runoff_mac_',num2str(yr,'%4u'),'.nc'] + + time = ncread(file_name,'time'); + ntime = size(time,1); + + date = datevec(time+datenum(1950,1,1,0,0,0)); + runoff_days_temp = ncread(file_name,'runoff'); + % files are padded with the first day of the following year, remove + runoff_days = runoff_days_temp(:,:,1:(ntime-1)); + month = date(1:(ntime-1),2); + + for m = 1:12 + aa = find(month == m); + runoff_temp = runoff_days(:,:,aa); + runoff_mc(:,:,m) = runoff_mc(:,:,m) + mean(runoff_temp,3)/num_years; + end + + clear runoff_days runoff_days_temp; + clear runoff_temp aa time date month; + +end + +% modify so that it is time, nlat, nlon +runoff = permute(runoff_mc,[3 2 1]); +area = permute(area,[2 1]); +lon = permute(lon,[2 1]); +lat = permute(lat,[2 1]); + +for m = 1:12; temp = squeeze(runoff(m,:,:)); total_runoff(m) = sum(temp(:).*area(:)); end +figure(1); clf; plot(total_runoff); xlabel('month'); ylabel('runoff, kg sec-1'); + +for m = 1:12; + figure(2); + clf + temp = squeeze(runoff(m,:,:)); + scatter3(lon(:),lat(:),log10(temp(:)),ones(size(temp(:)))*10,log10(temp(:))); view(2); caxis([-3 -1]); colorbar; + pause +end + +save glofas_hill_runoff_monthlyclim_NEP10k_05122023 runoff area lat lon readme; diff --git a/tools/rivers/bgc/NEP/mapriv_NEWS2.m b/tools/rivers/bgc/NEP/mapriv_NEWS2.m new file mode 100644 index 000000000..23bc481a2 --- /dev/null +++ b/tools/rivers/bgc/NEP/mapriv_NEWS2.m @@ -0,0 +1,859 @@ +% Routine to map Global NEWS nutrient data onto the MOM6 Arctic grid @ + +clear all; +addpath /home/cas/matlab +nc64startup; + +% name of netcdf file to be created +nc_file_name = 'RiverNutrients_GlobalNEWS2_plusFe_Q100_NEP10k.nc'; +% Arctic +%nc_file_name = 'RiverNutrients_GlobalNEWS2_plusFe_Q100_Arctic12k.nc'; + +% Parameters for the assignment algorithm. +Q_min = 100; % minimum flow in m3 sec +plot_width = 15; % width of window (in degrees) for inspecting locations + % of rivers and outflow points that have been assigned to + % them. +min_dist = 2.5; % minimum distance (degrees) of the closest outflow point + % for the river to be considered in the domain (useful + % for preventing the algorithm from trying to map rivers + % flowing to different ocean basins. +max_dist = 3; % maximum distance to search for a point +inspect_map = 'n'; % flag enabling you to pause and inspect each river + % mapping as it is being done. + +% set the bio-availability of phosphorus and the fractionation of dissolved +% organic; The code assumes that 30% of particulate phosphorus (PP) is +% bioavailable (e.g., Frolich et al., 1988). One also needs to divide the +% dissolved organic nitrogen and phosphorus components into fractions with +% different lability. The default values are based losely on Wiegner et +% al., (2006). If you have better information, feel free to modify. +frac_PP = 0.3; +frac_ldon = 0.3; +frac_sldon = 0.35; +frac_srdon = 0.35; +frac_ldop = 0.3; +frac_sldop = 0.35; +frac_srdop = 0.35; +% 40 nM dissolved iron concentration from De Baar and De Jong + 30nM +% Colloidal and nanoparticle flux as reported in Canfield and Raiswell +const_fed = 70.0e-6; + +% GlobalNEWS2 data obtained from Emilio Mayorga +filename = '/archive/cas/COBALT_EVAL/River_Data/GlobalNEWS2/GlobalNEWS2_RH2000Dataset-version1.0.xls' +basin = readtable(filename,'Sheet',2); +hydrology = readtable(filename,'Sheet',3); +load = readtable(filename,'Sheet',4); + +% find all the river basins that empty into "land", e.g., lakes +ocean = basin.ocean; +land_index = zeros(size(ocean)); +for n = 1:size(ocean,1); + test = strcmp('Land',ocean(n)); + land_index(n) = single(test); +end + +river_names_all = basin.basinname; + +% basin area in +area = basin.A; +lon_news_all = basin.mouth_lon; +lat_news_all = basin.mouth_lat; +% Depending on the model output, may need to adjust longitudes +lon_news_all(lon_news_all < 0) = lon_news_all(lon_news_all < 0) + 360; + +% Loads in Mg yr-1 converted to moles per sec +DIN_load_all = load.Ld_DIN*1e6/14/86400/365; +DIP_load_all = load.Ld_DIP*1e6/31/86400/365; +DON_load_all = load.Ld_DON*1e6/14/86400/365; +DOP_load_all = load.Ld_DOP*1e6/31/86400/365; +Si_load_all = load.Ld_DSi*1e6/28.1/86400/365; +PN_load_all = (load.Ld_PN*1e6/14/86400/365); +PP_load_all = (load.Ld_PP*1e6/31/86400/365)*frac_PP; + +% actual and natural discharge (convert from km3/yr to m3/sec) +% Used the actual hydrology to calculate concentrations +Qact_all = hydrology.Qact*1e9/(86400*365); +Qnat_all = hydrology.Qnat*1e9/(86400*365); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Load in monthly climatology of river forcing from the regional grid. % +% File contains: % +% runoff: monthly average runoff in kg m-2 sec-1 % +% area: area of grid cell in m-2 % +% lon: longitude (0-360 degrees) % +% lat: latitude % +% % +% The file was calculated from daily output using the routine % +% "make_climatology.m". This routine and all associated files can be % +% found in the same directory as the data file % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% NEP +load glofas_hill_runoff_monthlyclim_NEP10k_05122023.mat; +% Arctic +%load glofas_hill_runoff_monthlyclim_arctic12k_05112023.mat; +lon_mod = lon; +lat_mod = lat; +area_mod = area; +% convert runoff from kg m-2 sec-1 to m3 sec-1 +Q_mod_monthly = zeros(size(runoff)); +for m = 1:12 + Q_mod_monthly(m,:,:) = squeeze(runoff(m,:,:)).*area_mod./1000; +end +Q_mod_ann = squeeze(mean(Q_mod_monthly,1)); +clear lon lat runoff area; + +%grid_file = '/archive/cas/Regional_MOM6/NWA12/nwa12_ocean_static.nc'; +%temp = ncread(grid_file,'deptho'); +%depth = permute(temp,[2,1]); clear temp; +%depth(isnan(depth)) = -1; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Filter for rivers in the region, set thresholds for minimum river size, % +% set parameters for plotting routines. % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% use grid to filter rivers outside domain +lat_mod_max = max(lat_mod(:)); +lat_mod_min = min(lat_mod(:)); +lon_mod_max = max(lon_mod(:)); +lon_mod_min = min(lon_mod(:)); + +in_region = find( (lon_news_all <= lon_mod_max) & (lon_news_all >= lon_mod_min) & ... + (lat_news_all <= lat_mod_max) & (lat_news_all >= lat_mod_min) & ... + (isfinite(Qact_all)) & (Qact_all > Q_min) ); + +% If you are using a high threshold, grab one smaller river to constrain +% Carribean Islands +%if Q_min > 100 +% for n = 1:size(lon_news_all,1); +% test = strcmp('GHAASBasin1808',river_names_all{n}); +% if test == 1; +% cuba_ind = n; +% end +% end +% in_region = [in_region; cuba_ind]; +%end + +num_rivers = size(in_region,1); + +% Establish vectors of flow and nutrient loads for the NWA +Qact = Qact_all(in_region); +lon_news = lon_news_all(in_region); +lat_news = lat_news_all(in_region); +DIN_load = DIN_load_all(in_region); +DON_load = DON_load_all(in_region); +PN_load = PN_load_all(in_region); +DIP_load = DIP_load_all(in_region); +DOP_load = DOP_load_all(in_region); +PP_load = PP_load_all(in_region); +Si_load = Si_load_all(in_region); +river_names = river_names_all(in_region); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Following inspection of initial mapping, add any manual edits here to % +% prevent anomalous extrapolations, etc. % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %Add in 2 "dummy" rivers to handle Cuba, giving it properties like + %Jamaica or Haiti rather than Florida. + + %GHAASBasin1808 is in Haiti. Fluxes are characterized by particularly + %high particulate phosphorus inputs. +%for n = 1:num_rivers; +% n +% test = strcmp('GHAASBasin1808',river_names{n}) +% if test == 1; +% cuba_ind = n; +% end + + % Move the Susquehanna a bit south so that it catches the Chesapeake + % and not the Delaware. +% if strcmp('Susquehanna',river_names{n}) +% lat_news(n) = 38.5; +% lon_news(n) = -76.67; +% end +%end + +% Two "rivers" with properties like Haiti used to specify Cuba +%Qact(num_rivers+1) = Qact(cuba_ind)/2; +%lon_news(num_rivers+1) = -81; +%lat_news(num_rivers+1) = 22.6; +%DIN_load(num_rivers+1) = DIN_load(cuba_ind)/2; +%DON_load(num_rivers+1) = DON_load(cuba_ind)/2; +%PN_load(num_rivers+1) = PN_load(cuba_ind)/2; +%DIP_load(num_rivers+1) = DIP_load(cuba_ind)/2; +%DOP_load(num_rivers+1) = DOP_load(cuba_ind)/2; +%PP_load(num_rivers+1) = PP_load(cuba_ind)/2; +%Si_load(num_rivers+1) = Si_load(cuba_ind)/2; +%river_names(num_rivers+1) = river_names(cuba_ind); + +%Qact(num_rivers+2) = Qact(cuba_ind)/2; +%lon_news(num_rivers+2) = -83.25; +%lat_news(num_rivers+2) = 22.6; +%DIN_load(num_rivers+2) = DIN_load(cuba_ind)/2; +%DON_load(num_rivers+2) = DON_load(cuba_ind)/2; +%PN_load(num_rivers+2) = PN_load(cuba_ind)/2; +%DIP_load(num_rivers+2) = DIP_load(cuba_ind)/2; +%DOP_load(num_rivers+2) = DOP_load(cuba_ind)/2; +%PP_load(num_rivers+2) = PP_load(cuba_ind)/2; +%Si_load(num_rivers+2) = Si_load(cuba_ind)/2; +%river_names(num_rivers+2) = river_names(cuba_ind); + +%num_rivers = num_rivers + 1; + +% Adjust location of cfilename = '/archive/cas/COBALT_EVAL/River_Data/GlobalNEWS2/GlobalNEWS2_RH2000Dataset-version1.0.xls' +basin = readtable(filename,'Sheet',2); +hydrology = readtable(filename,'Sheet',3); +load = readtable(filename,'Sheet',4); +% Cuba; This has little effect on patterns in Florida. +%for n = 1:num_rivers; +% n +% test = strcmp('GHAASBasin448',river_names{n}) +% if test == 1; +% fla_ind = n; +% end +%end + +%lon_news(fla_ind) = -80.5; +%lat_news(fla_ind) = 26.6; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% END MANUAL EDITS % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Assigning outflow points to rivers. % +% 1. Assignment starts with the rivers with the smallest flow and works % +% to the largest, w/larger river characteristics taking precedence to % +% ensure the most significant rivers are well represented. % +% 2. The algorithm keeps choosing the closest points to each river mouth % +% until the assigned flow is as close as possible to that observed % +% 3. Once the outflow points are assigned using the mean flow values, % +% monthly concentrations are assigned to those points. % +% 4. A simple "nearest neighbor" algorithm is used to fill in the gaps % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Sort rivers by discharge +[Qact_sort,sort_ind] = sort(Qact,'ascend'); +lon_news_sort = lon_news(sort_ind); +lat_news_sort = lat_news(sort_ind); +DIN_load_sort = DIN_load(sort_ind); +DON_load_sort = DON_load(sort_ind); +PN_load_sort = PN_load(sort_ind); +DIP_load_sort = DIP_load(sort_ind); +DOP_load_sort = DOP_load(sort_ind); +PP_load_sort = PP_load(sort_ind); +Si_load_sort = Si_load(sort_ind); +river_names_sort = river_names(sort_ind); + +% Total N and P load diagnostics +N_load_sort = DIN_load_sort + DON_load_sort + PN_load_sort; +P_load_sort = DIP_load_sort + DOP_load_sort + PP_load_sort; + +% Calculate concentrations +% Loads are in moles N sec-1, Q in m3 s-1; conc in moles N m-3 +DIN_conc_sort = DIN_load_sort./Qact_sort; +DON_conc_sort = DON_load_sort./Qact_sort; +DIP_conc_sort = DIP_load_sort./Qact_sort; +DOP_conc_sort = DOP_load_sort./Qact_sort; +PN_conc_sort = PN_load_sort./Qact_sort; +PP_conc_sort = PP_load_sort./Qact_sort; +Si_conc_sort = Si_load_sort./Qact_sort; + +% initialize vectors to hold nutrient concentrations at eac runoff +% point. +aa = find(Q_mod_ann > 0); +Q_mod_vec = Q_mod_ann(aa); +din_conc_vec = zeros(size(Q_mod_vec)); +don_conc_vec = zeros(size(Q_mod_vec)); +pn_conc_vec = zeros(size(Q_mod_vec)); +dip_conc_vec = zeros(size(Q_mod_vec)); +dop_conc_vec = zeros(size(Q_mod_vec)); +pp_conc_vec = zeros(size(Q_mod_vec)); +si_conc_vec = zeros(size(Q_mod_vec)); + +lon_mod_runoff_vec = double(lon_mod(aa)); +lat_mod_runoff_vec = double(lat_mod(aa)); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Loop identifies points assigned to each river % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +for k=1:num_rivers + k + dist = pdist2([lon_news_sort(k) lat_news_sort(k)], ... + [lon_mod_runoff_vec lat_mod_runoff_vec]); + [dist_sort, dist_sort_ind] = sort(dist,'ascend'); + + if dist_sort(1) < min_dist; % filters out rivers lying outside the domain + Q_sum1 = 0; + Q_sum2 = 0; + n = 0; + while ((Q_sum2 < Qact_sort(k)) & (dist_sort(n+1) < max_dist)) + Q_sum1 = Q_sum2; + n = n+1; + Q_sum2 = Q_sum1 + Q_mod_vec(dist_sort_ind(n)); + end + % I generally find that the search algorithm works best when you keep + % grabbing points until the total flow captured exceeds the flow in + % the river. If you uncomment the the first part of the "if" + % statement here will pick the last below if it is closer than the + % first above. However, I've found that this option sometimes fails + % to map important rivers. + %if abs(Q_sum1 - Qact_sort(k)) < abs(Q_sum2 - Qact_sort(k)) + % nrp = n-1; % number of runoff points + % [Q_sum1 Qact_sort(k)] % a quick check for comparable flow + % din_conc_vec(dist_sort_ind(1:nrp)) = DIN_conc_sort(k); + % don_conc_vec(dist_sort_ind(1:nrp)) = DON_conc_sort(k); + % dip_conc_vec(dist_sort_ind(1:nrp)) = DIP_conc_sort(k); + % dop_conc_vec(dist_sort_ind(1:nrp)) = DOP_conc_sort(k); + % pn_conc_vec(dist_sort_ind(1:nrp)) = PN_conc_sort(k); + % pp_conc_vec(dist_sort_ind(1:nrp)) = PP_conc_sort(k); + % si_conc_vec(dist_sort_ind(1:nrp)) = Si_conc_sort(k); + %else + nrp = n; % number of runoff points + [Q_sum2 Qact_sort(k)] % a quick check for comparable flow + din_conc_vec(dist_sort_ind(1:nrp)) = DIN_conc_sort(k); + don_conc_vec(dist_sort_ind(1:nrp)) = DON_conc_sort(k); + dip_conc_vec(dist_sort_ind(1:nrp)) = DIP_conc_sort(k); + dop_conc_vec(dist_sort_ind(1:nrp)) = DOP_conc_sort(k); + pn_conc_vec(dist_sort_ind(1:nrp)) = PN_conc_sort(k); + pp_conc_vec(dist_sort_ind(1:nrp)) = PP_conc_sort(k); + si_conc_vec(dist_sort_ind(1:nrp)) = Si_conc_sort(k); + %end + + if inspect_map == 'y' + figure(1) + clf + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(Q_mod_vec),3,log10(Q_mod_vec)); + hold on + scatter3(lon_mod_runoff_vec(dist_sort_ind(1:nrp)),lat_mod_runoff_vec(dist_sort_ind(1:nrp)), ... + log10(Q_mod_vec(dist_sort_ind(1:nrp))),40, ... + log10(Q_mod_vec(dist_sort_ind(1:nrp))),'filled'); + view(2); + plot3(lon_news_sort(k),lat_news_sort(k),1e5,'k.','MarkerSize',20); + %contour(lon_mod,lat_mod,depth,[0 0],'k-'); + axis([lon_news_sort(k)-plot_width lon_news_sort(k)+plot_width ... + lat_news_sort(k)-plot_width lat_news_sort(k)+plot_width]); + caxis([-4 3]); + titl = ['river number: ',num2str(k),' name: ',river_names_sort{k}]; + title(titl); + colorbar; + + % provide a few diagnostics to ensure the calculation was done + % correctly (remove semicolon to inspect as they are mapped in + % the matlab output line. Feel free to add more here. + N_conc = DIN_conc_sort(k) + DON_conc_sort(k) + PN_conc_sort(k); + P_conc = DIP_conc_sort(k) + DOP_conc_sort(k) + PP_conc_sort(k); + Si_conc = Si_conc_sort(k); + river_names_sort(k) + [lon_news_sort(k) lat_news_sort(k)] + ind = dist_sort_ind(1:nrp); + 'lon lat' + [lon_news_sort(k) lat_news_sort(k)] + 'total flow in m3 sec' + [Qact_sort(k) sum(Q_mod_vec(ind))] + 'Nitrogen and phosphorus in Gg per year'; + [N_load_sort(k) sum(Q_mod_vec(ind))*N_conc]*14*86400*365/1e9; + [P_load_sort(k) sum(Q_mod_vec(ind))*P_conc]*31*86400*365/1e9; + 'N, P conc (mmoles m-3), DI, DO, P' + [DIN_conc_sort(k) DON_conc_sort(k) PN_conc_sort(k)]*1e3 + [DIP_conc_sort(k) DOP_conc_sort(k) PP_conc_sort(k)]*1e3 + 'Total N, Total P, Total N: Total P'; + [N_conc*1e3 P_conc*1e3 N_conc/P_conc] + 'DO:DI and P:DI ratios'; + [DON_conc_sort(k)/DIN_conc_sort(k) PN_conc_sort(k)/DIN_conc_sort(k)]; + [DOP_conc_sort(k)/DIP_conc_sort(k) PP_conc_sort(k)/DIP_conc_sort(k)]; + 'silica concentration (mmoles m-3)'; + [Si_conc_sort(k)]*1e3; + pause + end + + else + % This is for rivers that were captured by the coarse initial filter + % to determine if they were in the domain; but are actually outside + % the domain. Calibration suggested that a threshold of 0.75 degrees + % from the specified river mouth reliably identified these cases. + + if inspect_map == 'y' + figure(1) + clf + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(Q_mod_vec),3,log10(Q_mod_vec)); + hold on + view(2); + plot3(lon_news_sort(k),lat_news_sort(k),1e5,'k.','MarkerSize',20); + axis([lon_news_sort(k)-15 lon_news_sort(k)+15 ... + lat_news_sort(k)-15 lat_news_sort(k)+15]); + caxis([-4 3]); + titl = ['OUTSIDE: river number: ',num2str(k),' name: ',river_names_sort{k}]; + title(titl); + colorbar; + pause + end + + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% nearest neighbor search to fill in any 0 values left for each input field +% after the river mapping is done. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +aa = find(din_conc_vec == 0); +bb = find(din_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),din_conc_vec(bb), ... + 'nearest','nearest'); +din_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(don_conc_vec == 0); +bb = find(don_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),don_conc_vec(bb), ... + 'nearest','nearest'); +don_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(pn_conc_vec == 0); +bb = find(pn_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),pn_conc_vec(bb), ... + 'nearest','nearest'); +pn_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(dip_conc_vec == 0); +bb = find(dip_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),dip_conc_vec(bb), ... + 'nearest','nearest'); +dip_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(dop_conc_vec == 0); +bb = find(dop_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),dop_conc_vec(bb), ... + 'nearest','nearest'); +dop_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(pp_conc_vec == 0); +bb = find(pp_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),pp_conc_vec(bb), ... + 'nearest','nearest'); +pp_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +aa = find(si_conc_vec == 0); +bb = find(si_conc_vec > 0); +F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb),si_conc_vec(bb), ... + 'nearest','nearest'); +si_conc_vec(aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + +totn_conc_vec = din_conc_vec + don_conc_vec + pn_conc_vec; +totp_conc_vec = dip_conc_vec + dop_conc_vec + pp_conc_vec; + +% calculate ratios of dissolved and particulate to inorganic +din_flux_vec = din_conc_vec.*Q_mod_vec; +dip_flux_vec = dip_conc_vec.*Q_mod_vec; +don_flux_vec = don_conc_vec.*Q_mod_vec; +dop_flux_vec = dop_conc_vec.*Q_mod_vec; +pn_flux_vec = pn_conc_vec.*Q_mod_vec; +pp_flux_vec = pp_conc_vec.*Q_mod_vec; + +don_ratio = sum(don_flux_vec)/sum(din_flux_vec) +dop_ratio = sum(dop_flux_vec)/sum(dip_flux_vec) +pn_ratio = sum(pn_flux_vec)/sum(din_flux_vec) +pp_ratio = sum(pp_flux_vec)/sum(dip_flux_vec) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Produce plots to evaluate the mapping % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Use total fluxes to scale dots, m3 sec-1 * moles m-3 = moles sec-1 +totn_flux_vec = totn_conc_vec.*Q_mod_vec; +totp_flux_vec = totp_conc_vec.*Q_mod_vec; +totsi_flux_vec = si_conc_vec.*Q_mod_vec; + +% scale marker size with the total nitrogen flux +ms_vec = zeros(size(Q_mod_vec)); +ms_vec(log10(Q_mod_vec) < 0) = 1; +ms_vec((log10(Q_mod_vec) > 0) & (log10(Q_mod_vec) < 1)) = 1; +ms_vec((log10(Q_mod_vec) > 1) & (log10(Q_mod_vec) < 2)) = 1; +ms_vec((log10(Q_mod_vec) > 2) & (log10(Q_mod_vec) < 3)) = 1; +ms_vec(log10(Q_mod_vec) > 3) = 100; + +figure(1) +clf + +subplot(1,3,1); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,totn_conc_vec*1e3,ms_vec, ... + totn_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([30 150]); +colorbar +title('total nitrogen concentration, mmoles m-3'); + +subplot(1,3,2); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,totp_conc_vec*1e3,ms_vec, ... + totp_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 10]); +colorbar +title('total phosphorus concentration, mmoles m-3'); + +subplot(1,3,3) +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,totn_conc_vec./totp_conc_vec,ms_vec, ... + totn_conc_vec./totp_conc_vec,'filled'); +hold on +view(2); +caxis([10 100]); +colorbar +title('N:P ratio'); + +figure(2) +clf + +subplot(1,3,1); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,din_conc_vec*1e3,ms_vec, ... + din_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 100]); +colorbar +title('din concentration, mmoles m-3'); + +subplot(1,3,2); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,don_conc_vec*1e3,ms_vec, ... + don_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 100]); +colorbar +title('don concentration, mmoles m-3'); + +subplot(1,3,3) +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pn_conc_vec*1e3,ms_vec, ... + pn_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 100]); +colorbar +title('pn concentration, mmoles m-3'); + +figure(3) +clf + +subplot(1,3,1); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dip_conc_vec*1e3,ms_vec, ... + dip_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 5]); +colorbar +title('dip concentration, mmoles m-3'); + +subplot(1,3,2); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dop_conc_vec*1e3,ms_vec, ... + dop_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 3]); +colorbar +title('dop concentration, mmoles m-3'); + +subplot(1,3,3) +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pp_conc_vec*1e3,ms_vec, ... + pp_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([0 3]); +colorbar +title('pp concentration, mmoles m-3'); + +figure(10) +clf + +subplot(2,2,1); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,don_conc_vec./din_conc_vec,ms_vec, ... + don_conc_vec./din_conc_vec,'filled'); +hold on +view(2); +caxis([0 2]); +colorbar +title('don:din'); + +subplot(2,2,2); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pn_conc_vec./din_conc_vec,ms_vec, ... + pn_conc_vec./din_conc_vec,'filled'); +hold on +view(2); +caxis([0 2]); +colorbar +title('pn:din'); + +subplot(2,2,3) +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dop_conc_vec./dip_conc_vec,ms_vec, ... + dop_conc_vec./dip_conc_vec,'filled'); +hold on +view(2); +caxis([0 2]); +colorbar +title('dop:dip'); + +subplot(2,2,4) +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pp_conc_vec./dip_conc_vec,ms_vec, ... + pp_conc_vec./dip_conc_vec,'filled'); +hold on +view(2); +caxis([0 2]); +colorbar +title('pp:dip'); + +figure(5) +clf +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,si_conc_vec*1e3,ms_vec, ... + si_conc_vec*1e3,'filled'); +hold on +view(2); +caxis([30 300]); +colorbar +title('si concentration, mmoles m-3'); + +figure(6) +clf + +subplot(1,2,1); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(totn_flux_vec),ms_vec, ... + log10(totn_flux_vec),'filled'); +hold on +view(2); +caxis([-1 2]); +colorbar +title('N flux, log10(moles sec-1)'); + +subplot(1,2,2); +scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(totp_flux_vec),ms_vec, ... + log10(totp_flux_vec),'filled'); +hold on +view(2); +caxis([-3 1]); +colorbar +title('P flux, log10(moles sec-1)'); + +% output total fluxes in gigagrams of N +total_annual_n = sum(totn_flux_vec)*86400*365*14/1e12 +total_annual_p = sum(totp_flux_vec)*86400*365*31/1e12 +total_annual_si = sum(totsi_flux_vec)*86400*365*28.1/1e12 + +'press any key to continue' +pause + +% Initialize 2D concentration arrays; these are the ones read into MOM6 to +% specify the nutrient concentrations of river inputs. +din_conc = zeros(size(lon_mod)); +don_conc = zeros(size(lon_mod)); +dip_conc = zeros(size(lon_mod)); +dop_conc = zeros(size(lon_mod)); +pn_conc = zeros(size(lon_mod)); +pp_conc = zeros(size(lon_mod)); +si_conc = zeros(size(lon_mod)); + +% Map concentration vectors onto 2D arrays. +aa = find(Q_mod_ann > 0); +din_conc(aa) = din_conc_vec; +don_conc(aa) = don_conc_vec; +pn_conc(aa) = pn_conc_vec; +dip_conc(aa) = dip_conc_vec; +dop_conc(aa) = dop_conc_vec; +pp_conc(aa) = pp_conc_vec; +si_conc(aa) = si_conc_vec; + +NO3_CONC = din_conc; +LDON_CONC = frac_ldon*don_conc; +SLDON_CONC = frac_sldon*don_conc; +SRDON_CONC = frac_srdon*don_conc; +PO4_CONC = dip_conc; +LDOP_CONC = frac_ldop*dop_conc; +SLDOP_CONC = frac_sldop*dop_conc; +SRDOP_CONC = frac_srdop*dop_conc; +NDET_CONC = pn_conc; +PDET_CONC = pp_conc; % The bioavailability of particulate P has already + % been accounted for. + +SI_CONC = si_conc; + +% Add iron concentrations - initialize with nitrate and then overwrite +FED_CONC = NO3_CONC; +FEDET_CONC = NO3_CONC; +% 40 nM dissolved iron concentration from De Baar and De Jong + 30nM +% Colloidal and nanoparticle flux as reported in Canfield and Raiswell +FED_CONC(FED_CONC > 0) = const_fed; +FEDET_CONC(FEDET_CONC > 0) = 0.0; + +% series of quick figures to check the 2D nutrient input files. +ms = 8; +figure(7); +clf +subplot(3,2,1); +title('log10(NO3 CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(NO3_CONC(:)),ms,log10(NO3_CONC(:)),'filled'); +caxis([-4 -1]); colorbar; + +subplot(3,2,2); +title('log10(LDON CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(LDON_CONC(:)),ms,log10(LDON_CONC(:)),'filled'); +caxis([-4 -1]); colorbar; + +subplot(3,2,3); +title('log10(SLDON CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(SLDON_CONC(:)),ms,log10(SLDON_CONC(:)),'filled'); +caxis([-4 -1]); colorbar; + +subplot(3,2,4); +title('log10(SRDON CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(SRDON_CONC(:)),ms,log10(SRDON_CONC(:)),'filled'); +caxis([-4 -1]); colorbar; + +subplot(3,2,5); +title('log10(NDET CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(NDET_CONC(:)),ms,log10(NDET_CONC(:)),'filled'); +caxis([-4 -1]); colorbar; + +figure(8); +clf +subplot(3,2,1); +title('log10(PO4 CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(PO4_CONC(:)),ms,log10(PO4_CONC(:)),'filled'); +caxis([-4 -2]); colorbar; + +subplot(3,2,2); +title('log10(LDOP CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(LDOP_CONC(:)),ms,log10(LDOP_CONC(:)),'filled'); +caxis([-4 -2]); colorbar; + +subplot(3,2,3); +title('log10(SLDOP CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(SLDOP_CONC(:)),ms,log10(SLDOP_CONC(:)),'filled'); +caxis([-4 -2]); colorbar; + +subplot(3,2,4); +title('log10(SRDOP CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(SRDOP_CONC(:)),ms,log10(SRDOP_CONC(:)),'filled'); +caxis([-4 -2]); colorbar; + +subplot(3,2,5); +title('log10(PDET CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(PDET_CONC(:)),ms,log10(PDET_CONC(:)),'filled'); +caxis([-4 -2]); colorbar; + +figure(9) +clf +clf +subplot(3,2,1); +title('log10(FED CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(FED_CONC(:)),ms,log10(FED_CONC(:)),'filled'); +caxis([-5 -3]); colorbar; + +subplot(3,2,2); +title('log10(FEDET CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(FEDET_CONC(:)),ms,log10(FEDET_CONC(:)),'filled'); +caxis([-5 -3]); colorbar; + +subplot(3,2,3); +title('log10(SI CONC)'); hold on; +scatter3(lon_mod(:),lat_mod(:),log10(SI_CONC(:)),ms,log10(SI_CONC(:)),'filled'); +caxis([-3 -0]); colorbar; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Save Files % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% option to save matlab file +% save River_DIC_ALK_USGS_NWA ALK_CONC DIC_CONC + +% Construct netcdf file following format used by nutrient input files to +% MOM6 +time = 0; +nlat = size(lat_mod,1); +nlon = size(lat_mod,2); + +ncid = netcdf.create(nc_file_name,'CLOBBER'); +dimid0 = netcdf.defDim(ncid,'time',netcdf.getConstant('NC_UNLIMITED')); + +dimid1 = netcdf.defDim(ncid,'y',nlat); +dimid2 = netcdf.defDim(ncid,'x',nlon); + +%attributes inherited from the old river nutrient forcing file. The +%calendar needs to be specified even though the nutrient values are +%constant. Not sure how much of the rest is needed. +varid0 = netcdf.defVar(ncid,'time','double',dimid0); +netcdf.putAtt(ncid,varid0,'calendar','NOLEAP'); +netcdf.putAtt(ncid,varid0,'calendar_type','NOLEAP'); +netcdf.putAtt(ncid,varid0,'modulo','T'); +netcdf.putAtt(ncid,varid0,'units','days since 1900-1-1 0:00:00'); +netcdf.putAtt(ncid,varid0,'time_origin','01-JAN-1990 00:00:00'); + +varid1 = netcdf.defVar(ncid,'y','int',dimid1); +netcdf.putAtt(ncid,varid1,'cartesian_axis','Y'); +varid2 = netcdf.defVar(ncid,'x','int',dimid2); +netcdf.putAtt(ncid,varid2,'cartesian_axis','X'); +varid3 = netcdf.defVar(ncid,'lat','double',[dimid2 dimid1]); +netcdf.putAtt(ncid,varid3,'units','degrees north'); +varid4 = netcdf.defVar(ncid,'lon','double',[dimid2 dimid1]); +netcdf.putAtt(ncid,varid4,'units','degrees east'); +varid5 = netcdf.defVar(ncid,'NO3_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid5,'units','mol m-3'); +netcdf.putAtt(ncid,varid5,'long_name','DIN_CONC'); +varid6 = netcdf.defVar(ncid,'LDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid6,'units','mol m-3'); +netcdf.putAtt(ncid,varid6,'long_name','0.3*DON_CONC'); +varid7 = netcdf.defVar(ncid,'SLDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid7,'units','mol m-3'); +netcdf.putAtt(ncid,varid7,'long_name','0.35*DON_CONC'); +varid8 = netcdf.defVar(ncid,'SRDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid8,'units','mol m-3'); +netcdf.putAtt(ncid,varid8,'long_name','0.35*DON_CONC'); +varid9 = netcdf.defVar(ncid,'NDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid9,'units','mol m-3'); +netcdf.putAtt(ncid,varid9,'long_name','1.0*PN_CONC'); +varid10 = netcdf.defVar(ncid,'PO4_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid10,'units','mol m-3'); +netcdf.putAtt(ncid,varid10,'long_name','PO4_CONC'); +varid11 = netcdf.defVar(ncid,'LDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid11,'units','mol m-3'); +netcdf.putAtt(ncid,varid11,'long_name','0.3*DOP_CONC'); +varid12 = netcdf.defVar(ncid,'SLDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid12,'units','mol m-3'); +netcdf.putAtt(ncid,varid12,'long_name','0.35*DOP_CONC'); +varid13 = netcdf.defVar(ncid,'SRDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid13,'units','mol m-3'); +netcdf.putAtt(ncid,varid13,'long_name','0.35*DOP_CONC'); +varid14 = netcdf.defVar(ncid,'PDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid14,'units','mol m-3'); +netcdf.putAtt(ncid,varid14,'long_name','0.3*PP_CONC'); +varid15 = netcdf.defVar(ncid,'FED_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid15,'units','mol m-3'); +netcdf.putAtt(ncid,varid15,'long_name','FED_CONC'); +varid16 = netcdf.defVar(ncid,'FEDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid16,'units','mol m-3'); +netcdf.putAtt(ncid,varid16,'long_name','FEDET_CONC'); +varid17 = netcdf.defVar(ncid,'SI_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid17,'units','mol m-3'); +netcdf.putAtt(ncid,varid17,'long_name','SI_CONC'); +netcdf.close(ncid) + +ncid = netcdf.open(nc_file_name,'NC_WRITE'); +netcdf.putVar(ncid,varid0,0,1,time); +% nutrient input files appear seem to need dummy axes to be read in +% properly, but eventually do a grid by grid mapping that doesn't require +% these. +netcdf.putVar(ncid,varid1,1:nlat); +netcdf.putVar(ncid,varid2,1:nlon); +netcdf.putVar(ncid,varid3,permute(lon_mod,[2,1])); +netcdf.putVar(ncid,varid4,permute(lat_mod,[2,1])); +netcdf.putVar(ncid,varid5,permute(NO3_CONC,[3,2,1])); +netcdf.putVar(ncid,varid6,permute(LDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid7,permute(SLDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid8,permute(SRDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid9,permute(NDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid10,permute(PO4_CONC,[3,2,1])); +netcdf.putVar(ncid,varid11,permute(LDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid12,permute(SLDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid13,permute(SRDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid14,permute(PDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid15,permute(FED_CONC,[3,2,1])); +netcdf.putVar(ncid,varid16,permute(FEDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid17,permute(SI_CONC,[3,2,1])); +netcdf.close(ncid) diff --git a/tools/rivers/bgc/NEP/mapriv_combined_NEP10k.m b/tools/rivers/bgc/NEP/mapriv_combined_NEP10k.m new file mode 100644 index 000000000..9583ca265 --- /dev/null +++ b/tools/rivers/bgc/NEP/mapriv_combined_NEP10k.m @@ -0,0 +1,1342 @@ +% Routine to map USGS nutrient data onto the MOM6 Northwest Atlantic (NWA) +% grid. Run on matlab97 or above. + +clear all; +addpath /home/cas/matlab +nc64startup + +% name of netcdf file to be created +nc_file_name = 'RiverNutrients_Combined_Q100_NEP10k.nc'; + +% GLOBAL NEWS based map for filling in gaps +NEWS_file = 'RiverNutrients_GlobalNEWS2_plusFe_Q100_NEP10k.nc'; + +% load in monthly world ocean T, S climatology for saturated oxygen calculation +temp = ncread('Data/woa_sst_climo.nc','t_an'); +woa_temp = permute(temp,[3 2 1]); + +% Parameters for the assignment algorithm. +Q_min = 100; % minimum flow in m3 sec +plot_width = 15; % width of window (in degrees) for inspecting locations + % of rivers and outflow points that have been assigned to + % them. +min_dist = 1.5; % minimum distance (degrees) of the closest outflow point + % for the river to be considered in the domain (useful + % for preventing the algorithm from trying to map rivers + % flowing to different ocean basins. +max_dist = 2.0; % maximum distance (degrees) away that the algorithm + % looks for points for rivers that are in the domain +nutrient_option = 2; % option for deriving dissolved organic nutrients in RC4US +inspect_map = 'y'; % flag enabling you to pause and inspect each river + % mapping as it is being done. + +min_lon_ref = 0; % set to either 0 if model grid contains no negative + % values; set to -180 if model is on a -180-180 grid. + +% set the bio-availability of phosphorus and the fractionation of dissolved +% organic; PP is set to 30% based on Froelich; Partitioning of detritus +% between +frac_PP = 0.3; +frac_ldon = 0.3; +frac_sldon = 0.35; +frac_srdon = 0.35; +frac_ldop = 0.3; +frac_sldop = 0.35; +frac_srdop = 0.35; +% 40 nM dissolved iron concentration from De Baar and De Jong + 30nM +% Colloidal and nanoparticle flux as reported in Canfield and Raiswell +const_fed = 70.0e-6; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% USGS data compiled by Fabian Gomez % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +filename_chem = 'Data/RC4USCoast/mclim_19902022_chem.nc'; +alk_monthly_RC4US = ncread(filename_chem,'alk'); alk_monthly_RC4US = permute(alk_monthly_RC4US,[2 1]); +dic_monthly_RC4US = ncread(filename_chem,'dic'); dic_monthly_RC4US = permute(dic_monthly_RC4US,[2 1]); +no3_monthly_RC4US = ncread(filename_chem,'no3'); no3_monthly_RC4US = permute(no3_monthly_RC4US,[2 1]); +nh4_monthly_RC4US = ncread(filename_chem,'nh4'); nh4_monthly_RC4US = permute(nh4_monthly_RC4US,[2 1]); +din_monthly_RC4US = no3_monthly_RC4US + nh4_monthly_RC4US; +dip_monthly_RC4US = ncread(filename_chem,'po4'); dip_monthly_RC4US = permute(dip_monthly_RC4US,[2 1]); +si_monthly_RC4US = ncread(filename_chem,'sio2'); si_monthly_RC4US = permute(si_monthly_RC4US,[2 1]); +% The RC4US database seems to be in mmoles O m-3 rather than mmoles O2 m-3, +% divide by 2.0 for consistency with other O2 data sources +o2_monthly_RC4US = ncread(filename_chem,'do'); o2_monthly_RC4US = permute(o2_monthly_RC4US,[2 1])/2.0; +don_monthly_RC4US = ncread(filename_chem,'don'); don_monthly_RC4US = permute(don_monthly_RC4US,[2 1]); +temp_monthly_RC4US = ncread(filename_chem,'temp'); temp_monthly_RC4US = permute(temp_monthly_RC4US,[2 1]); +if nutrient_option == 1 + % This option will eventually set all river values for pn, dop and pp using GlobalNEWS + pn_monthly_RC4US = no3_monthly_RC4US*NaN; + dop_monthly_RC4US = no3_monthly_RC4US*NaN; + pp_monthly_RC4US = no3_monthly_RC4US*NaN; +elseif nutrient_option == 2 + % This option will use differences between total filtered and unfiltered + % and other properties to derive pn, dop and pp. This unfortunately + % generates negative values in some cases. + tnf_monthly_RC4US = ncread(filename_chem,'tnf'); tnf_monthly_RC4US = permute(tnf_monthly_RC4US,[2 1]); + don2_monthly_RC4US = tnf_monthly_RC4US - din_monthly_RC4US; + tnu_monthly_RC4US = ncread(filename_chem,'tnu'); tnu_monthly_RC4US = permute(tnu_monthly_RC4US,[2 1]); + pn_monthly_RC4US = tnu_monthly_RC4US - tnf_monthly_RC4US; + pn_monthly_RC4US(pn_monthly_RC4US < 0) = NaN; + tpf_monthly_RC4US = ncread(filename_chem,'tpf'); tpf_monthly_RC4US = permute(tpf_monthly_RC4US,[2 1]); + dop_monthly_RC4US = tpf_monthly_RC4US - dip_monthly_RC4US; + dop_monthly_RC4US(dop_monthly_RC4US < 0) = NaN; + tpu_monthly_RC4US = ncread(filename_chem,'tpu'); tpu_monthly_RC4US = permute(tpu_monthly_RC4US,[2 1]); + pp_monthly_RC4US = (tpu_monthly_RC4US - tpf_monthly_RC4US)*frac_PP; + pp_monthly_RC4US(pp_monthly_RC4US < 0) = NaN; +end +dfe_monthly_RC4US = no3_monthly_RC4US*NaN; +pfe_monthly_RC4US = no3_monthly_RC4US*NaN; + +filename_discharge = 'Data/RC4USCoast/mclim_19902022_disc.nc'; +Q_monthly_RC4US = ncread(filename_discharge,'disc'); Q_monthly_RC4US = permute(Q_monthly_RC4US,[2 1]); % m-3 sec-1 +station_names_RC4US = h5read(filename_discharge,'/river_name'); +lon_stations_RC4US = ncread(filename_discharge,'mouth_lon'); +lat_stations_RC4US = ncread(filename_discharge,'mouth_lat'); + +Q_ann_RC4US = mean(Q_monthly_RC4US,1,'native','omitnan')'; +dic_ann_RC4US = mean(dic_monthly_RC4US,1,'native','omitnan')'; +alk_ann_RC4US = mean(alk_monthly_RC4US,1,'native','omitnan')'; +no3_ann_RC4US = mean(no3_monthly_RC4US,1,'native','omitnan')'; +nh4_ann_RC4US = mean(nh4_monthly_RC4US,1,'native','omitnan')'; +o2_ann_RC4US = mean(o2_monthly_RC4US,1,'native','omitnan')'; +dip_ann_RC4US = mean(dip_monthly_RC4US,1,'native','omitnan')'; +si_ann_RC4US = mean(si_monthly_RC4US,1,'native','omitnan')'; +din_ann_RC4US = no3_ann_RC4US + nh4_ann_RC4US; +don_ann_RC4US = mean(don_monthly_RC4US,1,'native','omitnan')'; +if nutrient_option == 1 + don_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + pn_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + dop_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + pp_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; +elseif nutrient_option == 2 + don2_ann_RC4US = mean(don2_monthly_RC4US,1,'native','omitnan')'; + pn_ann_RC4US = mean(pn_monthly_RC4US,1,'native','omitnan')'; + dop_ann_RC4US = mean(dop_monthly_RC4US,1,'native','omitnan')'; + pp_ann_RC4US = mean(pp_monthly_RC4US,1,'native','omitnan')'; +end +dfe_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; +pfe_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + +for n = 1:size(lon_stations_RC4US,1) + + % Make any adjustments to the river locations here, e.g: + % + % Move the Susquehanna a bit south so that it catches the Chesapeake + % and not the Delaware. + %if strcmp('Susquehanna',station_names_RC4US{n}) + % lat_stations_RC4US(n) = 38.5; + % lon_stations_RC4US(n) = -77.5; + % %pause + %end + +end + +Q_ann_RC4US = mean(Q_monthly_RC4US,1,'native','omitnan')'; +dic_ann_RC4US = mean(dic_monthly_RC4US,1,'native','omitnan')'; +alk_ann_RC4US = mean(alk_monthly_RC4US,1,'native','omitnan')'; +no3_ann_RC4US = mean(no3_monthly_RC4US,1,'native','omitnan')'; +nh4_ann_RC4US = mean(nh4_monthly_RC4US,1,'native','omitnan')'; +o2_ann_RC4US = mean(o2_monthly_RC4US,1,'native','omitnan')'; +dip_ann_RC4US = mean(dip_monthly_RC4US,1,'native','omitnan')'; +si_ann_RC4US = mean(si_monthly_RC4US,1,'native','omitnan')'; +din_ann_RC4US = no3_ann_RC4US + nh4_ann_RC4US; +don_ann_RC4US = mean(don_monthly_RC4US,1,'native','omitnan')'; +if nutrient_option == 1 + don_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + pn_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + dop_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + pp_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; +elseif nutrient_option == 2 + don2_ann_RC4US = mean(don2_monthly_RC4US,1,'native','omitnan')'; + pn_ann_RC4US = mean(pn_monthly_RC4US,1,'native','omitnan')'; + dop_ann_RC4US = mean(dop_monthly_RC4US,1,'native','omitnan')'; + pp_ann_RC4US = mean(pp_monthly_RC4US,1,'native','omitnan')'; +end +dfe_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; +pfe_ann_RC4US = ones(size(lon_stations_RC4US))*NaN; + +% Created by Process_GLORICH_NEP.m, includes following variables: +%save NEP_GLORICH_data lon_stations_glorich lat_stations_glorich station_names_glorich ... +% Q_ann_glorich dic_ann_glorich alk_ann_glorich no3_ann_glorich nh4_ann_glorich ... +% din_ann_glorich pn_ann_glorich don_ann_glorich dip_ann_glorich dop_ann_glorich ... +% pp_ann_glorich si_ann_glorich o2_ann_glorich dfe_ann_glorich pfe_ann_glorich ... +% dic_monthly_glorich alk_monthly_glorich no3_monthly_glorich nh4_monthly_glorich ... +% din_monthly_glorich don_monthly_glorich pn_monthly_glorich dip_monthly_glorich ... +% dop_monthly_glorich pp_monthly_glorich dfe_monthly_glorich pfe_monthly_glorich ... +% si_monthly_glorich o2_monthly_glorich +load Data/GLORICH/NEP_GLORICH_data.mat; + +% Created by Process_ARCTICGRO.m, includes following variables: +%save arcticgro_data lon_stations_arcticgro lat_stations_arcticgro station_names_arcticgro ... +% Q_ann_arcticgro dic_ann_arcticgro alk_ann_arcticgro no3_ann_arcticgro nh4_ann_arcticgro ... +% din_ann_arcticgro pn_ann_arcticgro don_ann_arcticgro dip_ann_arcticgro dop_ann_arcticgro ... +% pp_ann_arcticgro si_ann_arcticgro o2_ann_arcticgro dfe_ann_arcticgro pfe_ann_arcticgro ... +% dic_monthly_arcticgro alk_monthly_arcticgro no3_monthly_arcticgro nh4_monthly_arcticgro ... +% din_monthly_arcticgro don_monthly_arcticgro pn_monthly_arcticgro dip_monthly_arcticgro ... +% dop_monthly_arcticgro pp_monthly_arcticgro dfe_monthly_arcticgro pfe_monthly_arcticgro ... +% si_monthly_arcticgro o2_monthly_arcticgro +load Data/ArcticGro/arcticgro_data.mat; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Combine all annual and monthly station data % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +station_names_all = [station_names_RC4US; string(station_names_glorich)'; string(station_names_arcticgro)']; +lon_stations_all = [lon_stations_RC4US; lon_stations_glorich'; lon_stations_arcticgro']; +if min_lon_ref == 0 + lon_stations_all(lon_stations_all < 0) = lon_stations_all(lon_stations_all < 0) + 360; +elseif min_lon_ref == -180 + lon_stations_all(lon_stations_all > 180) = lon_stations_all(lon_stations_all < 180) - 360; +end +lat_stations_all = [lat_stations_RC4US; lat_stations_glorich'; lat_stations_arcticgro']; +Q_ann_all = [Q_ann_RC4US; Q_ann_glorich'; Q_ann_arcticgro']; + +dic_ann_all = [dic_ann_RC4US; dic_ann_glorich'; dic_ann_arcticgro']; +alk_ann_all = [alk_ann_RC4US; alk_ann_glorich'; alk_ann_arcticgro']; +no3_ann_all = [no3_ann_RC4US; no3_ann_glorich'; no3_ann_arcticgro']; +nh4_ann_all = [nh4_ann_RC4US; nh4_ann_glorich'; nh4_ann_arcticgro']; +din_ann_all = [din_ann_RC4US; din_ann_glorich'; din_ann_arcticgro']; +don_ann_all = [don_ann_RC4US; don_ann_glorich'; don_ann_arcticgro']; +pn_ann_all = [pn_ann_RC4US; pn_ann_glorich'; pn_ann_arcticgro']; +dip_ann_all = [dip_ann_RC4US; dip_ann_glorich'; dip_ann_arcticgro']; +dop_ann_all = [dop_ann_RC4US; dop_ann_glorich'; dop_ann_arcticgro']; +pp_ann_all = [pp_ann_RC4US; pp_ann_glorich'; pp_ann_arcticgro']; +dfe_ann_all = [dfe_ann_RC4US; dfe_ann_glorich'; dfe_ann_arcticgro']; +pfe_ann_all = [pfe_ann_RC4US; pfe_ann_glorich'; pfe_ann_arcticgro']; +si_ann_all = [si_ann_RC4US; si_ann_glorich'; si_ann_arcticgro']; +o2_ann_all = [o2_ann_RC4US; o2_ann_glorich'; o2_ann_arcticgro']; + +dic_monthly_all = [dic_monthly_RC4US dic_monthly_glorich dic_monthly_arcticgro]; +alk_monthly_all = [alk_monthly_RC4US alk_monthly_glorich alk_monthly_arcticgro]; +no3_monthly_all = [no3_monthly_RC4US no3_monthly_glorich no3_monthly_arcticgro]; +nh4_monthly_all = [nh4_monthly_RC4US nh4_monthly_glorich nh4_monthly_arcticgro]; +din_monthly_all = [din_monthly_RC4US din_monthly_glorich din_monthly_arcticgro]; +don_monthly_all = [don_monthly_RC4US don_monthly_glorich don_monthly_arcticgro]; +pn_monthly_all = [pn_monthly_RC4US pn_monthly_glorich pn_monthly_arcticgro]; +dip_monthly_all = [dip_monthly_RC4US dip_monthly_glorich dip_monthly_arcticgro]; +dop_monthly_all = [dop_monthly_RC4US dop_monthly_glorich dop_monthly_arcticgro]; +pp_monthly_all = [pp_monthly_RC4US pp_monthly_glorich pp_monthly_arcticgro]; +dfe_monthly_all = [dfe_monthly_RC4US dfe_monthly_glorich dfe_monthly_arcticgro]; +pfe_monthly_all = [pfe_monthly_RC4US pfe_monthly_glorich pfe_monthly_arcticgro]; +si_monthly_all = [si_monthly_RC4US si_monthly_glorich si_monthly_arcticgro]; +o2_monthly_all = [o2_monthly_RC4US o2_monthly_glorich o2_monthly_arcticgro]; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Load in monthly climatology of river forcing from the regional grid. % +% File contains: % +% runoff: monthly average runoff in kg m-2 sec-1 % +% area: area of grid cell in m-2 % +% lon: longitude (0-360 degrees) % +% lat: latitude % % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +load glofas_hill_runoff_monthlyclim_NEP10k_05122023.mat; + +lon_mod = lon; +lat_mod = lat; +area_mod = area; +% convert runoff from kg m-2 sec-1 to m3 sec-1 +Q_mod_monthly = zeros(size(runoff)); +for m = 1:12 + Q_mod_monthly(m,:,:) = squeeze(runoff(m,:,:)).*area_mod./1000; +end +Q_mod_ann = squeeze(mean(Q_mod_monthly,1)); +clear lon lat runoff area; + +%grid_file = '/archive/cas/Regional_MOM6/NWA12/nwa12_ocean_static.nc'; +%temp = ncread(grid_file,'deptho'); +%depth = permute(temp,[2,1]); clear temp; +%depth(isnan(depth)) = -1; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Filter for rivers in the region, set thresholds for minimum river size, % +% set parameters for plotting routines. % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% use grid to filter rivers outside domain +lat_mod_max = max(lat_mod(:)); +lat_mod_min = min(lat_mod(:)); +lon_mod_max = max(lon_mod(:)); +lon_mod_min = min(lon_mod(:)); + +in_region = find( (lon_stations_all <= lon_mod_max) & (lon_stations_all >= lon_mod_min) & ... + (lat_stations_all <= lat_mod_max) & (lat_stations_all >= lat_mod_min) & ... + (isfinite(Q_ann_all)) & (Q_ann_all > Q_min) ); +num_rivers = size(in_region,1); + +station_names_reg = station_names_all(in_region); +lon_stations_reg = lon_stations_all(in_region); +lat_stations_reg = lat_stations_all(in_region); +Q_ann_reg = Q_ann_all(in_region); +dic_ann_reg = dic_ann_all(in_region); +alk_ann_reg = alk_ann_all(in_region); +no3_ann_reg = no3_ann_all(in_region); +nh4_ann_reg = nh4_ann_all(in_region); +din_ann_reg = din_ann_all(in_region); +don_ann_reg = don_ann_all(in_region); +pn_ann_reg = pn_ann_all(in_region); +dip_ann_reg = dip_ann_all(in_region); +dop_ann_reg = dop_ann_all(in_region); +pp_ann_reg = pp_ann_all(in_region); +dfe_ann_reg = dfe_ann_all(in_region); +pfe_ann_reg = pfe_ann_all(in_region); +si_ann_reg = si_ann_all(in_region); +o2_ann_reg = o2_ann_all(in_region); + +for m = 1:12 + dic_monthly_reg(m,:) = dic_monthly_all(m,in_region); + alk_monthly_reg(m,:) = alk_monthly_all(m,in_region); + no3_monthly_reg(m,:) = no3_monthly_all(m,in_region); + nh4_monthly_reg(m,:) = nh4_monthly_all(m,in_region); + din_monthly_reg(m,:) = din_monthly_all(m,in_region); + don_monthly_reg(m,:) = don_monthly_all(m,in_region); + pn_monthly_reg(m,:) = pn_monthly_all(m,in_region); + dip_monthly_reg(m,:) = dip_monthly_all(m,in_region); + dop_monthly_reg(m,:) = dop_monthly_all(m,in_region); + pp_monthly_reg(m,:) = pp_monthly_all(m,in_region); + dfe_monthly_reg(m,:) = dfe_monthly_all(m,in_region); + pfe_monthly_reg(m,:) = pfe_monthly_all(m,in_region); + si_monthly_reg(m,:) = si_monthly_all(m,in_region); + o2_monthly_reg(m,:) = o2_monthly_all(m,in_region); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Assigning outflow points to rivers. % +% 1. Assignment starts with the rivers with the smallest flow and works % +% to the largest, w/larger river characteristics taking precedence to % +% ensure the most significant rivers are well represented. % +% 2. The algorithm keeps choosing the closest points to each river mouth % +% until the assigned flow is as close as possible to that observed % +% 3. Once the outflow points are assigned using the mean flow values, % +% monthly concentrations are assigned to those points. % +% 4. A simple "nearest neighbor" algorithm is used to fill in the gaps % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Sort rivers by discharge +[Q_ann_sort,sort_ind] = sort(Q_ann_reg,'ascend'); + +station_names_sort = station_names_reg(sort_ind); +lon_stations_sort = lon_stations_reg(sort_ind); +lat_stations_sort = lat_stations_reg(sort_ind); +Q_ann_sort = Q_ann_reg(sort_ind); +dic_ann_sort = dic_ann_reg(sort_ind); +alk_ann_sort = alk_ann_reg(sort_ind); +no3_ann_sort = no3_ann_reg(sort_ind); +nh4_ann_sort = nh4_ann_reg(sort_ind); +din_ann_sort = din_ann_reg(sort_ind); +don_ann_sort = don_ann_reg(sort_ind); +pn_ann_sort = pn_ann_reg(sort_ind); +dip_ann_sort = dip_ann_reg(sort_ind); +dop_ann_sort = dop_ann_reg(sort_ind); +pp_ann_sort = pp_ann_reg(sort_ind); +dfe_ann_sort = dfe_ann_reg(sort_ind); +pfe_ann_sort = pfe_ann_reg(sort_ind); +si_ann_sort = si_ann_reg(sort_ind); +o2_ann_sort = o2_ann_reg(sort_ind); + +for m = 1:12 + dic_monthly_sort(m,:) = dic_monthly_reg(m,sort_ind); + alk_monthly_sort(m,:) = alk_monthly_reg(m,sort_ind); + no3_monthly_sort(m,:) = no3_monthly_reg(m,sort_ind); + nh4_monthly_sort(m,:) = nh4_monthly_reg(m,sort_ind); + din_monthly_sort(m,:) = din_monthly_reg(m,sort_ind); + don_monthly_sort(m,:) = don_monthly_reg(m,sort_ind); + pn_monthly_sort(m,:) = pn_monthly_reg(m,sort_ind); + dip_monthly_sort(m,:) = dip_monthly_reg(m,sort_ind); + dop_monthly_sort(m,:) = dop_monthly_reg(m,sort_ind); + pp_monthly_sort(m,:) = pp_monthly_reg(m,sort_ind); + dfe_monthly_sort(m,:) = dfe_monthly_reg(m,sort_ind); + pfe_monthly_sort(m,:) = pfe_monthly_reg(m,sort_ind); + si_monthly_sort(m,:) = si_monthly_reg(m,sort_ind); + o2_monthly_sort(m,:) = o2_monthly_reg(m,sort_ind); +end + +% Create vectors of values at the runoff points from the model grid. These +% are used to accelerate the mapping relative to wrangling the full grid +% with all the zeros included. "ind_ro" are the grid indexes with runoff +ind_ro = find(Q_mod_ann > 0); +Q_mod_vec = Q_mod_ann(ind_ro); +lon_mod_runoff_vec = lon_mod(ind_ro); +lat_mod_runoff_vec = lat_mod(ind_ro); +Q_mod_monthly_vecs = zeros(12,size(lon_mod_runoff_vec,1)); +for m = 1:12 + temp = squeeze(Q_mod_monthly(m,:,:)); + Q_mod_monthly_vecs(m,:) = temp(ind_ro); +end + +% Create a grid of saturated oxygen values using the world ocean atlas data +temp_woa_monthly_vecs = zeros(12,size(lon_mod_runoff_vec,1)); +o2sat_woa_monthly_vecs = zeros(12,size(lon_mod_runoff_vec,1)); + +% Constants for o2 saturation calculation (taken from COBALT) +a_0 = 2.00907; +a_1 = 3.22014; +a_2 = 4.05010; +a_3 = 4.94457; +a_4 = -2.56847e-1; +a_5 = 3.88767; +sal = 0; +b_0 = -6.24523e-3; +b_1 = -7.37614e-3; +b_2 = -1.03410e-2; +b_3 = -8.17083e-3; +c_0 = -4.88682e-7; + +for m = 1:12 + temp = squeeze(woa_temp(m,:,:)); + % limit for validity of o2sat calculation + temp(temp > 40) = 40; temp(temp < 0) = 0; + temp_woa_monthly_vecs(m,:) = temp(ind_ro); + % calculate the oxygen saturation at a given temperature and salinity = 0 + % code taken from COBALT with limits applied above + tt = 298.15 - temp_woa_monthly_vecs(m,:); + tkb = 273.15 + temp_woa_monthly_vecs(m,:); + ts = log(tt / tkb); + ts2 = ts * ts; + ts3 = ts2 * ts; + ts4 = ts3 * ts; + ts5 = ts4 * ts; + + o2sat_woa_monthly_vecs(m,:) = (1000.0/22391.6) * 1000 * ... %convert from ml/l to mmol m-3 + exp(a_0 + a_1*ts + a_2*ts2 + a_3*ts3 + a_4*ts4 + a_5*ts5 + ... + (b_0 + b_1*ts + b_2*ts2 + b_3*ts3 + c_0*sal)*sal); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Load in fields generated from global NEWS. Where necessary, the ratio % +% of constituents relative to DIN will be used to fill forcing gaps % +% The m-file used to generate the NEWS forcing file is included in this % +% directory and uses an analogous mapping algorithm to this one % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +din_ann_NEWS = nc_varget(NEWS_file,'NO3_CONC'); +aa = find(din_ann_NEWS > 0); +temp1 = nc_varget(NEWS_file,'LDON_CONC'); +temp2 = nc_varget(NEWS_file,'SLDON_CONC'); +temp3 = nc_varget(NEWS_file,'SRDON_CONC'); +don_ann_NEWS = temp1 + temp2 + temp3; +don_ratio_NEWS_vec = don_ann_NEWS(aa)./din_ann_NEWS(aa); +clear temp1 temp2 temp3; +pn_ann_NEWS = nc_varget(NEWS_file,'NDET_CONC'); +pn_ratio_NEWS_vec = pn_ann_NEWS(aa)./din_ann_NEWS(aa); + +dip_ann_NEWS = nc_varget(NEWS_file,'PO4_CONC'); +dip_ratio_NEWS_vec = dip_ann_NEWS(aa)./din_ann_NEWS(aa); +temp1 = nc_varget(NEWS_file,'LDOP_CONC'); +temp2 = nc_varget(NEWS_file,'SLDOP_CONC'); +temp3 = nc_varget(NEWS_file,'SRDOP_CONC'); +dop_ann_NEWS = temp1 + temp2 + temp3; +dop_ratio_NEWS_vec = dop_ann_NEWS(aa)./din_ann_NEWS(aa); +clear temp1 temp2 temp3; +pp_ann_NEWS = nc_varget(NEWS_file,'PDET_CONC'); +pp_ratio_NEWS_vec = pp_ann_NEWS(aa)./din_ann_NEWS(aa); +si_ann_NEWS = nc_varget(NEWS_file,'SI_CONC'); +si_ratio_NEWS_vec = si_ann_NEWS(aa)./din_ann_NEWS(aa); + +% Vectors to hold monthly values mapped onto model runoff points +dic_mod_monthly_vecs = zeros(12,size(lon_mod_runoff_vec,1)); +alk_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +no3_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +nh4_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +din_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +don_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +pn_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +dip_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +dop_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +pp_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +dfe_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +pfe_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +si_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); +o2_mod_monthly_vecs = zeros(12,size(lat_mod_runoff_vec,1)); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Loop identifies points assigned to each river % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +for k=1:num_rivers + dist = pdist2([lon_stations_sort(k) lat_stations_sort(k)], ... + [lon_mod_runoff_vec lat_mod_runoff_vec]); + [dist_sort, dist_sort_ind] = sort(dist,'ascend'); + + if dist_sort(1) < min_dist + Q_sum1 = 0; + Q_sum2 = 0; + n = 0; + while (Q_sum2 < Q_ann_sort(k) && (dist_sort(n+1) < max_dist)) + Q_sum1 = Q_sum2; + n = n+1; + Q_sum2 = Q_sum1 + Q_mod_vec(dist_sort_ind(n)); + end + %if abs(Q_sum1 - Q_ann_sort(k)) < abs(Q_sum2 - Q_ann_sort(k)) + % nrp = n-1; % number of runoff points + % [Q_sum1 Q_ann_sort(k)] % a quick check for comparable flow + %else + nrp = n; + [Q_sum2 Q_ann_sort(k)] + %end + + % enter monthly concentration values into an array of monthly values + % if no monthly value is available, use annuals. + for m = 1:12 + + if isnan(dic_monthly_sort(m,k)) + if isfinite(dic_ann_sort(k)) + dic_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dic_ann_sort(k); + else + dic_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = 0; + end + else + dic_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dic_monthly_sort(m,k); + end + + if isnan(alk_monthly_sort(m,k)) + if isfinite(dic_ann_sort(k)) + alk_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = alk_ann_sort(k); + else + alk_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = 0; + end + else + alk_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = alk_monthly_sort(m,k); + end + + % mapping assumes that DIN is defined for nutrient calculations since + % ratios relative to DIN are used to fill in other components. If + % DIN is not defined, values are left at 0 and eventually filled with + % a nearest neighbor filling (next section) + + if (isfinite(din_monthly_sort(m,k)) || isfinite(din_ann_sort(k)) ) + if isfinite(din_monthly_sort(m,k)) + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = din_monthly_sort(m,k); + elseif isfinite(din_ann_sort(k)) + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = din_ann_sort(k); + end + + if isfinite(no3_monthly_sort(m,k)) + no3_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = no3_monthly_sort(m,k); + elseif isfinite(no3_ann_sort(k)) + no3_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = no3_ann_sort(k); + end + + if isfinite(nh4_monthly_sort(m,k)) + nh4_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = nh4_monthly_sort(m,k); + elseif isfinite(nh4_ann_sort(k)) + nh4_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = nh4_ann_sort(k); + end + + if (isnan(don_monthly_sort(m,k)) && isnan(don_ann_sort(k))) + don_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + don_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(don_monthly_sort(m,k)) && ~isnan(don_ann_sort(k))) + don_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = don_ann_sort(k); + else + don_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = don_monthly_sort(m,k); + end + + if (isnan(pn_monthly_sort(m,k)) && isnan(pn_ann_sort(k))) + pn_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + pn_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(pn_monthly_sort(m,k)) && ~isnan(pn_ann_sort(k))) + pn_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = pn_ann_sort(k); + else + pn_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = pn_monthly_sort(m,k); + end + + if (isnan(dip_monthly_sort(m,k)) && isnan(dip_ann_sort(k))) + dip_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + dip_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(dip_monthly_sort(m,k)) && ~isnan(dip_ann_sort(k))) + dip_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dip_ann_sort(k); + else + dip_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dip_monthly_sort(m,k); + end + + if (isnan(dop_monthly_sort(m,k)) && isnan(dop_ann_sort(k))) + dop_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + dop_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(dop_monthly_sort(m,k)) && ~isnan(dop_ann_sort(k))) + dop_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dop_ann_sort(k); + else + dop_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = dop_monthly_sort(m,k); + end + + if (isnan(pp_monthly_sort(m,k)) && isnan(pp_ann_sort(k))) + pp_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + pp_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(pp_monthly_sort(m,k)) && ~isnan(pp_ann_sort(k))) + pp_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = pp_ann_sort(k); + else + pp_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = pp_monthly_sort(m,k); + end + + if (isnan(si_monthly_sort(m,k)) && isnan(si_ann_sort(k))) + si_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + din_mod_monthly_vecs(m,dist_sort_ind(1:nrp)).* ... + si_ratio_NEWS_vec(dist_sort_ind(1:nrp))'; + elseif (isnan(si_monthly_sort(m,k)) && ~isnan(si_ann_sort(k))) + si_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = si_ann_sort(k); + else + si_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = si_monthly_sort(m,k); + end + + if isnan(o2_monthly_sort(m,k)) + o2_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = ... + o2sat_woa_monthly_vecs(m,dist_sort_ind(1:nrp)); + else + o2_mod_monthly_vecs(m,dist_sort_ind(1:nrp)) = o2_monthly_sort(m,k); + end + + end + end + + % plot to check location if inspect_map == 'y'. The plot puts + % open circles at each runoff location in the model grid and fills + % those that are assigned to each river. Note that some of the smaller + % rivers may be replaced with larger ones as the fitting process + % continues. + + if inspect_map == 'y' + figure(1) + clf + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(Q_mod_vec),3,log10(Q_mod_vec)); + hold on + scatter3(lon_mod_runoff_vec(dist_sort_ind(1:nrp)),lat_mod_runoff_vec(dist_sort_ind(1:nrp)), ... + log10(Q_mod_vec(dist_sort_ind(1:nrp))),25, ... + log10(Q_mod_vec(dist_sort_ind(1:nrp))),'filled'); + view(2); + plot3(lon_stations_sort(k),lat_stations_sort(k),1e5,'k.','MarkerSize',20); + %contour(lon_mod,lat_mod,depth,[0 0],'k-'); + axis([lon_stations_sort(k)-plot_width/2 lon_stations_sort(k)+plot_width/2 ... + lat_stations_sort(k)-plot_width/2 lat_stations_sort(k)+plot_width/2]); + caxis([-4 3]); + titl = ['river number: ',num2str(k),' name: ',station_names_sort{k}]; + title(titl); + colorbar; + + % check the values of each mapping + N_check = mean(din_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all') + ... + mean(don_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all') + ... + mean(pn_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + P_check = mean(dip_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all') + ... + mean(dop_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all') + ... + mean(pp_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + SI_check = mean(si_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + DIN_check = mean(din_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + DON_check = mean(don_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + PN_check = mean(pn_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + DIP_check = mean(dip_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + DOP_check = mean(dop_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + PP_check = mean(pp_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + SI_check = mean(si_mod_monthly_vecs(:,dist_sort_ind(1:nrp)),'all'); + + station_names_sort(k) + ind = dist_sort_ind(1:nrp); + 'total flow in m3 sec' + [Q_ann_sort(k) sum(Q_mod_vec(dist_sort_ind(1:nrp)))] + 'N, P conc (mmoles m-3), DI, DO, P' + [DIN_check DON_check PN_check] + [DIP_check DOP_check PP_check] + 'Total N, Total P, Total N: Total P' + [N_check P_check N_check/P_check] + 'DO:DI and P:DI ratios'; + [DON_check/DIN_check PN_check/DIN_check]; + [DOP_check/DIP_check PP_check/DIP_check]; + 'silica concentration (mmoles m-3)'; + SI_check; + + pause + end + + % If river is outside the domain, skip all of the calculations above and + % just plot for inspection/evaluation + else + % This is for rivers that were outside of the domain + if inspect_map == 'y' + figure(1) + clf + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,log10(Q_mod_vec),3,log10(Q_mod_vec)); + hold on + view(2); + plot3(lon_stations_sort(k),lat_stations_sort(k),1e5,'k.','MarkerSize',20); + axis([lon_stations_sort(k)-10 lon_stations_sort(k)+10 ... + lat_stations_sort(k)-10 lat_stations_sort(k)+10]); + caxis([-4 3]); + titl = ['OUTSIDE: river number: ',num2str(k),' name: ',station_names_sort{k}]; + title(titl); + colorbar; + + pause + end + + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% nearest neighbor search to fill in any runoff points that were not % +% assigned after the runoff mapping step % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +lon_mod_runoff_vec = double(lon_mod_runoff_vec); +lat_mod_runoff_vec = double(lat_mod_runoff_vec); + +for m = 1:12 + aa = find(dic_mod_monthly_vecs(m,:) == 0); + bb = find(dic_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + dic_mod_monthly_vecs(m,bb)','nearest','nearest'); + dic_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(alk_mod_monthly_vecs(m,:) == 0); + bb = find(alk_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + alk_mod_monthly_vecs(m,bb)','nearest','nearest'); + alk_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(no3_mod_monthly_vecs(m,:) == 0); + bb = find(no3_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + no3_mod_monthly_vecs(m,bb)','nearest','nearest'); + no3_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(nh4_mod_monthly_vecs(m,:) == 0); + bb = find(nh4_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + nh4_mod_monthly_vecs(m,bb)','nearest','nearest'); + nh4_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(din_mod_monthly_vecs(m,:) == 0); + bb = find(din_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + din_mod_monthly_vecs(m,bb)','nearest','nearest'); + din_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(don_mod_monthly_vecs(m,:) == 0); + bb = find(don_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + don_mod_monthly_vecs(m,bb)','nearest','nearest'); + don_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(pn_mod_monthly_vecs(m,:) == 0); + bb = find(pn_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + pn_mod_monthly_vecs(m,bb)','nearest','nearest'); + pn_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(dip_mod_monthly_vecs(m,:) == 0); + bb = find(dip_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + dip_mod_monthly_vecs(m,bb)','nearest','nearest'); + dip_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(dop_mod_monthly_vecs(m,:) == 0); + bb = find(dop_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + dop_mod_monthly_vecs(m,bb)','nearest','nearest'); + dop_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(pp_mod_monthly_vecs(m,:) == 0); + bb = find(pp_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + pp_mod_monthly_vecs(m,bb)','nearest','nearest'); + pp_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); + + aa = find(si_mod_monthly_vecs(m,:) == 0); + bb = find(si_mod_monthly_vecs(m,:) > 0); + F = scatteredInterpolant(lon_mod_runoff_vec(bb),lat_mod_runoff_vec(bb), ... + si_mod_monthly_vecs(m,bb)','nearest','nearest'); + si_mod_monthly_vecs(m,aa) = F(lon_mod_runoff_vec(aa),lat_mod_runoff_vec(aa)); +end + +% For o2sat, fill in any 0 values with saturated o2 at the world ocean +% atlas climatology +for m = 1:12 + aa = find(o2_mod_monthly_vecs(m,:) == 0); + o2_mod_monthly_vecs(m,aa) = o2sat_woa_monthly_vecs(m,aa); +end + +totn_mod_monthly_vecs = din_mod_monthly_vecs + don_mod_monthly_vecs + pn_mod_monthly_vecs; +totp_mod_monthly_vecs = dip_mod_monthly_vecs + dop_mod_monthly_vecs + pp_mod_monthly_vecs; + +dicflux_mod_monthly_vecs = dic_mod_monthly_vecs.*Q_mod_monthly_vecs; +alkflux_mod_monthly_vecs = alk_mod_monthly_vecs.*Q_mod_monthly_vecs; +dinflux_mod_monthly_vecs = din_mod_monthly_vecs.*Q_mod_monthly_vecs; +no3flux_mod_monthly_vecs = no3_mod_monthly_vecs.*Q_mod_monthly_vecs; +nh4flux_mod_monthly_vecs = nh4_mod_monthly_vecs.*Q_mod_monthly_vecs; +dipflux_mod_monthly_vecs = dip_mod_monthly_vecs.*Q_mod_monthly_vecs; +donflux_mod_monthly_vecs = don_mod_monthly_vecs.*Q_mod_monthly_vecs; +dopflux_mod_monthly_vecs = dop_mod_monthly_vecs.*Q_mod_monthly_vecs; +pnflux_mod_monthly_vecs = pn_mod_monthly_vecs.*Q_mod_monthly_vecs; +ppflux_mod_monthly_vecs = pp_mod_monthly_vecs.*Q_mod_monthly_vecs; +siflux_mod_monthly_vecs = si_mod_monthly_vecs.*Q_mod_monthly_vecs; +totnflux_mod_monthly_vecs = totn_mod_monthly_vecs.*Q_mod_monthly_vecs; +totpflux_mod_monthly_vecs = totp_mod_monthly_vecs.*Q_mod_monthly_vecs; +o2flux_mod_monthly_vecs = o2_mod_monthly_vecs.*Q_mod_monthly_vecs; + +dicflux_mod_ann_vec = mean(dicflux_mod_monthly_vecs,1); +alkflux_mod_ann_vec = mean(alkflux_mod_monthly_vecs,1); +dinflux_mod_ann_vec = mean(dinflux_mod_monthly_vecs,1); +no3flux_mod_ann_vec = mean(no3flux_mod_monthly_vecs,1); +nh4flux_mod_ann_vec = mean(nh4flux_mod_monthly_vecs,1); +dipflux_mod_ann_vec = mean(dipflux_mod_monthly_vecs,1); +donflux_mod_ann_vec = mean(donflux_mod_monthly_vecs,1); +dopflux_mod_ann_vec = mean(dopflux_mod_monthly_vecs,1); +pnflux_mod_ann_vec = mean(pnflux_mod_monthly_vecs,1); +ppflux_mod_ann_vec = mean(ppflux_mod_monthly_vecs,1); +siflux_mod_ann_vec = mean(siflux_mod_monthly_vecs,1); +totnflux_mod_ann_vec = mean(totnflux_mod_monthly_vecs,1); +totpflux_mod_ann_vec = mean(totpflux_mod_monthly_vecs,1); +o2flux_mod_ann_vec = mean(o2flux_mod_monthly_vecs,1); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Produce plots to evaluate the mapping % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% scale marker size with the freshwater flux +ms_vec = zeros(size(Q_mod_vec)); +ms_vec(log10(Q_mod_vec) < 0) = 1; +ms_vec((log10(Q_mod_vec) > 0) & (log10(Q_mod_vec) < 1)) = 2.5; +ms_vec((log10(Q_mod_vec) > 1) & (log10(Q_mod_vec) < 2)) = 10; +ms_vec((log10(Q_mod_vec) > 2) & (log10(Q_mod_vec) < 3)) = 25; +ms_vec(log10(Q_mod_vec) > 3) = 100; + +% DIC, Alk concentrations and DIC:Alk +for m = 1:12 + + figure(1) + clf + + subplot(1,3,1); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dic_mod_monthly_vecs(m,:),ms_vec, ... + dic_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 2500]); + colorbar + titlestr = ['DIC, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(1,3,2); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,alk_mod_monthly_vecs(m,:),ms_vec, ... + alk_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 2500]); + colorbar + titlestr = ['Alk, meq m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(1,3,3) + dic_alk_ratio = dic_mod_monthly_vecs(m,:)./alk_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dic_alk_ratio, ... + ms_vec,dic_alk_ratio,'filled'); + hold on + view(2); + caxis([0.8 1.2]); + colorbar + title('DIC:Alk ratio'); + + pause +end + +% Nitrogen Concentrations +for m = 1:12 + + figure(1) + clf + + subplot(2,2,1); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,din_mod_monthly_vecs(m,:),ms_vec, ... + din_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 100]); + colorbar + titlestr = ['DIN, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,2,2); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,no3_mod_monthly_vecs(m,:),ms_vec, ... + din_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 100]); + colorbar + titlestr = ['no3, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,2,3); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,nh4_mod_monthly_vecs(m,:),ms_vec, ... + nh4_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 20]); + colorbar + titlestr = ['nh4, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,2,4); + no3_din_ratio = no3_mod_monthly_vecs(m,:)./din_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,no3_din_ratio, ... + ms_vec,no3_din_ratio,'filled'); + hold on + view(2); + caxis([0 1.0]); + colorbar + title('NO3:DIN ratio'); + + pause +end + +for m = 1:12 + + subplot(2,3,1); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,din_mod_monthly_vecs(m,:),ms_vec, ... + din_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 100]); + colorbar + titlestr = ['DIN, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,2); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,don_mod_monthly_vecs(m,:),ms_vec, ... + don_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 100]); + colorbar + titlestr = ['DON, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,3); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pn_mod_monthly_vecs(m,:),ms_vec, ... + pn_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 100]); + colorbar + titlestr = ['PN, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,5) + don_din_ratio = don_mod_monthly_vecs(m,:)./din_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,don_din_ratio, ... + ms_vec,don_din_ratio,'filled'); + hold on + view(2); + caxis([0 2]); + colorbar + title('DON:DIN ratio'); + + subplot(2,3,6) + pn_din_ratio = pn_mod_monthly_vecs(m,:)./din_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pn_din_ratio, ... + ms_vec,pn_din_ratio,'filled'); + hold on + view(2); + caxis([0 2]); + colorbar + title('PN:DIN ratio'); + + pause +end + +% Phosphorus Concentrations +for m = 1:12 + + figure(1) + clf + + subplot(2,3,1); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dip_mod_monthly_vecs(m,:),ms_vec, ... + dip_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 3]); + colorbar + titlestr = ['DIP, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,2); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dop_mod_monthly_vecs(m,:),ms_vec, ... + dop_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 3]); + colorbar + titlestr = ['DOP, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,3); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pp_mod_monthly_vecs(m,:),ms_vec, ... + pp_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 3]); + colorbar + titlestr = ['PP, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,3,5) + dop_dip_ratio = dop_mod_monthly_vecs(m,:)./dip_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,dop_dip_ratio, ... + ms_vec,dop_dip_ratio,'filled'); + hold on + view(2); + caxis([0 3]); + colorbar + title('DOP:DIP ratio'); + + subplot(2,3,6) + pp_dip_ratio = pp_mod_monthly_vecs(m,:)./dip_mod_monthly_vecs(m,:); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,pp_dip_ratio, ... + ms_vec,pp_dip_ratio,'filled'); + hold on + view(2); + caxis([0 2]); + colorbar + title('PP:DIP ratio'); + + pause +end + +% silica and oxygen concentrations +for m = 1:12 + figure(1) + clf + + subplot(2,1,1); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,si_mod_monthly_vecs(m,:),ms_vec, ... + si_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 200]); + colorbar + titlestr = ['Si, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + subplot(2,1,2); + scatter3(lon_mod_runoff_vec,lat_mod_runoff_vec,o2_mod_monthly_vecs(m,:),ms_vec, ... + o2_mod_monthly_vecs(m,:),'filled'); + hold on + view(2); + caxis([0 350]); + colorbar + titlestr = ['o2, mmoles m-3; month = ',num2str(m)]; + title(titlestr); + + pause +end + +% Initialize 2D concentration arrays; these are the ones read into MOM6 to +% specify the nutrient concentrations of river inputs. +DIC_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +ALK_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +NO3_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +NH4_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +LDON_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +SLDON_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +SRDON_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +PO4_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +LDOP_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +SLDOP_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +SRDOP_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +NDET_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +PDET_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +SI_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); +O2_CONC = zeros(12,size(lon_mod,1),size(lon_mod,2)); + +% Map concentration vectors onto 2D arrays +temp = zeros(size(lon_mod)); +for m = 1:12 + temp(ind_ro) = dic_mod_monthly_vecs(m,:); + DIC_CONC(m,:,:) = temp; + temp(ind_ro) = alk_mod_monthly_vecs(m,:); + ALK_CONC(m,:,:) = temp; + + temp(ind_ro) = no3_mod_monthly_vecs(m,:); % contains all NEWS DIN + %temp(ind_ro) = din_mod_monthly_vecs(m,:); + NO3_CONC(m,:,:) = temp; + temp(ind_ro) = nh4_mod_monthly_vecs(m,:); + NH4_CONC(m,:,:) = temp; + temp(ind_ro) = don_mod_monthly_vecs(m,:); + LDON_CONC(m,:,:) = frac_ldon*temp; + SLDON_CONC(m,:,:) = frac_sldon*temp; + SRDON_CONC(m,:,:) = frac_srdon*temp; + temp(ind_ro) = pn_mod_monthly_vecs(m,:); + NDET_CONC(m,:,:) = temp; + + temp(ind_ro) = dip_mod_monthly_vecs(m,:); + PO4_CONC(m,:,:) = temp; + temp(ind_ro) = dop_mod_monthly_vecs(m,:); + LDOP_CONC(m,:,:) = frac_ldop*temp; + SLDOP_CONC(m,:,:) = frac_sldop*temp; + SRDOP_CONC(m,:,:) = frac_srdop*temp; + temp(ind_ro) = pp_mod_monthly_vecs(m,:); + PDET_CONC(m,:,:) = temp; + + temp(ind_ro) = si_mod_monthly_vecs(m,:); + SI_CONC(m,:,:) = temp; + temp(ind_ro) = o2_mod_monthly_vecs(m,:); + O2_CONC(m,:,:) = temp; +end + +% MOM6 is taking river values in moles m-3 for other constituents. Change +% for consistency across river constituents +DIC_CONC = DIC_CONC./1e3; +ALK_CONC = ALK_CONC./1e3; +NO3_CONC = NO3_CONC./1e3; +NH4_CONC = NH4_CONC./1e3; +LDON_CONC = LDON_CONC./1e3; +SLDON_CONC = SLDON_CONC./1e3; +SRDON_CONC = SRDON_CONC./1e3; +NDET_CONC = NDET_CONC./1e3; +PO4_CONC = PO4_CONC./1e3; +LDOP_CONC = LDON_CONC./1e3; +SLDOP_CONC = SLDON_CONC./1e3; +SRDOP_CONC = SRDON_CONC./1e3; +PDET_CONC = PDET_CONC./1e3; +SI_CONC = SI_CONC./1e3; +O2_CONC = O2_CONC./1e3; + +% Add iron concentrations - initialize with nitrate and then overwrite +FED_CONC = NO3_CONC; +FEDET_CONC = NO3_CONC; +% 40 nM dissolved iron concentration from De Baar and De Jong + 30nM +% Colloidal and nanoparticle flux as reported in Canfield and Raiswell +FED_CONC(FED_CONC > 0) = const_fed; +FEDET_CONC(FEDET_CONC > 0) = 0.0; + +ms = 8; + +% quick set of plots to ensure the mapping to the model grid was successful +for m = 1:12 + % DIC and alkalinity + figure(1) + clf + subplot(2,1,1); + title('log10(DIC CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(DIC_CONC(m,:)),ms,log10(DIC_CONC(m,:)),'filled'); + caxis([-1 1]); colorbar; + + subplot(2,1,2); + title('log10(ALK CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(ALK_CONC(m,:)),ms,log10(ALK_CONC(m,:)),'filled'); + caxis([-1 1]); colorbar; + + %pause +end + +% Nitrogen +for m = 1:12 + figure(1); + clf + subplot(3,2,1); + title('log10(NO3 CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(NO3_CONC(m,:)),ms,log10(NO3_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + subplot(3,2,2); + title('log10(NH4 CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(NH4_CONC(m,:)),ms,log10(NH4_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + subplot(3,2,3); + title('log10(LDON CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(LDON_CONC(m,:)),ms,log10(LDON_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + subplot(3,2,4); + title('log10(SLDON CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(SLDON_CONC(m,:)),ms,log10(SLDON_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + subplot(3,2,5); + title('log10(SRDON CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(SRDON_CONC(m,:)),ms,log10(SRDON_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + subplot(3,2,6); + title('log10(NDET CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(NDET_CONC(m,:)),ms,log10(NDET_CONC(m,:)),'filled'); + caxis([-4 -1]); colorbar; + + %pause +end + +% Phosphorus +for m = 1:12 + figure(1); + clf + subplot(3,2,1); + title('log10(PO4 CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(PO4_CONC(m,:)),ms,log10(PO4_CONC(m,:)),'filled'); + caxis([-4 -2]); colorbar; + + subplot(3,2,2); + title('log10(LDOP CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(LDOP_CONC(m,:)),ms,log10(LDOP_CONC(m,:)),'filled'); + caxis([-4 -2]); colorbar; + + subplot(3,2,3); + title('log10(SLDOP CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(SLDOP_CONC(m,:)),ms,log10(SLDOP_CONC(m,:)),'filled'); + caxis([-4 -2]); colorbar; + + subplot(3,2,4); + title('log10(SRDOP CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(SRDOP_CONC(m,:)),ms,log10(SRDOP_CONC(m,:)),'filled'); + caxis([-4 -2]); colorbar; + + subplot(3,2,5); + title('log10(PDET CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(PDET_CONC(m,:)),ms,log10(PDET_CONC(m,:)),'filled'); + caxis([-4 -2]); colorbar; + + %pause; +end + +% Iron, Silica, Oxygen +for m = 1:12 + figure(1) + clf + subplot(3,2,1); + title('log10(FED CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(FED_CONC(m,:)),ms,log10(FED_CONC(m,:)),'filled'); + caxis([-5 -3]); colorbar; + + subplot(3,2,2); + title('log10(FEDET CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(FEDET_CONC(m,:)),ms,log10(FEDET_CONC(m,:)),'filled'); + caxis([-5 -3]); colorbar; + + subplot(3,2,3); + title('log10(SI CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(SI_CONC(m,:)),ms,log10(SI_CONC(m,:)),'filled'); + caxis([-3 0]); colorbar; + + subplot(3,2,4); + title('log10(O2 CONC)'); hold on; + scatter3(lon_mod(:),lat_mod(:),log10(O2_CONC(m,:)),ms,log10(O2_CONC(m,:)),'filled'); + caxis([-3 0]); colorbar; + + %pause; +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Save Files % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% option to save matlab file +% save River_DIC_ALK_RC4US_NWA ALK_CONC DIC_CONC + +% Construct netcdf file following format used by nutrient input files to +% MOM6 + +% Make reference dates for standard non-leap year +dates = [1990 1 16 12 0 0; 1990 2 15 0 0 0; 1990 3 16 12 0 0; ... + 1990 4 16 0 0 0; 1990 5 16 12 0 0; 1990 6 16 0 0 0; ... + 1990 7 16 12 0 0; 1990 8 16 12 0 0; 1990 9 16 0 0 0; ... + 1990 10 16 12 0 0; 1990 11 16 0 0 0; 1990 12 16 12 0 0]; +time = datenum(dates) - datenum([1990 1 1 0 0 0]); + +nlat = size(lat_mod,1); +nlon = size(lat_mod,2); + +ncid = netcdf.create(nc_file_name,'CLOBBER'); +dimid0 = netcdf.defDim(ncid,'time',netcdf.getConstant('NC_UNLIMITED')); +dimid1 = netcdf.defDim(ncid,'y',nlat); +dimid2 = netcdf.defDim(ncid,'x',nlon); + +varid0 = netcdf.defVar(ncid,'time','double',dimid0); +netcdf.putAtt(ncid,varid0,'calendar','NOLEAP'); +netcdf.putAtt(ncid,varid0,'calendar_type','NOLEAP'); +netcdf.putAtt(ncid,varid0,'modulo','T'); +netcdf.putAtt(ncid,varid0,'units','days since 1990-1-1 0:00:00'); +netcdf.putAtt(ncid,varid0,'time_origin','01-JAN-1990 00:00:00'); +varid1 = netcdf.defVar(ncid,'y','int',dimid1); +netcdf.putAtt(ncid,varid1,'cartesian_axis','Y'); +varid2 = netcdf.defVar(ncid,'x','int',dimid2); +netcdf.putAtt(ncid,varid2,'cartesian_axis','X'); +varid3 = netcdf.defVar(ncid,'lat','double',[dimid2 dimid1]); +netcdf.putAtt(ncid,varid3,'units','degrees north'); +varid4 = netcdf.defVar(ncid,'lon','double',[dimid2 dimid1]); +netcdf.putAtt(ncid,varid4,'units','degrees east'); +varid5 = netcdf.defVar(ncid,'DIC_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid5,'units','mol m-3'); +netcdf.putAtt(ncid,varid5,'long_name','DIC_CONC'); +varid6 = netcdf.defVar(ncid,'ALK_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid6,'units','mole Eq. m-3'); +netcdf.putAtt(ncid,varid6,'long_name','ALK_CONC'); +varid7 = netcdf.defVar(ncid,'NO3_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid7,'units','mol m-3'); +netcdf.putAtt(ncid,varid7,'long_name','NO3_CONC'); +varid8 = netcdf.defVar(ncid,'NH4_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid8,'units','mol m-3'); +netcdf.putAtt(ncid,varid8,'long_name','NH4_CONC'); +varid9 = netcdf.defVar(ncid,'LDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid9,'units','mol m-3'); +netcdf.putAtt(ncid,varid9,'long_name','0.3*DON_CONC'); +varid10 = netcdf.defVar(ncid,'SLDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid10,'units','mol m-3'); +netcdf.putAtt(ncid,varid10,'long_name','0.35*DON_CONC'); +varid11 = netcdf.defVar(ncid,'SRDON_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid11,'units','mol m-3'); +netcdf.putAtt(ncid,varid11,'long_name','0.35*DON_CONC'); +varid12 = netcdf.defVar(ncid,'NDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid12,'units','mol m-3'); +netcdf.putAtt(ncid,varid12,'long_name','1.0*PN_CONC'); +varid13 = netcdf.defVar(ncid,'PO4_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid13,'units','mol m-3'); +netcdf.putAtt(ncid,varid13,'long_name','PO4_CONC+0.3*PP_CONC'); +varid14 = netcdf.defVar(ncid,'LDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid14,'units','mol m-3'); +netcdf.putAtt(ncid,varid14,'long_name','0.3*DOP_CONC'); +varid15 = netcdf.defVar(ncid,'SLDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid15,'units','mol m-3'); +netcdf.putAtt(ncid,varid15,'long_name','0.35*DOP_CONC'); +varid16 = netcdf.defVar(ncid,'SRDOP_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid16,'units','mol m-3'); +netcdf.putAtt(ncid,varid16,'long_name','0.35*DOP_CONC'); +varid17 = netcdf.defVar(ncid,'PDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid17,'units','mol m-3'); +netcdf.putAtt(ncid,varid17,'long_name','0*PP_CONC'); +varid18 = netcdf.defVar(ncid,'FED_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid18,'units','mol m-3'); +netcdf.putAtt(ncid,varid18,'long_name','FED_CONC'); +varid19 = netcdf.defVar(ncid,'FEDET_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid19,'units','mol m-3'); +netcdf.putAtt(ncid,varid19,'long_name','FEDET_CONC'); +varid20 = netcdf.defVar(ncid,'O2_CONC','double',[dimid2,dimid1,dimid0]); +netcdf.putAtt(ncid,varid20,'units','mol m-3'); +netcdf.putAtt(ncid,varid20,'long_name','O2_CONC'); +netcdf.close(ncid) + +ncid = netcdf.open(nc_file_name,'NC_WRITE'); +netcdf.putVar(ncid,varid0,0,12,time); +% nutrient input files appear seem to need dummy axes to be read in +% properly, but eventually do a grid by grid mapping that doesn't require +% these. +netcdf.putVar(ncid,varid1,1:nlat); +netcdf.putVar(ncid,varid2,1:nlon); +netcdf.putVar(ncid,varid3,permute(lat_mod,[2,1])); +netcdf.putVar(ncid,varid4,permute(lon_mod,[2,1])); +netcdf.putVar(ncid,varid5,permute(DIC_CONC,[3,2,1])); +netcdf.putVar(ncid,varid6,permute(ALK_CONC,[3,2,1])); +netcdf.putVar(ncid,varid7,permute(NO3_CONC,[3,2,1])); +netcdf.putVar(ncid,varid8,permute(NH4_CONC,[3,2,1])); +netcdf.putVar(ncid,varid9,permute(LDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid10,permute(SLDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid11,permute(SRDON_CONC,[3,2,1])); +netcdf.putVar(ncid,varid12,permute(NDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid13,permute(PO4_CONC,[3,2,1])); +netcdf.putVar(ncid,varid14,permute(LDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid15,permute(SLDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid16,permute(SRDOP_CONC,[3,2,1])); +netcdf.putVar(ncid,varid17,permute(PDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid18,permute(FED_CONC,[3,2,1])); +netcdf.putVar(ncid,varid19,permute(FEDET_CONC,[3,2,1])); +netcdf.putVar(ncid,varid20,permute(O2_CONC,[3,2,1])); +netcdf.close(ncid) diff --git a/tools/rivers/bgc/GlobalNEWS2_RH2000Dataset-version1.0.xls b/tools/rivers/bgc/NWA/GlobalNEWS2_RH2000Dataset-version1.0.xls similarity index 100% rename from tools/rivers/bgc/GlobalNEWS2_RH2000Dataset-version1.0.xls rename to tools/rivers/bgc/NWA/GlobalNEWS2_RH2000Dataset-version1.0.xls diff --git a/tools/rivers/bgc/README.md b/tools/rivers/bgc/NWA/README.md similarity index 80% rename from tools/rivers/bgc/README.md rename to tools/rivers/bgc/NWA/README.md index 434413a4f..6a14aa436 100644 --- a/tools/rivers/bgc/README.md +++ b/tools/rivers/bgc/NWA/README.md @@ -1,6 +1,6 @@ -## Example Scripts for BGC runoff generation +## Example Scripts for NWA BGC runoff generation -This folder contains example scripts for BGC runoff file generation. Users can follow the following instructions to generate BGC runoff file: +This folder contains example scripts for NWA BGC runoff file generation. Users can follow the following instructions to generate BGC runoff file: ``` # 1. First run write_glofas_ave.py to get climatological monthly runoff ./write_glofas_ave.py diff --git a/tools/rivers/bgc/mapriv_NEWS2_NWA12_GLOFAS.m b/tools/rivers/bgc/NWA/mapriv_NEWS2_NWA12_GLOFAS.m similarity index 100% rename from tools/rivers/bgc/mapriv_NEWS2_NWA12_GLOFAS.m rename to tools/rivers/bgc/NWA/mapriv_NEWS2_NWA12_GLOFAS.m diff --git a/tools/rivers/bgc/mapriv_combined_NWA12_GLOFAS.m b/tools/rivers/bgc/NWA/mapriv_combined_NWA12_GLOFAS.m similarity index 100% rename from tools/rivers/bgc/mapriv_combined_NWA12_GLOFAS.m rename to tools/rivers/bgc/NWA/mapriv_combined_NWA12_GLOFAS.m diff --git a/tools/rivers/bgc/nc64startup.m b/tools/rivers/bgc/NWA/nc64startup.m similarity index 100% rename from tools/rivers/bgc/nc64startup.m rename to tools/rivers/bgc/NWA/nc64startup.m diff --git a/tools/rivers/bgc/write_glofas_ave.py b/tools/rivers/bgc/NWA/write_glofas_ave.py similarity index 100% rename from tools/rivers/bgc/write_glofas_ave.py rename to tools/rivers/bgc/NWA/write_glofas_ave.py diff --git a/xmls/NEP10/CEFI_NEP_cobalt.xml b/xmls/NEP10/CEFI_NEP_cobalt.xml index 7d6e97c84..4a1a2618f 100644 --- a/xmls/NEP10/CEFI_NEP_cobalt.xml +++ b/xmls/NEP10/CEFI_NEP_cobalt.xml @@ -21,6 +21,9 @@ C6: fremake -x CEFI_NEP_cobalt.xml -p ncrc6.intel23 -t repro MOM6_SIS2_GENERIC_4P_compile_symm frerun -x CEFI_NEP_cobalt.xml -p ncrc6.intel23 -t repro CEFI_NEP_COBALT_V1 +pp (on PPAN): +frepp -t 2019 -c ocean_cobalt_tracers_month_z -d /archive/e1n/fre/cefi/NEP/2024_08/NEP10k_082024_clean_spinup/gfdl.ncrc5-intel22-repro/history/ -x CEFI_NEP_cobalt.xml -p gfdl.ncrc6-intel23 -T repro CEFI_NEP_COBALT_V1 + Regression test frerun -x CEFI_NEP_cobalt.xml -p ncrc6.intel23 -q debug -r test -t repro CEFI_NEP_COBALT_V1 @@ -63,7 +66,7 @@ frecheck -v -r restart -p ncrc6.intel23 -t repro -x CEFI_NEP_cobalt.xml CEFI_NE - + @@ -599,12 +602,7 @@ cat > $work/INPUT/MOM_override << MOM_OVERRIDE_EOF #override PHA_MLD_DRHO = 0.03 #override HREF_FOR_MLD = 5.0 #override DIAG_MLD_DENSITY_DIFF = 0.03 -#override OBC_REMAPPING_USE_OM4_SUBCELLS = True -#override REMAPPING_USE_OM4_SUBCELLS = True -#override DIAG_REMAPPING_USE_OM4_SUBCELLS = True #override VISC_REM_BUG = True -#override FRICTWORK_BUG = True -#override INTWAVE_REMAPPING_USE_OM4_SUBCELLS = True #override GENERIC_TRACER_IC_FILE = "NEP_COBALT_spinup_2003.nc" MOM_OVERRIDE_EOF @@ -704,6 +702,7 @@ COBALT_INPUT_EOF ]]> + @@ -986,8 +985,10 @@ endif - - + + + + diff --git a/xmls/NEP10/MOM_inputs/2024_08/MOM_input b/xmls/NEP10/MOM_inputs/2024_08/MOM_input index 432c9030f..417b356cd 100644 --- a/xmls/NEP10/MOM_inputs/2024_08/MOM_input +++ b/xmls/NEP10/MOM_inputs/2024_08/MOM_input @@ -189,10 +189,6 @@ OBC_REMAPPING_USE_OM4_SUBCELLS = False ! [Boolean] default = True ! If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. ! See REMAPPING_USE_OM4_SUBCELLS for more details. We recommend setting this ! option to false. -OBC_REMAPPING_USE_OM4_SUBCELLS = False ! [Boolean] default = True - ! If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. - ! See REMAPPING_USE_OM4_SUBCELLS for more details. We recommend setting this - ! option to false. MASKING_DEPTH = 1.0 ! [m] default = -9999.0 ! The depth below which to mask points as land points, for which all ! fluxes are zeroed out. MASKING_DEPTH is ignored if negative.