| [24717] | 1 | %Test Name: Earth_Antarctica_GIA
|
|---|
| 2 |
|
|---|
| 3 | %Data paths: {{{
|
|---|
| 4 | modeldatapath='/Users/larour/ModelData';
|
|---|
| 5 | shppath='/Users/larour/issm-jpl/proj-group/qgis/Slr';
|
|---|
| 6 | gshhsshapefile=[modeldatapath '/Gshhg/Shp/GSHHS_shp/c/GSHHS_c_L1.shp'];
|
|---|
| 7 | %}}}
|
|---|
| 8 |
|
|---|
| 9 | %create sealevel model to hold our information:
|
|---|
| 10 | sl=sealevelmodel();
|
|---|
| 11 |
|
|---|
| 12 | %Create basins using boundaries from shapefile: %{{{
|
|---|
| 13 | sl.addbasin(basin('continent','southamerica','name','southamerica','epsg',5530,'boundaries',{... %{{{
|
|---|
| 14 | boundary('shppath',shppath,'shpfilename','PanamaWest','epsg',4326),...
|
|---|
| 15 | boundary('shppath',shppath,'shpfilename','SouthAmericaNorthWest','epsg',4326,'orientation','reverse'),...
|
|---|
| 16 | boundary('shppath',shppath,'shpfilename','NorthAmericaSouthWest','epsg',4326),...
|
|---|
| 17 | boundary('shppath',shppath,'shpfilename','SouthAmericaWest','epsg',4326),...
|
|---|
| 18 | boundary('shppath',shppath,'shpfilename','MarieByrd1Summit','epsg',3031),...
|
|---|
| 19 | boundary('shppath',shppath,'shpfilename','Ellsworth2','epsg',3031,'orientation','reverse'),...
|
|---|
| 20 | boundary('shppath',shppath,'shpfilename','Ellsworth1Summit','epsg',3031),...
|
|---|
| 21 | boundary('shppath',shppath,'shpfilename','Bellingshausen2','epsg',3031,'orientation','reverse'),...
|
|---|
| 22 | boundary('shppath',shppath,'shpfilename','Bellingshausen1Summit','epsg',3031),...
|
|---|
| 23 | boundary('shppath',shppath,'shpfilename','Wilkins1','epsg',3031),...
|
|---|
| 24 | boundary('shppath',shppath,'shpfilename','Wilkins2Summit','epsg',3031),...
|
|---|
| 25 | boundary('shppath',shppath,'shpfilename','Peninsula','epsg',3031),...
|
|---|
| 26 | boundary('shppath',shppath,'shpfilename','Palmer1Summit','epsg',3031),...
|
|---|
| 27 | boundary('shppath',shppath,'shpfilename','Palmer1','epsg',3031),...
|
|---|
| 28 | boundary('shppath',shppath,'shpfilename','Ronne2Summit','epsg',3031),...
|
|---|
| 29 | boundary('shppath',shppath,'shpfilename','Ronne2b','epsg',3031),...
|
|---|
| 30 | boundary('shppath',shppath,'shpfilename','SlessorSummit','epsg',3031),...
|
|---|
| 31 | boundary('shppath',shppath,'shpfilename','Slessor1','epsg',3031,'orientation','reverse'),...
|
|---|
| 32 | boundary('shppath',shppath,'shpfilename','BruntSummit','epsg',3031),...
|
|---|
| 33 | boundary('shppath',shppath,'shpfilename','Brunt1','epsg',3031),...
|
|---|
| 34 | boundary('shppath',shppath,'shpfilename','DronningMaud1Summit','epsg',3031),...
|
|---|
| 35 | boundary('shppath',shppath,'shpfilename','SouthAmericaEast','epsg',4326),...
|
|---|
| 36 | boundary('shppath',shppath,'shpfilename','MidAtlantic','epsg',4326),...
|
|---|
| 37 | boundary('shppath',shppath,'shpfilename','SouthAmericaNorthEast','epsg',4326,'orientation','reverse'),...
|
|---|
| 38 | boundary('shppath',shppath,'shpfilename','PanamaEast','epsg',4326),...
|
|---|
| 39 | boundary('shppath',shppath,'shpfilename','Panama','epsg',4326)}));
|
|---|
| 40 | %}}}
|
|---|
| 41 | sl.addbasin(basin('continent','northamerica','name','northamerica','epsg',3628,'boundaries',{... %{{{
|
|---|
| 42 | boundary('shppath',shppath,'shpfilename','NorthAmericaSouthWest','epsg',4326),...
|
|---|
| 43 | boundary('shppath',shppath,'shpfilename','SouthAmericaNorthWest','epsg',4326),...
|
|---|
| 44 | boundary('shppath',shppath,'shpfilename','PanamaWest','epsg',4326),...
|
|---|
| 45 | boundary('shppath',shppath,'shpfilename','Panama','epsg',4326,'orientation','reverse'),...
|
|---|
| 46 | boundary('shppath',shppath,'shpfilename','PanamaEast','epsg',4326),...
|
|---|
| 47 | boundary('shppath',shppath,'shpfilename','SouthAmericaNorthEast','epsg',4326),...
|
|---|
| 48 | boundary('shppath',shppath,'shpfilename','MidAtlantic','epsg',4326),...
|
|---|
| 49 | boundary('shppath',shppath,'shpfilename','Atlantic','epsg',4326),...
|
|---|
| 50 | boundary('shppath',shppath,'shpfilename','South2Summit','epsg',3413),...
|
|---|
| 51 | boundary('shppath',shppath,'shpfilename','South2','epsg',3413,'orientation','reverse'),...
|
|---|
| 52 | boundary('shppath',shppath,'shpfilename','South1Summit','epsg',3413),...
|
|---|
| 53 | boundary('shppath',shppath,'shpfilename','SouthWest1','epsg',3413,'orientation','reverse'),...
|
|---|
| 54 | boundary('shppath',shppath,'shpfilename','SouthWest1Summit','epsg',3413),...
|
|---|
| 55 | boundary('shppath',shppath,'shpfilename','Russell1','epsg',3413,'orientation','reverse'),...
|
|---|
| 56 | boundary('shppath',shppath,'shpfilename','Russell1Summit','epsg',3413),...
|
|---|
| 57 | boundary('shppath',shppath,'shpfilename','Jakobshavn1','epsg',3413,'orientation','reverse'),...
|
|---|
| 58 | boundary('shppath',shppath,'shpfilename','Jakobshavn1Summit','epsg',3413),...
|
|---|
| 59 | boundary('shppath',shppath,'shpfilename','Store1','epsg',3413,'orientation','reverse'),...
|
|---|
| 60 | boundary('shppath',shppath,'shpfilename','Store1Summit','epsg',3413),...
|
|---|
| 61 | boundary('shppath',shppath,'shpfilename','NorthWest1','epsg',3413,'orientation','reverse'),...
|
|---|
| 62 | boundary('shppath',shppath,'shpfilename','NorthWest1Summit','epsg',3413),...
|
|---|
| 63 | boundary('shppath',shppath,'shpfilename','Heilprin1','epsg',3413,'orientation','reverse'),...
|
|---|
| 64 | boundary('shppath',shppath,'shpfilename','Petermann2Summit','epsg',3413),...
|
|---|
| 65 | boundary('shppath',shppath,'shpfilename','Petermann1','epsg',3413,'orientation','reverse'),...
|
|---|
| 66 | boundary('shppath',shppath,'shpfilename','Petermann1Summit','epsg',3413),...
|
|---|
| 67 | boundary('shppath',shppath,'shpfilename','Ostenfeld1','epsg',3413,'orientation','reverse'),...
|
|---|
| 68 | boundary('shppath',shppath,'shpfilename','Ostenfeld1Summit','epsg',3413),...
|
|---|
| 69 | boundary('shppath',shppath,'shpfilename','Academy1','epsg',3413,'orientation','reverse'),...
|
|---|
| 70 | boundary('shppath',shppath,'shpfilename','Academy1Summit','epsg',3413),...
|
|---|
| 71 | boundary('shppath',shppath,'shpfilename','ArcticCanada','epsg',4326),...
|
|---|
| 72 | boundary('shppath',shppath,'shpfilename','NorthAmericaNorthWest','epsg',4326),...
|
|---|
| 73 | boundary('shppath',shppath,'shpfilename','NorthAmericaWest','epsg',4326)}));
|
|---|
| 74 | %}}}
|
|---|
| 75 | sl.addbasin(basin('continent','australia','name','australia','epsg',4462,'boundaries',{... %{{{
|
|---|
| 76 | boundary('shppath',shppath,'shpfilename','Ross3Summit','epsg',3031),...
|
|---|
| 77 | boundary('shppath',shppath,'shpfilename','Australia','epsg',4326,'orientation','reverse'),...
|
|---|
| 78 | boundary('shppath',shppath,'shpfilename','SouthAsiaSummit','epsg',4326),...
|
|---|
| 79 | boundary('shppath',shppath,'shpfilename','SouthAsia','epsg',4326,'orientation','reverse'),...
|
|---|
| 80 | boundary('shppath',shppath,'shpfilename','IndianOceanSummit','epsg',4326),...
|
|---|
| 81 | boundary('shppath',shppath,'shpfilename','IndianOcean','epsg',4326,'orientation','reverse'),...
|
|---|
| 82 | boundary('shppath',shppath,'shpfilename','AmerySummit','epsg',3031),...
|
|---|
| 83 | boundary('shppath',shppath,'shpfilename','QueenMary2','epsg',3031),...
|
|---|
| 84 | boundary('shppath',shppath,'shpfilename','QueenMary1Summit','epsg',3031),...
|
|---|
| 85 | boundary('shppath',shppath,'shpfilename','Wilkes3','epsg',3031,'orientation','reverse'),...
|
|---|
| 86 | boundary('shppath',shppath,'shpfilename','Wilkes2Summit','epsg',3031),...
|
|---|
| 87 | boundary('shppath',shppath,'shpfilename','Adelie2','epsg',3031,'orientation','reverse'),...
|
|---|
| 88 | boundary('shppath',shppath,'shpfilename','Adelie1Summit','epsg',3031),...
|
|---|
| 89 | boundary('shppath',shppath,'shpfilename','Victoria2','epsg',3031,'orientation','reverse'),...
|
|---|
| 90 | boundary('shppath',shppath,'shpfilename','Victoria1Summit','epsg',3031),...
|
|---|
| 91 | boundary('shppath',shppath,'shpfilename','Ross4','epsg',3031,'orientation','reverse')}));
|
|---|
| 92 | %}}}
|
|---|
| 93 | sl.addbasin(basin('continent','eurasia','name','eurasia','epsg',3416,'boundaries',{... %{{{
|
|---|
| 94 | boundary('shppath',shppath,'shpfilename','South2Summit','epsg',3413),...
|
|---|
| 95 | boundary('shppath',shppath,'shpfilename','Atlantic','epsg',4326,'orientation','reverse'),...
|
|---|
| 96 | boundary('shppath',shppath,'shpfilename','MidAtlantic','epsg',4326),...
|
|---|
| 97 | boundary('shppath',shppath,'shpfilename','SouthAmericaEast','epsg',4326,'orientation','reverse'),...
|
|---|
| 98 | boundary('shppath',shppath,'shpfilename','DronningMaud1Summit','epsg',3031),...
|
|---|
| 99 | boundary('shppath',shppath,'shpfilename','DronningMaud1','epsg',3031,'orientation','reverse'),...
|
|---|
| 100 | boundary('shppath',shppath,'shpfilename','QueenMaud2Summit','epsg',3031),...
|
|---|
| 101 | boundary('shppath',shppath,'shpfilename','QueenMaud2','epsg',3031,'orientation','reverse'),...
|
|---|
| 102 | boundary('shppath',shppath,'shpfilename','Kemp2Summit','epsg',3031),...
|
|---|
| 103 | boundary('shppath',shppath,'shpfilename','Kemp2','epsg',3031,'orientation','reverse'),...
|
|---|
| 104 | boundary('shppath',shppath,'shpfilename','KempSummit','epsg',3031),...
|
|---|
| 105 | boundary('shppath',shppath,'shpfilename','Amery5','epsg',3031),...
|
|---|
| 106 | boundary('shppath',shppath,'shpfilename','AmerySummit','epsg',3031),...
|
|---|
| 107 | boundary('shppath',shppath,'shpfilename','IndianOcean','epsg',4326),...
|
|---|
| 108 | boundary('shppath',shppath,'shpfilename','IndianOceanSummit','epsg',4326),...
|
|---|
| 109 | boundary('shppath',shppath,'shpfilename','SouthAsia','epsg',4326),...
|
|---|
| 110 | boundary('shppath',shppath,'shpfilename','SouthAsiaSummit','epsg',4326),...
|
|---|
| 111 | boundary('shppath',shppath,'shpfilename','ArcticRussia','epsg',4326,'orientation','reverse'),...
|
|---|
| 112 | boundary('shppath',shppath,'shpfilename','NorthAmericaNorthWest','epsg',4326),...
|
|---|
| 113 | boundary('shppath',shppath,'shpfilename','ArcticCanada','epsg',4326,'orientation','reverse'),...
|
|---|
| 114 | boundary('shppath',shppath,'shpfilename','Academy1Summit','epsg',3413),...
|
|---|
| 115 | boundary('shppath',shppath,'shpfilename','NEGIS1','epsg',3413,'orientation','reverse'),...
|
|---|
| 116 | boundary('shppath',shppath,'shpfilename','NEGIS1Summit','epsg',3413),...
|
|---|
| 117 | boundary('shppath',shppath,'shpfilename','Bistrup2','epsg',3413,'orientation','reverse'),...
|
|---|
| 118 | boundary('shppath',shppath,'shpfilename','Bistrup1Summit','epsg',3413),...
|
|---|
| 119 | boundary('shppath',shppath,'shpfilename','Daugaard-Jensen5','epsg',3413,'orientation','reverse'),...
|
|---|
| 120 | boundary('shppath',shppath,'shpfilename','Daugaard-Jensen-Summit4','epsg',3413),...
|
|---|
| 121 | boundary('shppath',shppath,'shpfilename','Geikie2','epsg',3413,'orientation','reverse'),...
|
|---|
| 122 | boundary('shppath',shppath,'shpfilename','Geikie1Summit','epsg',3413),...
|
|---|
| 123 | boundary('shppath',shppath,'shpfilename','Kangerlussuaq3','epsg',3413,'orientation','reverse'),...
|
|---|
| 124 | boundary('shppath',shppath,'shpfilename','Kangerlussuaq2Summit','epsg',3413),...
|
|---|
| 125 | boundary('shppath',shppath,'shpfilename','Helheim3','epsg',3413,'orientation','reverse'),...
|
|---|
| 126 | boundary('shppath',shppath,'shpfilename','Helheim2Summit','epsg',3413),...
|
|---|
| 127 | boundary('shppath',shppath,'shpfilename','KogeBugt3','epsg',3413,'orientation','reverse'),...
|
|---|
| 128 | boundary('shppath',shppath,'shpfilename','KogeBugt2Summit','epsg',3413),...
|
|---|
| 129 | boundary('shppath',shppath,'shpfilename','SouthEast3','epsg',3413,'orientation','reverse'),...
|
|---|
| 130 | boundary('shppath',shppath,'shpfilename','SouthEast2Summit','epsg',3413),...
|
|---|
| 131 | boundary('shppath',shppath,'shpfilename','South3','epsg',3413,'orientation','reverse')}));
|
|---|
| 132 | %}}}
|
|---|
| 133 | sl.addbasin(basin('continent','pacific','name','pacific','epsg',3031,'boundaries',{... %{{{
|
|---|
| 134 | boundary('shppath',shppath,'shpfilename','Ross3Summit','epsg',3031),...
|
|---|
| 135 | boundary('shppath',shppath,'shpfilename','MarieByrd2','epsg',3031,'orientation','reverse'),...
|
|---|
| 136 | boundary('shppath',shppath,'shpfilename','MarieByrd1Summit','epsg',3031),...
|
|---|
| 137 | boundary('shppath',shppath,'shpfilename','SouthAmericaWest','epsg',4326,'orientation','reverse'),...
|
|---|
| 138 | boundary('shppath',shppath,'shpfilename','NorthAmericaSouthWest','epsg',4326),...
|
|---|
| 139 | boundary('shppath',shppath,'shpfilename','NorthAmericaWest','epsg',4326,'orientation','reverse'),...
|
|---|
| 140 | boundary('shppath',shppath,'shpfilename','NorthAmericaNorthWest','epsg',4326),...
|
|---|
| 141 | boundary('shppath',shppath,'shpfilename','ArcticRussia','epsg',4326),...
|
|---|
| 142 | boundary('shppath',shppath,'shpfilename','SouthAsiaSummit','epsg',4326),...
|
|---|
| 143 | boundary('shppath',shppath,'shpfilename','Australia','epsg',4326)}));
|
|---|
| 144 | %}}}
|
|---|
| 145 | sl.addbasin(basin('continent','greenland','name','greenland','epsg',3413,'boundaries',{... %{{{
|
|---|
| 146 | boundary('shppath',shppath,'shpfilename','South2Summit','epsg',3413),...
|
|---|
| 147 | boundary('shppath',shppath,'shpfilename','South2','epsg',3413,'orientation','reverse'),...
|
|---|
| 148 | boundary('shppath',shppath,'shpfilename','South1Summit','epsg',3413),...
|
|---|
| 149 | boundary('shppath',shppath,'shpfilename','SouthWest1','epsg',3413,'orientation','reverse'),...
|
|---|
| 150 | boundary('shppath',shppath,'shpfilename','SouthWest1Summit','epsg',3413),...
|
|---|
| 151 | boundary('shppath',shppath,'shpfilename','Russell1','epsg',3413,'orientation','reverse'),...
|
|---|
| 152 | boundary('shppath',shppath,'shpfilename','Russell1Summit','epsg',3413),...
|
|---|
| 153 | boundary('shppath',shppath,'shpfilename','Jakobshavn1','epsg',3413,'orientation','reverse'),...
|
|---|
| 154 | boundary('shppath',shppath,'shpfilename','Jakobshavn1Summit','epsg',3413),...
|
|---|
| 155 | boundary('shppath',shppath,'shpfilename','Store1','epsg',3413,'orientation','reverse'),...
|
|---|
| 156 | boundary('shppath',shppath,'shpfilename','Store1Summit','epsg',3413),...
|
|---|
| 157 | boundary('shppath',shppath,'shpfilename','NorthWest1','epsg',3413,'orientation','reverse'),...,
|
|---|
| 158 | boundary('shppath',shppath,'shpfilename','NorthWest1Summit','epsg',3413),...
|
|---|
| 159 | boundary('shppath',shppath,'shpfilename','Heilprin1','epsg',3413,'orientation','reverse'),...,
|
|---|
| 160 | boundary('shppath',shppath,'shpfilename','Petermann2Summit','epsg',3413),...
|
|---|
| 161 | boundary('shppath',shppath,'shpfilename','Petermann1','epsg',3413,'orientation','reverse'),...,
|
|---|
| 162 | boundary('shppath',shppath,'shpfilename','Petermann1Summit','epsg',3413),...
|
|---|
| 163 | boundary('shppath',shppath,'shpfilename','Ostenfeld1','epsg',3413,'orientation','reverse'),...,
|
|---|
| 164 | boundary('shppath',shppath,'shpfilename','Ostenfeld1Summit','epsg',3413),...
|
|---|
| 165 | boundary('shppath',shppath,'shpfilename','Academy1','epsg',3413,'orientation','reverse'),...,
|
|---|
| 166 | boundary('shppath',shppath,'shpfilename','Academy1Summit','epsg',3413),...
|
|---|
| 167 | boundary('shppath',shppath,'shpfilename','NEGIS1','epsg',3413,'orientation','reverse'),...,
|
|---|
| 168 | boundary('shppath',shppath,'shpfilename','NEGIS1Summit','epsg',3413),...
|
|---|
| 169 | boundary('shppath',shppath,'shpfilename','Bistrup2','epsg',3413,'orientation','reverse'),...,
|
|---|
| 170 | boundary('shppath',shppath,'shpfilename','Bistrup1Summit','epsg',3413),...
|
|---|
| 171 | boundary('shppath',shppath,'shpfilename','Daugaard-Jensen5','epsg',3413,'orientation','reverse'),...,
|
|---|
| 172 | boundary('shppath',shppath,'shpfilename','Daugaard-Jensen-Summit4','epsg',3413),...
|
|---|
| 173 | boundary('shppath',shppath,'shpfilename','Geikie2','epsg',3413,'orientation','reverse'),...,
|
|---|
| 174 | boundary('shppath',shppath,'shpfilename','Geikie1Summit','epsg',3413),...
|
|---|
| 175 | boundary('shppath',shppath,'shpfilename','Kangerlussuaq3','epsg',3413,'orientation','reverse'),...,
|
|---|
| 176 | boundary('shppath',shppath,'shpfilename','Kangerlussuaq2Summit','epsg',3413),...
|
|---|
| 177 | boundary('shppath',shppath,'shpfilename','Helheim3','epsg',3413,'orientation','reverse'),...,
|
|---|
| 178 | boundary('shppath',shppath,'shpfilename','Helheim2Summit','epsg',3413),...
|
|---|
| 179 | boundary('shppath',shppath,'shpfilename','KogeBugt3','epsg',3413,'orientation','reverse'),...,
|
|---|
| 180 | boundary('shppath',shppath,'shpfilename','KogeBugt2Summit','epsg',3413),...
|
|---|
| 181 | boundary('shppath',shppath,'shpfilename','SouthEast3','epsg',3413,'orientation','reverse'),...,
|
|---|
| 182 | boundary('shppath',shppath,'shpfilename','SouthEast2Summit','epsg',3413),...
|
|---|
| 183 | boundary('shppath',shppath,'shpfilename','South3','epsg',3413,'orientation','reverse')}));
|
|---|
| 184 | %}}}
|
|---|
| 185 | sl.addbasin(basin('continent','antarctica','name','antarctica','epsg',3031,'boundaries',{... %{{{
|
|---|
| 186 | boundary('shppath',shppath,'shpfilename','DronningMaud1Summit','epsg',3031),...
|
|---|
| 187 | boundary('shppath',shppath,'shpfilename','DronningMaud1','epsg',3031,'orientation','reverse'),...
|
|---|
| 188 | boundary('shppath',shppath,'shpfilename','QueenMaud2Summit','epsg',3031),...
|
|---|
| 189 | boundary('shppath',shppath,'shpfilename','QueenMaud2','epsg',3031,'orientation','reverse'),...
|
|---|
| 190 | boundary('shppath',shppath,'shpfilename','Kemp2Summit','epsg',3031),...
|
|---|
| 191 | boundary('shppath',shppath,'shpfilename','Kemp2','epsg',3031,'orientation','reverse'),...
|
|---|
| 192 | boundary('shppath',shppath,'shpfilename','KempSummit','epsg',3031),...
|
|---|
| 193 | boundary('shppath',shppath,'shpfilename','Amery5','epsg',3031),...
|
|---|
| 194 | boundary('shppath',shppath,'shpfilename','AmerySummit','epsg',3031),...
|
|---|
| 195 | boundary('shppath',shppath,'shpfilename','QueenMary2','epsg',3031),...
|
|---|
| 196 | boundary('shppath',shppath,'shpfilename','QueenMary1Summit','epsg',3031),...
|
|---|
| 197 | boundary('shppath',shppath,'shpfilename','Wilkes3','epsg',3031,'orientation','reverse'),...
|
|---|
| 198 | boundary('shppath',shppath,'shpfilename','Wilkes2Summit','epsg',3031),...
|
|---|
| 199 | boundary('shppath',shppath,'shpfilename','Adelie2','epsg',3031,'orientation','reverse'),...
|
|---|
| 200 | boundary('shppath',shppath,'shpfilename','Adelie1Summit','epsg',3031),...
|
|---|
| 201 | boundary('shppath',shppath,'shpfilename','Victoria2','epsg',3031,'orientation','reverse'),...
|
|---|
| 202 | boundary('shppath',shppath,'shpfilename','Victoria1Summit','epsg',3031),...
|
|---|
| 203 | boundary('shppath',shppath,'shpfilename','Ross4','epsg',3031,'orientation','reverse'),...
|
|---|
| 204 | boundary('shppath',shppath,'shpfilename','Ross3Summit','epsg',3031),...
|
|---|
| 205 | boundary('shppath',shppath,'shpfilename','MarieByrd2','epsg',3031,'orientation','reverse'),...
|
|---|
| 206 | boundary('shppath',shppath,'shpfilename','MarieByrd1Summit','epsg',3031),...
|
|---|
| 207 | boundary('shppath',shppath,'shpfilename','Ellsworth2','epsg',3031,'orientation','reverse'),...
|
|---|
| 208 | boundary('shppath',shppath,'shpfilename','Ellsworth1Summit','epsg',3031),...
|
|---|
| 209 | boundary('shppath',shppath,'shpfilename','Bellingshausen2','epsg',3031,'orientation','reverse'),...
|
|---|
| 210 | boundary('shppath',shppath,'shpfilename','Bellingshausen1Summit','epsg',3031),...
|
|---|
| 211 | boundary('shppath',shppath,'shpfilename','Wilkins1','epsg',3031),...
|
|---|
| 212 | boundary('shppath',shppath,'shpfilename','Wilkins2Summit','epsg',3031),...
|
|---|
| 213 | boundary('shppath',shppath,'shpfilename','Peninsula','epsg',3031),...
|
|---|
| 214 | boundary('shppath',shppath,'shpfilename','Palmer1Summit','epsg',3031),...
|
|---|
| 215 | boundary('shppath',shppath,'shpfilename','Palmer1','epsg',3031),...
|
|---|
| 216 | boundary('shppath',shppath,'shpfilename','Ronne2Summit','epsg',3031),...
|
|---|
| 217 | boundary('shppath',shppath,'shpfilename','Ronne2b','epsg',3031),...
|
|---|
| 218 | boundary('shppath',shppath,'shpfilename','SlessorSummit','epsg',3031),...
|
|---|
| 219 | boundary('shppath',shppath,'shpfilename','Slessor1','epsg',3031,'orientation','reverse'),...
|
|---|
| 220 | boundary('shppath',shppath,'shpfilename','BruntSummit','epsg',3031),...
|
|---|
| 221 | boundary('shppath',shppath,'shpfilename','Brunt1','epsg',3031)}));
|
|---|
| 222 | %}}}
|
|---|
| 223 | %}}}
|
|---|
| 224 | %Go through basins and mesh: %{{{
|
|---|
| 225 | %meshing parameters: {{{
|
|---|
| 226 | hmin=500; hmax=2000; hmin=hmin*1000; hmax=hmax*1000;
|
|---|
| 227 | tolerance=100; %tolerance of 100m on Earth position when mergin 3d meshes
|
|---|
| 228 | threshold=1;
|
|---|
| 229 | defaultoptions={'KeepVertices',0,'MaxCornerAngle',0.0000000001,'NoBoundaryRefinment',1};
|
|---|
| 230 | alreadyloaded=0;
|
|---|
| 231 | %}}}
|
|---|
| 232 | for i=sl.basinindx('basin','all'),
|
|---|
| 233 |
|
|---|
| 234 | bas=sl.basins{i};
|
|---|
| 235 | disp(sprintf('Meshing basin %s\n',bas.name));
|
|---|
| 236 |
|
|---|
| 237 | %recover basin domain:
|
|---|
| 238 | domain=bas.contour();
|
|---|
| 239 |
|
|---|
| 240 | %recover coastline inside basin, using GSHHS_c_L1. It's a lat,long file, hence epsg 4326
|
|---|
| 241 | coastline=bas.shapefilecrop('shapefile',gshhsshapefile,'epsgshapefile',4326,'threshold',threshold);
|
|---|
| 242 |
|
|---|
| 243 | %mesh:
|
|---|
| 244 | md=bamg(model,'domain',domain,'subdomains',coastline,'hmin',hmin,'hmax',hmax,defaultoptions{:});
|
|---|
| 245 | plotmodel(md,'data','mesh');pause(1);
|
|---|
| 246 | error;
|
|---|
| 247 |
|
|---|
| 248 | %miscellaneous:
|
|---|
| 249 | md.mesh.epsg=bas.epsg; md.miscellaneous.name=bas.name;
|
|---|
| 250 |
|
|---|
| 251 | %recover mask where we have land:
|
|---|
| 252 | md.private.bamg.landmask=double(md.private.bamg.mesh.Triangles(:,4)>=1);
|
|---|
| 253 |
|
|---|
| 254 | %vertex connectivity:
|
|---|
| 255 | md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
|
|---|
| 256 |
|
|---|
| 257 | %add model to sl icecaps:
|
|---|
| 258 | sl.addicecap(md);
|
|---|
| 259 | end
|
|---|
| 260 | %}}}
|
|---|
| 261 | error;
|
|---|
| 262 | %Parameterize ice sheets : {{{
|
|---|
| 263 |
|
|---|
| 264 | for ind=sl.basinindx('continent',{'greenland','antarctica'}),
|
|---|
| 265 | disp(sprintf('Parameterizing basin %s\n', sl.icecaps{ind}.miscellaneous.name));
|
|---|
| 266 |
|
|---|
| 267 | md=sl.icecaps{ind}; bas=sl.basins{ind};
|
|---|
| 268 |
|
|---|
| 269 | %basin/continent specific code: {{{
|
|---|
| 270 | if bas.iscontinentany('antarctica'),
|
|---|
| 271 | openstreet=openstreetantarctica;
|
|---|
| 272 | domainoutline=domainoutlinea;
|
|---|
| 273 | end;
|
|---|
| 274 | if bas.iscontinentany('greenland'),
|
|---|
| 275 | openstreet=openstreetgreenland;
|
|---|
| 276 | domainoutline=domainoutlineg;
|
|---|
| 277 | end
|
|---|
| 278 | %}}}
|
|---|
| 279 | %masks : %{{{
|
|---|
| 280 | %new type of mask:
|
|---|
| 281 | md.mask=maskpsl;
|
|---|
| 282 |
|
|---|
| 283 | %land and ocean from open street:
|
|---|
| 284 | flags=1-(~ContourToNodes(md.mesh.x,md.mesh.y,openstreet,0));
|
|---|
| 285 | pos=find(flags==0); flags(pos)=-1;
|
|---|
| 286 | md.mask.land_levelset=flags;
|
|---|
| 287 | md.mask.ocean_levelset=-flags;
|
|---|
| 288 |
|
|---|
| 289 | %ice levelset from domain outlines:
|
|---|
| 290 | flags= -(~ContourToNodes(md.mesh.x,md.mesh.y,[shppath '/' domainoutline],1));
|
|---|
| 291 | mds=extract(md,[shppath '/' domainoutline]);
|
|---|
| 292 | flags(mds.mesh.extractedvertices)=~mds.mesh.vertexonboundary;
|
|---|
| 293 | pos=find(flags==0 & md.mesh.vertexonboundary); flags(pos)=1;
|
|---|
| 294 | md.mask.ice_levelset=-flags;
|
|---|
| 295 |
|
|---|
| 296 | %on ice front, we are not on ocean;
|
|---|
| 297 | pos=find(md.mask.ice_levelset==0); md.mask.ocean_levelset(pos)=-1; md.mask.land_levelset(pos)=1;
|
|---|
| 298 |
|
|---|
| 299 | %now, glaciers:
|
|---|
| 300 | mesh_glacier_levelset=md.private.bamg.landmask; %see meshing for why
|
|---|
| 301 | md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1);
|
|---|
| 302 | pos=find(mesh_glacier_levelset);
|
|---|
| 303 | md.mask.glacier_levelset(md.mesh.elements(pos,:))=1;
|
|---|
| 304 |
|
|---|
| 305 | pos=find(md.mask.glacier_levelset);
|
|---|
| 306 | if ~isempty(pos),
|
|---|
| 307 | mds=extract(md,md.mask.glacier_levelset);
|
|---|
| 308 | md.mask.ice_levelset(mds.mesh.extractedvertices)=-(~mds.mesh.vertexonboundary);
|
|---|
| 309 | end
|
|---|
| 310 | %}}}
|
|---|
| 311 | %initial grounding line position: {{{
|
|---|
| 312 | if bas.iscontinentany('antarctica'), % {{{
|
|---|
| 313 |
|
|---|
| 314 | %figure out initial grounding line position from the bedmap2 dataset:
|
|---|
| 315 | land_type=interpBedmap2(md.mesh.x,md.mesh.y,'icemask_grounded_and_shelves');
|
|---|
| 316 |
|
|---|
| 317 | pos=find(isnan(land_type));
|
|---|
| 318 | pos2=find(~isnan(land_type));
|
|---|
| 319 | for i=1:length(pos),
|
|---|
| 320 | in=find_point(md.mesh.x(pos2),md.mesh.y(pos2),md.mesh.x(pos(i)),md.mesh.y(pos(i)));
|
|---|
| 321 | land_type(pos(i))=land_type(pos2(in));
|
|---|
| 322 | end
|
|---|
| 323 |
|
|---|
| 324 | gridoniceshelf=zeros(md.mesh.numberofvertices,1);
|
|---|
| 325 | gridoniceshelf(land_type>0.9)=land_type(land_type>0.9);
|
|---|
| 326 | gridoniceshelf(gridoniceshelf>0)=-1;
|
|---|
| 327 | gridoniceshelf(gridoniceshelf~=-1)=1;
|
|---|
| 328 |
|
|---|
| 329 | md.mask.groundedice_levelset=gridoniceshelf;
|
|---|
| 330 |
|
|---|
| 331 | %correction to land and ocean levelset: ice shelf is not on land!
|
|---|
| 332 | pos=find(md.mask.ice_levelset<=0 & md.mask.groundedice_levelset<=0);
|
|---|
| 333 | md.mask.ocean_levelset(pos)=1;
|
|---|
| 334 | md.mask.land_levelset(pos)=-1;
|
|---|
| 335 | %}}}
|
|---|
| 336 | elseif bas.iscontinentany('greenland'), % {{{
|
|---|
| 337 |
|
|---|
| 338 | mask=interpBedmachine(md.mesh.x,md.mesh.y,'mask','greenland');
|
|---|
| 339 |
|
|---|
| 340 | gridoniceshelf=ones(md.mesh.numberofvertices,1);
|
|---|
| 341 | pos=find(md.mask.ice_levelset>0); groundedice_levelset(pos)=-1;
|
|---|
| 342 | pos=find(mask<=0); gridoniceshelf(pos)=0;
|
|---|
| 343 |
|
|---|
| 344 | md.mask.groundedice_levelset=gridoniceshelf;
|
|---|
| 345 | %}}}
|
|---|
| 346 | end
|
|---|
| 347 | %}}}
|
|---|
| 348 | %latlong: % {{{
|
|---|
| 349 | [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326');
|
|---|
| 350 | %}}}
|
|---|
| 351 | %flowequation: % {{{
|
|---|
| 352 | md=setflowequation(md,'SSA','all');
|
|---|
| 353 | %}}}
|
|---|
| 354 | %geometry: {{{
|
|---|
| 355 | if bas.iscontinentany('antarctica'),% {{{
|
|---|
| 356 | di=md.materials.rho_ice/md.materials.rho_water;
|
|---|
| 357 |
|
|---|
| 358 | disp(' reading dem');
|
|---|
| 359 | surface=interpBedmachine(md.mesh.x,md.mesh.y,'surface','antarctica');
|
|---|
| 360 | surface=griddata(md.mesh.x(~isnan(surface)),md.mesh.y(~isnan(surface)),surface(~isnan(surface)),md.mesh.x,md.mesh.y,'nearest');
|
|---|
| 361 | surfaceold=surface;
|
|---|
| 362 |
|
|---|
| 363 | disp(' reading firn layer');
|
|---|
| 364 | load(firnpath)
|
|---|
| 365 | firn_layer=InterpFromGridToMesh(x1,y1,firn,md.mesh.x,md.mesh.y,0);
|
|---|
| 366 |
|
|---|
| 367 | disp(' reading bedrock');
|
|---|
| 368 | base=interpBedmachine(md.mesh.x,md.mesh.y,'bed','antarctica');
|
|---|
| 369 | bedmap=interpBedmap2(md.mesh.x,md.mesh.y,'bed');
|
|---|
| 370 | base(base<-8000 | isnan(base))=bedmap(base<-8000 | isnan(base));
|
|---|
| 371 |
|
|---|
| 372 | %thickness:
|
|---|
| 373 | thickness=surface-base;
|
|---|
| 374 |
|
|---|
| 375 | %hydrostatic equilibrium on ice shelves:
|
|---|
| 376 | pos=find(md.mask.groundedice_levelset<0); thickness(pos)=surface(pos)/(1-di)-firn_layer(pos);
|
|---|
| 377 | surface(pos)=(1-di)*thickness(pos);
|
|---|
| 378 | base(pos)=-di*thickness(pos);
|
|---|
| 379 |
|
|---|
| 380 | %firn layer over grounded ice:
|
|---|
| 381 | pos=find(md.mask.groundedice_levelset>=0); thickness(pos)=thickness(pos)-firn_layer(pos);
|
|---|
| 382 |
|
|---|
| 383 | %check again:
|
|---|
| 384 | pos=find(thickness<1); thickness(pos)=1;
|
|---|
| 385 |
|
|---|
| 386 | %reset:
|
|---|
| 387 | surface=base+thickness;
|
|---|
| 388 |
|
|---|
| 389 | %make sure we are at hydrostatic at the grounding line:
|
|---|
| 390 | bed=base;bed(md.mask.groundedice_levelset<0)=base(md.mask.groundedice_levelset<0)-2000;
|
|---|
| 391 |
|
|---|
| 392 | %set:
|
|---|
| 393 | md.geometry.surface=surface;
|
|---|
| 394 | md.geometry.bed=bed;
|
|---|
| 395 | md.geometry.base=base;
|
|---|
| 396 | md.geometry.thickness=thickness;
|
|---|
| 397 | %}}}
|
|---|
| 398 | elseif bas.iscontinentany('greenland'),% {{{
|
|---|
| 399 | mask=interpBedmachine(md.mesh.x,md.mesh.y,'mask','greenland');
|
|---|
| 400 | surface=interpBedmachine(md.mesh.x,md.mesh.y,'surface','greenland');
|
|---|
| 401 | thickness=interpBedmachine(md.mesh.x,md.mesh.y,'thickness','greenland');
|
|---|
| 402 | base=surface-thickness;
|
|---|
| 403 | di=md.materials.rho_ice/md.materials.rho_water;
|
|---|
| 404 |
|
|---|
| 405 | pos=find(isnan(base) | isnan(surface));
|
|---|
| 406 | thickness(pos)=1;
|
|---|
| 407 | surface(pos)=(1-di)*thickness(pos);
|
|---|
| 408 | base(pos)=-di*thickness(pos);
|
|---|
| 409 |
|
|---|
| 410 | %places where thickness is 0:
|
|---|
| 411 | pos=find(thickness==0);
|
|---|
| 412 | thickness(pos)=.1; surface(pos)=base(pos)+thickness(pos);
|
|---|
| 413 |
|
|---|
| 414 | %adjust bed and base:
|
|---|
| 415 | bed=base;
|
|---|
| 416 | pos=find(md.mask.groundedice_levelset>=0); bed(pos)=base(pos);
|
|---|
| 417 | pos=find(md.mask.groundedice_levelset<0); bed(pos)=base(pos)-100;
|
|---|
| 418 |
|
|---|
| 419 | md.geometry.surface=surface;
|
|---|
| 420 | md.geometry.bed=bed;
|
|---|
| 421 | md.geometry.base=bed;
|
|---|
| 422 | md.geometry.thickness=thickness;
|
|---|
| 423 |
|
|---|
| 424 | end %}}}
|
|---|
| 425 | % }}}
|
|---|
| 426 | %velocities: %{{{
|
|---|
| 427 | if bas.iscontinentany('antarctica'),% {{{
|
|---|
| 428 |
|
|---|
| 429 | velnc =[jplsvn '/database/antarctica/velocity_ref_v28Apr2011bamber_900m_median.nc'];
|
|---|
| 430 |
|
|---|
| 431 | %read and process netcdf
|
|---|
| 432 | xmin = double(ncreadatt(velnc,'/','xmin'));
|
|---|
| 433 | ymax = double(ncreadatt(velnc,'/','ymax'));
|
|---|
| 434 | nx = double(ncreadatt(velnc,'/','nx'));
|
|---|
| 435 | ny = double(ncreadatt(velnc,'/','ny'));
|
|---|
| 436 | spacing = double(ncreadatt(velnc,'/','spacing'));
|
|---|
| 437 | vx = double(ncread(velnc,'vx'));
|
|---|
| 438 | vy = double(ncread(velnc,'vy'));
|
|---|
| 439 | x=xmin+(0:1:nx)'*spacing;
|
|---|
| 440 | y=(ymax-ny*spacing)+(0:1:ny)'*spacing;
|
|---|
| 441 | vx=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0);
|
|---|
| 442 | vy=InterpFromGridToMesh(x,y,flipud(vy'),md.mesh.x,md.mesh.y,0);
|
|---|
| 443 | vel=sqrt(vx.^2+vy.^2);
|
|---|
| 444 |
|
|---|
| 445 | %}}}
|
|---|
| 446 | elseif bas.iscontinentany('greenland'),% {{{
|
|---|
| 447 |
|
|---|
| 448 | velocities=refinementmetric(md.mesh,'greenland','greenland','nsidcvxvy');
|
|---|
| 449 | vx=velocities(:,1); vy=velocities(:,2); vel=sqrt(vx.^2+vy.^2);
|
|---|
| 450 |
|
|---|
| 451 | velocities=refinementmetric(md.mesh,'greenland','greenland','joughinvxvy');
|
|---|
| 452 | vxfar=velocities(:,1); vyfar=velocities(:,2); velfar=sqrt(vxfar.^2+vyfar.^2);
|
|---|
| 453 |
|
|---|
| 454 | %use velfar inland (from mosaic), and 2008 velocities nearcoastline:
|
|---|
| 455 | pos=find(vel==0 | isnan(vel)); vel(pos)=velfar(pos); vx(pos)=vxfar(pos);vy(pos)=vyfar(pos);
|
|---|
| 456 |
|
|---|
| 457 | %set water velocity to 0:
|
|---|
| 458 | pos=find(isnan(vel)); vel(pos)=0; vx(pos)=0; vy(pos)=0;
|
|---|
| 459 |
|
|---|
| 460 | end % }}}
|
|---|
| 461 |
|
|---|
| 462 | md.inversion.vx_obs=vx;
|
|---|
| 463 | md.inversion.vy_obs=vy;
|
|---|
| 464 | md.inversion.vel_obs=vel;
|
|---|
| 465 | md.initialization.vx=md.inversion.vx_obs;
|
|---|
| 466 | md.initialization.vy=md.inversion.vy_obs;
|
|---|
| 467 | md.initialization.vz=zeros(md.mesh.numberofvertices,1);
|
|---|
| 468 | md.initialization.vel=md.inversion.vel_obs;
|
|---|
| 469 |
|
|---|
| 470 | % }}}
|
|---|
| 471 | %friction: {{{
|
|---|
| 472 | [sx,sy,s]=slope(md); sslope=averaging(md,s,0);
|
|---|
| 473 |
|
|---|
| 474 | pos=find(md.inversion.vel_obs==0);
|
|---|
| 475 | vel=max(md.inversion.vel_obs,0.1); vel(pos)=1;
|
|---|
| 476 |
|
|---|
| 477 | md.friction.coefficient=sqrt(md.materials.rho_ice*md.geometry.thickness.*(sslope)./((md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.base).*vel/md.constants.yts));
|
|---|
| 478 | md.friction.coefficient=min(md.friction.coefficient,200);
|
|---|
| 479 |
|
|---|
| 480 | min_drag_coeff=35;
|
|---|
| 481 | min_drag_coeff_outlets=5;
|
|---|
| 482 | threshhold_drag_coeff=70;
|
|---|
| 483 | background_drag_coeff=200;
|
|---|
| 484 |
|
|---|
| 485 | maxvel=30;
|
|---|
| 486 | pos=find(md.friction.coefficient<min_drag_coeff);
|
|---|
| 487 | md.friction.coefficient(pos)=min_drag_coeff;
|
|---|
| 488 |
|
|---|
| 489 | %Take care of iceshelves: no drag md.drag
|
|---|
| 490 | pos=find(md.mask.groundedice_levelset<0);
|
|---|
| 491 | md.friction.coefficient(pos)=0;
|
|---|
| 492 | md.friction.p=ones(md.mesh.numberofelements,1);
|
|---|
| 493 | md.friction.q=ones(md.mesh.numberofelements,1);
|
|---|
| 494 | % }}}
|
|---|
| 495 | %Temperatures and surface mass balance: {{{
|
|---|
| 496 | if bas.iscontinentany('antarctica'),% {{{
|
|---|
| 497 | %smb {{{
|
|---|
| 498 | smbpath = [modeldatapath2 '/RACMO2Antarctica/'];
|
|---|
| 499 |
|
|---|
| 500 | Years = 1979:2010; %we are going to need 1980 -> 1990 as average, and perturbation from 2005 -> 2010
|
|---|
| 501 |
|
|---|
| 502 | ncdata = sprintf('%s/%s-%i%s',smbpath,'smb.KNMI',1979,'.ANT3K27.DRIFTALB.nc');
|
|---|
| 503 | lat = ncread(ncdata,'lat'); lon = ncread(ncdata,'lon');
|
|---|
| 504 | lat = double(lat); lon = double(lon);
|
|---|
| 505 | lat = lat(:); lon = lon(:);
|
|---|
| 506 | [x,y]=ll2xy(lat,lon,-1);
|
|---|
| 507 | index=BamgTriangulate(x,y);
|
|---|
| 508 | di=md.materials.rho_ice/md.materials.rho_water;
|
|---|
| 509 |
|
|---|
| 510 | smbs=zeros(length(x),numel(Years));
|
|---|
| 511 | for i=1:length(Years),
|
|---|
| 512 | year=Years(i);
|
|---|
| 513 | ncdata = sprintf('%s/%s-%i%s',smbpath,'smb.KNMI',year,'.ANT3K27.DRIFTALB.nc');
|
|---|
| 514 | smb=ncread(ncdata,'smb'); smb=sum(smb(:,:,:),3)/1000/di; smb=smb(:);
|
|---|
| 515 | smbs(:,i)=smb;
|
|---|
| 516 | end
|
|---|
| 517 | smb=InterpFromMeshToMesh2d(index,x,y,smbs,md.mesh.x,md.mesh.y);
|
|---|
| 518 |
|
|---|
| 519 | %mean from 1979 to 1990:
|
|---|
| 520 | pos=find(Years>=1979 & Years<=1990);
|
|---|
| 521 | smbmean7990=mean(smb(:,pos),2);
|
|---|
| 522 | smbmean7990=griddata(md.mesh.x(~isnan(smbmean7990)),md.mesh.y(~isnan(smbmean7990)),smbmean7990(~isnan(smbmean7990)),md.mesh.x,md.mesh.y,'nearest');
|
|---|
| 523 |
|
|---|
| 524 | %mean from 2005-2010:
|
|---|
| 525 | pos=find(Years>=2005 & Years<=2010);
|
|---|
| 526 | smbmean0510=mean(smb(:,pos),2);
|
|---|
| 527 | smbmean0510=griddata(md.mesh.x(~isnan(smbmean0510)),md.mesh.y(~isnan(smbmean0510)),smbmean0510(~isnan(smbmean0510)),md.mesh.x,md.mesh.y,'nearest');
|
|---|
| 528 |
|
|---|
| 529 | md.smb.mass_balance = smbmean0510;;
|
|---|
| 530 |
|
|---|
| 531 | pos=find(md.mask.glacier_levelset);
|
|---|
| 532 | md.smb.mass_balance(pos)=smbmean0510(pos)-smbmean7990(pos);
|
|---|
| 533 |
|
|---|
| 534 | %}}}
|
|---|
| 535 | %temperature: {{{
|
|---|
| 536 | temperaturepath=['/Users/larour/ModelData/RACMO2Antarctica/annt33.mat'];
|
|---|
| 537 | load(temperaturepath);
|
|---|
| 538 | index=delaunay(x1,y1);
|
|---|
| 539 | md.initialization.temperature=InterpFromMeshToMesh2d(index,x1(:),y1(:),annt33(:),md.mesh.x,md.mesh.y,'default',273.15);
|
|---|
| 540 | md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
|
|---|
| 541 | % }}}
|
|---|
| 542 | %}}}
|
|---|
| 543 | elseif bas.iscontinentany('greenland'),% {{{
|
|---|
| 544 | %% Annual mean for the period given by Years (smb0)
|
|---|
| 545 | Years = 1960:2014; %we are going to need 1960 -> 1990 as average, and perturbation from 2010 -> 2014.
|
|---|
| 546 |
|
|---|
| 547 | smb2 = nan(md.mesh.numberofvertices,numel(Years));
|
|---|
| 548 | st2 = nan(md.mesh.numberofvertices,numel(Years));
|
|---|
| 549 |
|
|---|
| 550 | %triangulate first:
|
|---|
| 551 | ncdata = [MARname '/MARv3.5-15km-monthly-ERA-40-1971.nc'];
|
|---|
| 552 | lat = ncread(ncdata,'LAT'); lon = ncread(ncdata,'LON');
|
|---|
| 553 | lat = double(lat); lon = double(lon);
|
|---|
| 554 | lat = lat(:); lon = lon(:);
|
|---|
| 555 | [x,y]=ll2xy(lat,lon,+1,45,70);
|
|---|
| 556 | index=BamgTriangulate(x,y);
|
|---|
| 557 |
|
|---|
| 558 | smb1s=zeros(length(x),numel(Years));
|
|---|
| 559 | st1s=zeros(length(x),numel(Years));
|
|---|
| 560 |
|
|---|
| 561 | for ii = 1:numel(Years)
|
|---|
| 562 | if Years(ii)<=1978
|
|---|
| 563 | ncdata = [MARname '/MARv3.5-15km-monthly-ERA-40-' num2str(Years(ii)) '.nc'];
|
|---|
| 564 | elseif Years(ii)>1978
|
|---|
| 565 | ncdata = [MARname '/MARv3.5-15km-monthly-ERA-Interim-' num2str(Years(ii)) '.nc'];
|
|---|
| 566 | end
|
|---|
| 567 |
|
|---|
| 568 | st1 = ncread(ncdata,'STcorr'); % Surface temperature
|
|---|
| 569 | st1(st1<-0.9000e30)=nan;
|
|---|
| 570 | st1 = double(mean(st1,3)); % Annual mean surface temp
|
|---|
| 571 |
|
|---|
| 572 | st1s(:,ii)=st1(:);
|
|---|
| 573 |
|
|---|
| 574 | smb1 = ncread(ncdata,'SMBcorr'); % in mmWe/month
|
|---|
| 575 | smb1(smb1<-0.900e30)=nan;
|
|---|
| 576 | smb1 = double(sum(smb1,3)/1000); %mWe/year
|
|---|
| 577 |
|
|---|
| 578 | smb1s(:,ii)=smb1(:);
|
|---|
| 579 |
|
|---|
| 580 | end
|
|---|
| 581 |
|
|---|
| 582 | st2s=InterpFromMeshToMesh2d(index,x,y,st1s,md.mesh.x,md.mesh.y);
|
|---|
| 583 | smb2s=InterpFromMeshToMesh2d(index,x,y,smb1s,md.mesh.x,md.mesh.y);
|
|---|
| 584 |
|
|---|
| 585 | smb19601990=mean(smb2s(:,1:31),2);
|
|---|
| 586 | smb20102014=mean(smb2s(:,end-4:end),2);
|
|---|
| 587 |
|
|---|
| 588 | st20102014 = mean(st2s(:,end-4:end),2);
|
|---|
| 589 |
|
|---|
| 590 | md.initialization.temperature=273.25+st20102014;
|
|---|
| 591 | md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
|
|---|
| 592 | md.smb.mass_balance = smb20102014;
|
|---|
| 593 |
|
|---|
| 594 | pos=find(md.mask.glacier_levelset);
|
|---|
| 595 | md.smb.mass_balance(pos)=smb20102014(pos)-smb19601990(pos);
|
|---|
| 596 |
|
|---|
| 597 | end % }}}
|
|---|
| 598 | %}}}
|
|---|
| 599 | %Mechanical boundary conditions: {{{
|
|---|
| 600 | %intialize spc boundary conditions:
|
|---|
| 601 | md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
|
|---|
| 602 | md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
|
|---|
| 603 | md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
|
|---|
| 604 | md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
|
|---|
| 605 |
|
|---|
| 606 | %boundary conditions according to masks:
|
|---|
| 607 |
|
|---|
| 608 | %water (freeze velocity to 0)
|
|---|
| 609 | pos=find(md.mask.ocean_levelset>=0);
|
|---|
| 610 | md.stressbalance.spcvx(pos)=0;
|
|---|
| 611 | md.stressbalance.spcvy(pos)=0;
|
|---|
| 612 | md.stressbalance.spcvz(pos)=0;
|
|---|
| 613 |
|
|---|
| 614 | %glaciers (freeeze velocity to observed)
|
|---|
| 615 | pos=find(md.mask.glacier_levelset);
|
|---|
| 616 | md.stressbalance.spcvx(pos)=0; %flux divergence will be equated to smb perturbation.
|
|---|
| 617 | md.stressbalance.spcvy(pos)=0;
|
|---|
| 618 |
|
|---|
| 619 | md.stressbalance.spcvz(pos)=0;
|
|---|
| 620 |
|
|---|
| 621 | %land with no ice (freeeze velocity to 0)
|
|---|
| 622 | pos=find(md.mask.land_levelset>=0 & md.mask.ice_levelset>0 & ~md.mask.glacier_levelset);
|
|---|
| 623 | md.stressbalance.spcvx(pos)=0;
|
|---|
| 624 | md.stressbalance.spcvy(pos)=0;
|
|---|
| 625 | md.stressbalance.spcvz(pos)=0;
|
|---|
| 626 |
|
|---|
| 627 | %if strictly ice (and not glacier), then no boundary condition:
|
|---|
| 628 | pos=find(md.mask.ice_levelset<=0 & ~md.mask.glacier_levelset);
|
|---|
| 629 | md.stressbalance.spcvx(pos)=NaN;
|
|---|
| 630 | md.stressbalance.spcvy(pos)=NaN;
|
|---|
| 631 | md.stressbalance.spcvz(pos)=NaN;
|
|---|
| 632 |
|
|---|
| 633 | %constrain boundaries of basin to observed velocities:
|
|---|
| 634 | pos=find(md.mesh.vertexonboundary & md.mask.ice_levelset<=0);
|
|---|
| 635 | md.stressbalance.spcvx(pos)=md.inversion.vx_obs(pos);
|
|---|
| 636 | md.stressbalance.spcvy(pos)=md.inversion.vy_obs(pos);
|
|---|
| 637 | md.stressbalance.spcvz(pos)=0;
|
|---|
| 638 |
|
|---|
| 639 | %fixing quirks here and there:
|
|---|
| 640 | if strcmpi(bas.name,'Ross'),
|
|---|
| 641 | flags = ContourToNodes(mdband.mesh.x,mdband.mesh.y,[antarcticashppath '/Murdo.exp'],1);
|
|---|
| 642 | pos=find(flags);
|
|---|
| 643 |
|
|---|
| 644 | md.stressbalance.spcvx(pos)=md.inversion.vx_obs(pos);
|
|---|
| 645 | md.stressbalance.spcvy(pos)=md.inversion.vy_obs(pos);
|
|---|
| 646 | md.stressbalance.spcvz(pos)=0;
|
|---|
| 647 | end
|
|---|
| 648 |
|
|---|
| 649 | %prognostic
|
|---|
| 650 | md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
|
|---|
| 651 | % }}}
|
|---|
| 652 | %Thermal boundary conditions: {{{
|
|---|
| 653 | if bas.iscontinentany('antarctica'),% {{{
|
|---|
| 654 | searisepath =[modeldatapath '/SeaRISE/Antarctica5km_shelves_v1.0'];
|
|---|
| 655 | heatfluxpath =[searisepath '/bheatflx_fox.mat'];
|
|---|
| 656 |
|
|---|
| 657 | if strcmpi(bas.name,'ellsworth'), % {{{
|
|---|
| 658 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 659 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 660 | md.basalforcings.deepwater_melting_rate=79.3397;
|
|---|
| 661 | md.basalforcings.deepwater_elevation=-1000;
|
|---|
| 662 | md.basalforcings.upperwater_elevation=-141.59;
|
|---|
| 663 | %}}}
|
|---|
| 664 | elseif strcmpi(bas.name,'amery'), % {{{
|
|---|
| 665 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 666 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 667 | md.basalforcings.deepwater_melting_rate=8;
|
|---|
| 668 | md.basalforcings.deepwater_elevation=-1500;
|
|---|
| 669 | md.basalforcings.upperwater_elevation=-400; % }}}
|
|---|
| 670 | elseif strcmpi(bas.name,'ross'), % {{{
|
|---|
| 671 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 672 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 673 | md.basalforcings.deepwater_melting_rate=5;
|
|---|
| 674 | md.basalforcings.deepwater_elevation=-1100;
|
|---|
| 675 | md.basalforcings.upperwater_elevation=-200; % }}}
|
|---|
| 676 | elseif strcmpi(bas.name,'ronne'), % {{{
|
|---|
| 677 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 678 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 679 | md.basalforcings.deepwater_melting_rate=4;
|
|---|
| 680 | md.basalforcings.deepwater_elevation=-1800;
|
|---|
| 681 | md.basalforcings.upperwater_elevation=-400; % }}}
|
|---|
| 682 | elseif strcmpi(bas.name,'brunt'), % {{{
|
|---|
| 683 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 684 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 685 | md.basalforcings.deepwater_melting_rate=.3;
|
|---|
| 686 | md.basalforcings.deepwater_elevation=-900;
|
|---|
| 687 | md.basalforcings.upperwater_elevation=-100;
|
|---|
| 688 | % }}}
|
|---|
| 689 | elseif strcmpi(bas.name,'kemp'), % {{{
|
|---|
| 690 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 691 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 692 | md.basalforcings.deepwater_melting_rate=.2;
|
|---|
| 693 | md.basalforcings.deepwater_elevation=-800;
|
|---|
| 694 | md.basalforcings.upperwater_elevation=-100; % }}}
|
|---|
| 695 | elseif strcmpi(bas.name,'queenMaud'), % {{{
|
|---|
| 696 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 697 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 698 | md.basalforcings.deepwater_melting_rate=.5;
|
|---|
| 699 | md.basalforcings.deepwater_elevation=-650;
|
|---|
| 700 | md.basalforcings.upperwater_elevation=-180; % }}}
|
|---|
| 701 | elseif strcmpi(bas.name,'dronningMaud'), % {{{
|
|---|
| 702 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 703 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 704 | md.basalforcings.deepwater_melting_rate=3;
|
|---|
| 705 | md.basalforcings.deepwater_elevation=-800;
|
|---|
| 706 | md.basalforcings.upperwater_elevation=-300; % }}}
|
|---|
| 707 | elseif strcmpi(bas.name,'slessor'), % {{{
|
|---|
| 708 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 709 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 710 | md.basalforcings.deepwater_melting_rate=1.4;
|
|---|
| 711 | md.basalforcings.deepwater_elevation=-1200;
|
|---|
| 712 | md.basalforcings.upperwater_elevation=-500; % }}}
|
|---|
| 713 | elseif strcmpi(bas.name,'victoria'), % {{{
|
|---|
| 714 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 715 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 716 | md.basalforcings.deepwater_melting_rate=10;
|
|---|
| 717 | md.basalforcings.deepwater_elevation=-1000;
|
|---|
| 718 | md.basalforcings.upperwater_elevation=-200; % }}}
|
|---|
| 719 | elseif strcmpi(bas.name,'adelie'), % {{{
|
|---|
| 720 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 721 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 722 | md.basalforcings.deepwater_melting_rate=4;
|
|---|
| 723 | md.basalforcings.deepwater_elevation=-1000;
|
|---|
| 724 | md.basalforcings.upperwater_elevation=-200; % }}}
|
|---|
| 725 | elseif strcmpi(bas.name,'wilkes'), % {{{
|
|---|
| 726 | md.basalforcings=linerbasalforcings(md);
|
|---|
| 727 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 728 | md.basalforcings.deepwater_melting_rate=14;
|
|---|
| 729 | md.basalforcings.deepwater_elevation=-1500;
|
|---|
| 730 | md.basalforcings.upperwater_elevation=-50; % }}}
|
|---|
| 731 | elseif strcmpi(bas.name,'wilkins'), % {{{
|
|---|
| 732 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 733 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 734 | md.basalforcings.deepwater_melting_rate=6;
|
|---|
| 735 | md.basalforcings.deepwater_elevation=-600;
|
|---|
| 736 | md.basalforcings.upperwater_elevation=-50; % }}}
|
|---|
| 737 | elseif strcmpi(bas.name,'queenmary'), % {{{
|
|---|
| 738 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 739 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 740 | md.basalforcings.deepwater_melting_rate=10;
|
|---|
| 741 | md.basalforcings.deepwater_elevation=-1200;
|
|---|
| 742 | md.basalforcings.upperwater_elevation=-200; % }}}
|
|---|
| 743 | elseif strcmpi(bas.name,'mariebyrd'), % {{{
|
|---|
| 744 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 745 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 746 | md.basalforcings.deepwater_melting_rate=8;
|
|---|
| 747 | md.basalforcings.deepwater_elevation=-700;
|
|---|
| 748 | md.basalforcings.upperwater_elevation=-200; % }}}
|
|---|
| 749 | elseif strcmpi(bas.name,'bellingshausen'), % {{{
|
|---|
| 750 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 751 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 752 | md.basalforcings.deepwater_melting_rate=7;
|
|---|
| 753 | md.basalforcings.deepwater_elevation=-550;
|
|---|
| 754 | md.basalforcings.upperwater_elevation=-150; % }}}
|
|---|
| 755 | elseif strcmpi(bas.name,'palmer'), % {{{
|
|---|
| 756 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 757 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 758 | md.basalforcings.deepwater_melting_rate=.1;
|
|---|
| 759 | md.basalforcings.deepwater_elevation=-1000;
|
|---|
| 760 | md.basalforcings.upperwater_elevation=-50; % }}}
|
|---|
| 761 | elseif strcmpi(bas.name,'peninsula'), % {{{
|
|---|
| 762 | md.basalforcings=linearbasalforcings(md);
|
|---|
| 763 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 764 | md.basalforcings.deepwater_melting_rate=2;
|
|---|
| 765 | md.basalforcings.deepwater_elevation=-550;
|
|---|
| 766 | md.basalforcings.upperwater_elevation=-50; % }}}
|
|---|
| 767 | else
|
|---|
| 768 | %melt rates: {{{
|
|---|
| 769 | md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 770 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 771 | load('./Data/monthly_mean_melt_adj09.mat'); fw2=fw_flux;
|
|---|
| 772 | load('./Data/monthly_mean_melt_adj04.mat'); fw1=fw_flux;
|
|---|
| 773 | fw_flux(:,:,end+1:end+size(fw2,3))=fw2;
|
|---|
| 774 | load('./Data/grid_info_cp09.mat')
|
|---|
| 775 | [xi,yi]= ll2xy(lat,lon,-1,0,71); index=delaunay(xi,yi);
|
|---|
| 776 | mrate_mesh=InterpFromMeshToMesh2d(index,xi(:),yi(:),-1*reshape(squeeze(mean(fw_flux,3)),prod(size(xi)),1),md.mesh.x,md.mesh.y,'default',0);
|
|---|
| 777 | pos=find(md.mask.groundedice_levelset<0 & mrate_mesh~=0);
|
|---|
| 778 | md.basalforcings.floatingice_melting_rate(pos)=mrate_mesh(pos);
|
|---|
| 779 | % }}}
|
|---|
| 780 | end
|
|---|
| 781 | md.thermal.spctemperature=[md.initialization.temperature;1]; %impose observed temperature on surface
|
|---|
| 782 | load(heatfluxpath)
|
|---|
| 783 | md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,bheatflx_fox,md.mesh.x,md.mesh.y,0);
|
|---|
| 784 | % }}}
|
|---|
| 785 | elseif bas.iscontinentany('greenland'),% {{{
|
|---|
| 786 | %Initialize melt rates: {{{
|
|---|
| 787 | md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 788 | md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 789 |
|
|---|
| 790 | %Specificy melt rates and distances: {{{
|
|---|
| 791 | %Seroussi et al., 2011 - melt rates at grounding line of 79N are about 50 m/yr
|
|---|
| 792 | %Rignot and Steffen, 2008 - Petermann high rates up to 20 km out
|
|---|
| 793 | if strcmpi(bas.name,'negis'),
|
|---|
| 794 | iceshelves={'gl_zac_1996','gl_79n_1996','gl_StorBistrup96'};
|
|---|
| 795 | distances=[10,10,10];
|
|---|
| 796 | meltrates=[7,2,4];
|
|---|
| 797 | meltratesgl=[25,22,4]; %35,35
|
|---|
| 798 | elseif strcmpi(bas.name,'ostenfeld'),
|
|---|
| 799 | distances=[10,10,10,10,10];
|
|---|
| 800 | iceshelves={'gl_Ryder96','gl_Ostenfeld96','gl_Harder96','gl_Steensby96','gl_Brikkerne96'};
|
|---|
| 801 | meltrates=[7,11,4,5,4];
|
|---|
| 802 | meltratesgl=[25,26,4,5,4];
|
|---|
| 803 | elseif strcmpi(bas.name,'petermann'),
|
|---|
| 804 | distances=[20,10];
|
|---|
| 805 | iceshelves={'gl_Peter96','gl_Humb96'};
|
|---|
| 806 | meltrates=[6,4]; %5
|
|---|
| 807 | meltratesgl=[22,4]; %26
|
|---|
| 808 | elseif strcmpi(bas.name,'academy'),
|
|---|
| 809 | distances=[10];
|
|---|
| 810 | iceshelves={'gl_Hagen96'};
|
|---|
| 811 | meltrates=[4];
|
|---|
| 812 | meltratesgl=[4];
|
|---|
| 813 | elseif strcmpi(bas.name,'jakobshavn'),
|
|---|
| 814 | distances=[10];
|
|---|
| 815 | iceshelves={'gl_Jak'};
|
|---|
| 816 | meltrates=[4];
|
|---|
| 817 | meltratesgl=[109]; %Prescott et al, 2003
|
|---|
| 818 | else
|
|---|
| 819 | iceshelves={};
|
|---|
| 820 | end
|
|---|
| 821 | %}}}
|
|---|
| 822 |
|
|---|
| 823 | for i=1:length(iceshelves),
|
|---|
| 824 | iceshelf=iceshelves{i};
|
|---|
| 825 | distance=distances(i);
|
|---|
| 826 | meltrate=meltrates(i);
|
|---|
| 827 | meltrategl=meltratesgl(i);
|
|---|
| 828 |
|
|---|
| 829 | in= ContourToNodes(md.mesh.x,md.mesh.y,[shppath '/' iceshelf '_closed_polygon.shp'],1);
|
|---|
| 830 | segs=MeshProfileIntersection(md.mesh.elements,md.mesh.x,md.mesh.y,[shppath '/' iceshelf '_closed_polygon.shp']);
|
|---|
| 831 | pos=find(~isnan(segs(:,1)));
|
|---|
| 832 | segs=[segs(pos,1:2); segs(pos,3:4)];
|
|---|
| 833 |
|
|---|
| 834 | %figure out nodes within ~10 km of grounding line:
|
|---|
| 835 | v=find(in);
|
|---|
| 836 | glpt=zeros(md.mesh.numberofvertices,1);
|
|---|
| 837 | a=zeros(md.mesh.numberofvertices,1);
|
|---|
| 838 | for i=v'
|
|---|
| 839 | x=md.mesh.x(i); y=md.mesh.y(i);
|
|---|
| 840 | dist=sqrt((segs(:,1)-x).^2+(segs(:,2)-y).^2);
|
|---|
| 841 | if min(dist)<distance*1000; glpt(i)=1; else glpt(i)=0; end
|
|---|
| 842 | end
|
|---|
| 843 |
|
|---|
| 844 | md.basalforcings.floatingice_melting_rate(find(in))=meltrate;
|
|---|
| 845 | md.basalforcings.floatingice_melting_rate(find(in&glpt))=meltrategl;
|
|---|
| 846 |
|
|---|
| 847 | end
|
|---|
| 848 |
|
|---|
| 849 | %Expand to the rest of the basin, in case of retreat of the glaciers:
|
|---|
| 850 | if ~isempty(iceshelves),
|
|---|
| 851 | nomelt=md.basalforcings.floatingice_melting_rate==0;
|
|---|
| 852 | pos=find(nomelt); pos2=find(~nomelt);
|
|---|
| 853 | if ~isempty(pos2),
|
|---|
| 854 | if length(pos2)>3,
|
|---|
| 855 | md.basalforcings.floatingice_melting_rate(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),md.basalforcings.floatingice_melting_rate(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');
|
|---|
| 856 | end
|
|---|
| 857 | end
|
|---|
| 858 | end
|
|---|
| 859 | %}}}
|
|---|
| 860 | %Thermal spcs and heatflux: %{{{
|
|---|
| 861 | md.thermal.spctemperature=md.initialization.temperature; %impose observed temperature on surface
|
|---|
| 862 |
|
|---|
| 863 | load(seariseg_heatflux)
|
|---|
| 864 | [xi,yi]= ll2xy(md.mesh.lat,md.mesh.long,+1,39,71);
|
|---|
| 865 | md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,bheatflx,xi,yi,0);
|
|---|
| 866 | %}}}
|
|---|
| 867 | end % }}}
|
|---|
| 868 | %}}}
|
|---|
| 869 | %Slr: {{{
|
|---|
| 870 | if bas.iscontinentany('antarctica'),
|
|---|
| 871 |
|
|---|
| 872 | delH=textread([modeldatapath '/AdhikariSciAdvTrends/AIS_delH_trend.txt']);
|
|---|
| 873 | longAIS=delH(:,1);
|
|---|
| 874 | latAIS=delH(:,2);
|
|---|
| 875 | delHAIS=delH(:,3);
|
|---|
| 876 | index=delaunay(latAIS,longAIS);
|
|---|
| 877 |
|
|---|
| 878 | lat=md.mesh.lat;
|
|---|
| 879 | long=md.mesh.long+360;
|
|---|
| 880 | pos=find(long>360);long(pos)=long(pos)-360;
|
|---|
| 881 |
|
|---|
| 882 | delHAIS=InterpFromMesh2d(index,longAIS,latAIS,delHAIS,long,lat);
|
|---|
| 883 | northpole=find_point(md.mesh.long,md.mesh.lat,0,90); delHAIS(northpole)=0;
|
|---|
| 884 |
|
|---|
| 885 | md.slr.deltathickness=delHAIS(md.mesh.elements)*[1;1;1]/3/100;
|
|---|
| 886 | else
|
|---|
| 887 | delH=textread([modeldatapath '/AdhikariSciAdvTrends/GIS_delH_trend.txt']);
|
|---|
| 888 | longGIS=delH(:,1);
|
|---|
| 889 | latGIS=delH(:,2);
|
|---|
| 890 | delHGIS=delH(:,3);
|
|---|
| 891 |
|
|---|
| 892 | index=delaunay(latGIS,longGIS);
|
|---|
| 893 |
|
|---|
| 894 | lat=md.mesh.lat;
|
|---|
| 895 | long=md.mesh.long+360;
|
|---|
| 896 | pos=find(long>360);long(pos)=long(pos)-360;
|
|---|
| 897 |
|
|---|
| 898 | delHGIS=InterpFromMesh2d(index,longGIS,latGIS,delHGIS,long,lat);
|
|---|
| 899 | md.slr.deltathickness=delHGIS(md.mesh.elements)*[1;1;1]/3/100;
|
|---|
| 900 | end
|
|---|
| 901 | md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
|
|---|
| 902 | md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
|
|---|
| 903 | md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
|
|---|
| 904 | md.slr.Ngia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ngia,md.mesh.long,md.mesh.lat);
|
|---|
| 905 | md.slr.Ugia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ugia,md.mesh.long,md.mesh.lat);
|
|---|
| 906 | %}}}
|
|---|
| 907 | % material properties: {{{
|
|---|
| 908 | disp(' creating flow law paramter');
|
|---|
| 909 | md.materials.rheology_B=paterson(md.initialization.temperature);
|
|---|
| 910 | pos=find(md.materials.rheology_B<=0);
|
|---|
| 911 | md.materials.rheology_B(pos)=-md.materials.rheology_B(pos);
|
|---|
| 912 | md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
|
|---|
| 913 | %}}}
|
|---|
| 914 | %diverse: {{{
|
|---|
| 915 | md.miscellaneous.name=bas.name;
|
|---|
| 916 | % }}}
|
|---|
| 917 |
|
|---|
| 918 | sl.icecaps{ind}=md;
|
|---|
| 919 | end
|
|---|
| 920 | %}}}
|
|---|
| 921 | % ParameterizeContinents {{{
|
|---|
| 922 |
|
|---|
| 923 | %glacier load: {{{
|
|---|
| 924 | delH=textread([modeldatapath '/AdhikariSciAdvTrends/GLA_delH_trend_15regions.txt']);
|
|---|
| 925 | longglaciers=delH(:,1); latglaciers=delH(:,2); delHglaciers=sum(delH(:,3:end),2);
|
|---|
| 926 | indexglaciers=delaunay(latglaciers,longglaciers);
|
|---|
| 927 | %}}}
|
|---|
| 928 |
|
|---|
| 929 | for ind=sl.basinindx('continent',{'southamerica','northamerica','australia','eurasia','pacific'}),
|
|---|
| 930 | disp(sprintf('Masks for basin %s\n', sl.icecaps{ind}.miscellaneous.name));
|
|---|
| 931 | md=sl.icecaps{ind}; bas=sl.basins{ind};
|
|---|
| 932 |
|
|---|
| 933 | %recover lat,long:
|
|---|
| 934 | [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326');
|
|---|
| 935 |
|
|---|
| 936 | %interpolate glacier loads onto mesh: %{{{
|
|---|
| 937 | lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360;
|
|---|
| 938 | mdelHglaciers=InterpFromMesh2d(indexglaciers,longglaciers,latglaciers,delHglaciers,long,lat);
|
|---|
| 939 | mdelHglacierse=mdelHglaciers(md.mesh.elements)*[1;1;1]/3;
|
|---|
| 940 | pos=find(~md.private.bamg.landmask); mdelHglacierse(pos)=0;
|
|---|
| 941 | % }}}
|
|---|
| 942 | %mask: %{{{
|
|---|
| 943 | %new type of mask:
|
|---|
| 944 | md.mask=maskpsl;
|
|---|
| 945 | %Figure out mask from initial mesh: deal with land and ocean masks (land include grounded ice). %{{{
|
|---|
| 946 | %first, transform land element mask into vertex driven one.
|
|---|
| 947 | land=md.private.bamg.landmask;
|
|---|
| 948 | land_mask=-ones(md.mesh.numberofvertices,1);
|
|---|
| 949 |
|
|---|
| 950 | landels=find(land);
|
|---|
| 951 | land_mask(md.mesh.elements(landels,:))=1;
|
|---|
| 952 |
|
|---|
| 953 | %gothrough edges of each land element
|
|---|
| 954 | connectedels=md.mesh.elementconnectivity(landels,:);
|
|---|
| 955 | connectedisonocean=~land(connectedels);
|
|---|
| 956 | sumconnectedisonocean=sum(connectedisonocean,2);
|
|---|
| 957 |
|
|---|
| 958 | %figure out which land elements are connected to the ocean:
|
|---|
| 959 | landelsconocean=landels(find(sumconnectedisonocean));
|
|---|
| 960 |
|
|---|
| 961 | ind1=[md.mesh.elements(landelsconocean,1); md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3)];
|
|---|
| 962 | ind2=[md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3); md.mesh.elements(landelsconocean,1)];
|
|---|
| 963 |
|
|---|
| 964 | h=waitbar(0,'Starting mask processing');
|
|---|
| 965 |
|
|---|
| 966 | %edge ind1 and ind2:
|
|---|
| 967 | for i=1:length(ind1),
|
|---|
| 968 | if ~mod(i,100),
|
|---|
| 969 | waitbar(i/length(ind1),h,sprintf('%.2i%% along...',i/length(ind1)*100));
|
|---|
| 970 | end
|
|---|
| 971 |
|
|---|
| 972 | els1=md.mesh.vertexconnectivity(ind1(i),1: md.mesh.vertexconnectivity(ind1(i),end));
|
|---|
| 973 | els2=md.mesh.vertexconnectivity(ind2(i),1: md.mesh.vertexconnectivity(ind2(i),end));
|
|---|
| 974 | els=intersect(els1,els2);
|
|---|
| 975 |
|
|---|
| 976 | if length(find(land(els)))==1,
|
|---|
| 977 | %this edge is on the beach, 0 the edge:
|
|---|
| 978 | land_mask(ind1(i))=0;
|
|---|
| 979 | land_mask(ind2(i))=0;
|
|---|
| 980 | end
|
|---|
| 981 | end
|
|---|
| 982 | md.mask.land_levelset=land_mask;
|
|---|
| 983 | close(h);
|
|---|
| 984 | %}}}
|
|---|
| 985 | %work on ocean, glaciers and ice: {{{
|
|---|
| 986 | %ocean: opposite of land:
|
|---|
| 987 | md.mask.ocean_levelset=-md.mask.land_levelset;
|
|---|
| 988 |
|
|---|
| 989 | %initliaze:
|
|---|
| 990 | md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
|
|---|
| 991 |
|
|---|
| 992 | %%cross check that whereever we have an ice load, the mask is <0 on each vertex: if it's not th ecase, zero the ice load:
|
|---|
| 993 | pos=find(mdelHglacierse);
|
|---|
| 994 | vertices=md.mesh.elements(pos,:); %potentially having an ice load.
|
|---|
| 995 | pos=find(md.mask.land_levelset(vertices)>=0);
|
|---|
| 996 | icevertices=vertices(pos);
|
|---|
| 997 | md.mask.ice_levelset(icevertices)=-md.mask.land_levelset(icevertices);
|
|---|
| 998 |
|
|---|
| 999 | %then take care of glaciers: don't! this is determined by the land mask.
|
|---|
| 1000 | %pos=find(md.slr.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
|
|---|
| 1001 | %pos=find(mdelHglacierse); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
|
|---|
| 1002 |
|
|---|
| 1003 |
|
|---|
| 1004 | %grounded ice:
|
|---|
| 1005 | md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
|
|---|
| 1006 | pos=find(mdelHglaciers); md.mask.groundedice_levelset(pos)=1;
|
|---|
| 1007 |
|
|---|
| 1008 |
|
|---|
| 1009 | %now, glaciers:
|
|---|
| 1010 | md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1);
|
|---|
| 1011 | pos=find(mdelHglacierse); md.mask.glacier_levelset(md.mesh.elements(pos,:))=1;
|
|---|
| 1012 | % }}}
|
|---|
| 1013 | %}}}
|
|---|
| 1014 | %slr loading/calibration: {{{
|
|---|
| 1015 | %load glaciers and ice caps from GRACE:
|
|---|
| 1016 | md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
|
|---|
| 1017 | pos=find(mdelHglacierse); md.slr.deltathickness(pos)=mdelHglacierse(pos)/100;
|
|---|
| 1018 |
|
|---|
| 1019 | %initialize:
|
|---|
| 1020 | md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
|
|---|
| 1021 |
|
|---|
| 1022 | %love numbers:
|
|---|
| 1023 | md.slr.steric_rate=(1.1+.38)*ones(md.mesh.numberofvertices,1); %steric + water storage.
|
|---|
| 1024 | md.slr.ocean_area_scaling=1;
|
|---|
| 1025 | md.slr.Ngia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ngia,md.mesh.long,md.mesh.lat);
|
|---|
| 1026 | md.slr.Ugia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ugia,md.mesh.long,md.mesh.lat);
|
|---|
| 1027 | %}}}
|
|---|
| 1028 | %geometry: {{{
|
|---|
| 1029 | di=md.materials.rho_ice/md.materials.rho_water;
|
|---|
| 1030 |
|
|---|
| 1031 | md.slr.spcthickness=ones(md.mesh.numberofvertices+1,2);
|
|---|
| 1032 | %time tags: make sure they are far in the past and future.
|
|---|
| 1033 | t1=-10000; t2=+10000; md.slr.spcthickness(end,1)=-10000;
|
|---|
| 1034 | md.slr.spcthickness(end,2)=+10000;
|
|---|
| 1035 |
|
|---|
| 1036 | %impart the glacier load:
|
|---|
| 1037 | pos=find(mdelHglaciers);
|
|---|
| 1038 | md.slr.spcthickness(pos,2)=md.slr.spcthickness(pos,1)+(t2-t1)*mdelHglaciers(pos)/100*di;
|
|---|
| 1039 | %mdelHglaciers is cm/yr water equivalent
|
|---|
| 1040 |
|
|---|
| 1041 | md.geometry.thickness=ones(md.mesh.numberofvertices,1);
|
|---|
| 1042 | md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
|
|---|
| 1043 | md.geometry.base=md.geometry.surface-md.geometry.thickness;
|
|---|
| 1044 | md.geometry.bed=md.geometry.base;
|
|---|
| 1045 | % }}}
|
|---|
| 1046 | %materials: {{{
|
|---|
| 1047 | md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
|
|---|
| 1048 | md.materials.rheology_B=paterson(md.initialization.temperature);
|
|---|
| 1049 | md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
|
|---|
| 1050 | % }}}
|
|---|
| 1051 | %Solution: {{{
|
|---|
| 1052 | md.transient.issmb=0; md.transient.isstressbalance=0;
|
|---|
| 1053 | md.transient.ismasstransport=0; md.transient.isthermal=0;
|
|---|
| 1054 | md.transient.isslr=0; md.slr.loop_increment=300;
|
|---|
| 1055 | % }}}
|
|---|
| 1056 | sl.icecaps{ind}=md;
|
|---|
| 1057 | end
|
|---|
| 1058 | % }}}
|
|---|
| 1059 | %Assemble Earth in 3D {{{
|
|---|
| 1060 |
|
|---|
| 1061 | %parameters:
|
|---|
| 1062 | plotting=1;
|
|---|
| 1063 | tolerance=100;
|
|---|
| 1064 | loneedgesdetect=0;
|
|---|
| 1065 |
|
|---|
| 1066 | %create earth model by concatenating all the icecaps in 3d:
|
|---|
| 1067 | sl.caticecaps('tolerance',tolerance,'loneedgesdetect',loneedgesdetect);
|
|---|
| 1068 |
|
|---|
| 1069 | %figure out how each icecap's mesh connects to the larger earth mesh:
|
|---|
| 1070 | sl.intersections();
|
|---|
| 1071 |
|
|---|
| 1072 | %figure out connectivity:
|
|---|
| 1073 | disp('Mesh connectivity');
|
|---|
| 1074 | sl.earth.mesh.vertexconnectivity=NodeConnectivity(sl.earth.mesh.elements,sl.earth.mesh.numberofvertices);
|
|---|
| 1075 |
|
|---|
| 1076 | %areas:
|
|---|
| 1077 | disp('Mesh nodal areas');
|
|---|
| 1078 | sl.earth.mesh.area=averaging(sl.earth,GetAreas3DTria(sl.earth.mesh.elements,sl.earth.mesh.x,sl.earth.mesh.y,sl.earth.mesh.z),4);
|
|---|
| 1079 |
|
|---|
| 1080 | %transfer a list of fields from each icecap and continent back to Earth:
|
|---|
| 1081 | sl.transfer('mask.groundedice_levelset');
|
|---|
| 1082 | sl.transfer('mask.ice_levelset');
|
|---|
| 1083 | sl.transfer('mask.ocean_levelset');
|
|---|
| 1084 | sl.transfer('mask.land_levelset');
|
|---|
| 1085 | sl.transfer('mask.glacier_levelset');
|
|---|
| 1086 | sl.transfer('geometry.thickness');
|
|---|
| 1087 | sl.transfer('geometry.surface');
|
|---|
| 1088 | sl.transfer('geometry.base');
|
|---|
| 1089 | sl.transfer('geometry.bed');
|
|---|
| 1090 | sl.transfer('initialization.temperature');
|
|---|
| 1091 | sl.transfer('materials.rheology_B');
|
|---|
| 1092 | sl.transfer('mesh.lat');
|
|---|
| 1093 | sl.transfer('mesh.long');
|
|---|
| 1094 | sl.transfer('slr.deltathickness');
|
|---|
| 1095 | sl.transfer('slr.Ngia');
|
|---|
| 1096 | sl.transfer('slr.Ugia');
|
|---|
| 1097 |
|
|---|
| 1098 | %radius:
|
|---|
| 1099 | sl.earth.mesh.r=sqrt(sl.earth.mesh.x.^2+sl.earth.mesh.y.^2+sl.earth.mesh.z.^2);
|
|---|
| 1100 |
|
|---|
| 1101 | %check on the mesh transitions: {{{
|
|---|
| 1102 | if plotting,
|
|---|
| 1103 | flags=ones(sl.earth.mesh.numberofvertices,1);
|
|---|
| 1104 | for i=1:length(sl.transitions),
|
|---|
| 1105 | flags(sl.transitions{i})=i;
|
|---|
| 1106 | end
|
|---|
| 1107 | plotmodel(sl.earth,'data',flags,'shading','faceted','coastline','on','coast_color','g')
|
|---|
| 1108 | end
|
|---|
| 1109 | %}}}}
|
|---|
| 1110 |
|
|---|
| 1111 | % }}}
|
|---|
| 1112 | %Solve Sea-level equation on Earth only: {{{
|
|---|
| 1113 | md=sl.earth; %we don't do computations on ice sheets or land.
|
|---|
| 1114 |
|
|---|
| 1115 | %elastic loading from love numbers:
|
|---|
| 1116 | nlov=1001;
|
|---|
| 1117 | md.slr.love_h = love_numbers('h'); md.slr.love_h(nlov+1:end)=[];
|
|---|
| 1118 | md.slr.love_k = love_numbers('k'); md.slr.love_k(nlov+1:end)=[];
|
|---|
| 1119 | md.slr.love_l = love_numbers('l'); md.slr.love_l(nlov+1:end)=[];
|
|---|
| 1120 | md.slr.ocean_area_scaling = 0;
|
|---|
| 1121 |
|
|---|
| 1122 | %Miscellaneous
|
|---|
| 1123 | md.miscellaneous.name='test2004';
|
|---|
| 1124 |
|
|---|
| 1125 | %New stuff
|
|---|
| 1126 | md.slr.spcthickness = NaN(md.mesh.numberofvertices,1);
|
|---|
| 1127 | md.slr.Ngia = zeros(md.mesh.numberofvertices,1);
|
|---|
| 1128 | md.slr.Ugia = zeros(md.mesh.numberofvertices,1);
|
|---|
| 1129 |
|
|---|
| 1130 | %Solution parameters
|
|---|
| 1131 | md.slr.reltol=NaN;
|
|---|
| 1132 | md.slr.abstol=1e-3;
|
|---|
| 1133 | md.slr.geodetic=1;
|
|---|
| 1134 |
|
|---|
| 1135 | %eustatic + rigid + elastic run:
|
|---|
| 1136 | md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=0;
|
|---|
| 1137 | md.cluster=generic('name',oshostname(),'np',3);
|
|---|
| 1138 | %md.verbose=verbose('111111111');
|
|---|
| 1139 | md=solve(md,'Sealevelrise');
|
|---|
| 1140 | SnoRotation=md.results.SealevelriseSolution.Sealevel;
|
|---|
| 1141 |
|
|---|
| 1142 | %eustatic + rigid + elastic + rotation run:
|
|---|
| 1143 | md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=1;
|
|---|
| 1144 | md.cluster=generic('name',oshostname(),'np',3);
|
|---|
| 1145 | %md.verbose=verbose('111111111');
|
|---|
| 1146 | md=solve(md,'Sealevelrise');
|
|---|
| 1147 | SRotation=md.results.SealevelriseSolution.Sealevel;
|
|---|
| 1148 |
|
|---|
| 1149 | %Fields and tolerances to track changes
|
|---|
| 1150 | field_names ={'noRotation','Rotation'};
|
|---|
| 1151 | field_tolerances={1e-13,1e-13};
|
|---|
| 1152 | field_values={SnoRotation,SRotation};
|
|---|
| 1153 |
|
|---|
| 1154 | %}}}
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
