![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 24763
- Timestamp:
- 04/30/20 09:27:29 (5 years ago)
- File:
-
- 1 edited
-
issm/trunk-jpl/test/SandBox/test2004.m (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/test/SandBox/test2004.m
r24722 r24763 3 3 %Data paths: {{{ 4 4 modeldatapath='/Users/larour/ModelData'; 5 shppath='/Users/larour/issm-jpl/proj-group/qgis/ Slr';5 shppath='/Users/larour/issm-jpl/proj-group/qgis/NightlyRun'; 6 6 gshhsshapefile=[modeldatapath '/Gshhg/Shp/GSHHS_shp/c/GSHHS_c_L1.shp']; 7 7 %}}} 8 8 9 skipmesh=1; 10 9 11 %create sealevel model to hold our information: 10 sl=sealevelmodel(); 11 12 if ~skipmesh, 13 sl=sealevelmodel(); 12 14 %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)})); 15 %HemisphereWest: {{{ 16 sl.addbasin(basin('continent','hemispherewest','name','hemispherewest','epsg',3587,'boundaries',{... %Peru projection 3587 17 boundary('shppath',shppath,'shpfilename','HemisphereSplit','epsg',4326,'orientation','reverse'),... 18 boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031),... 19 boundary('shppath',shppath,'shpfilename','RonneBrunt','epsg',3031,'orientation','reverse'),... 20 boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031),... 21 boundary('shppath',shppath,'shpfilename','RonneFront','epsg',3031,'orientation','reverse'),... 22 boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031),... 23 boundary('shppath',shppath,'shpfilename','WestAntarctica2','epsg',3031,'orientation','reverse'),... 24 boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031)... 25 })); 26 %}}} 27 %HemisphereEast: {{{ 28 sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','epsg',4462,'boundaries',{... %Australian projection lat,long 29 boundary('shppath',shppath,'shpfilename','HemisphereSplit','epsg',4326),... 30 boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031),... 31 boundary('shppath',shppath,'shpfilename','EastAntarctica2','epsg',3031,'orientation','reverse'),... 32 boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031)... 33 })); 34 %}}} 35 %Antarctica excluding Ronne: {{{ 36 sl.addbasin(basin('continent','antarctica','name','antarctica-grounded','epsg',3031,'boundaries',{... 37 boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031),... 38 boundary('shppath',shppath,'shpfilename','EastAntarctica2','epsg',3031),... 39 boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031),... 40 boundary('shppath',shppath,'shpfilename','WestAntarctica2','epsg',3031)... 41 boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031)... 42 boundary('shppath',shppath,'shpfilename','RonneIceShelf','epsg',3031)... 43 boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031)... 44 boundary('shppath',shppath,'shpfilename','RonneBrunt','epsg',3031)... 45 })); 46 %}}} 47 %Ronne: {{{ 48 sl.addbasin(basin('continent','antarctica','name','ronne','epsg',3031,'boundaries',{... 49 boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031),... 50 boundary('shppath',shppath,'shpfilename','RonneIceShelf','epsg',3031),... 51 boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031),... 52 boundary('shppath',shppath,'shpfilename','RonneFront','epsg',3031,'orientation','reverse')... 53 })); 222 54 %}}} 223 55 %}}} … … 244 76 md=bamg(model,'domain',domain,'subdomains',coastline,'hmin',hmin,'hmax',hmax,defaultoptions{:}); 245 77 plotmodel(md,'data','mesh');pause(1); 246 error;247 78 248 79 %miscellaneous: … … 259 90 end 260 91 %}}} 261 e rror;92 end 262 93 %Parameterize ice sheets : {{{ 263 94 264 for ind=sl.basinindx('continent',{' greenland','antarctica'}),95 for ind=sl.basinindx('continent',{'antarctica'}), 265 96 disp(sprintf('Parameterizing basin %s\n', sl.icecaps{ind}.miscellaneous.name)); 266 97 267 98 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 end278 %}}}279 99 %masks : %{{{ 280 100 %new type of mask: 281 101 md.mask=maskpsl; 282 102 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; 103 %land and ocean: 104 md.mask.land_levelset=ones(md.mesh.numberofvertices,1); 105 md.mask.ocean_levelset=zeros(md.mesh.numberofvertices,1); 288 106 289 107 %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; 108 md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1); 109 pos=find(md.mesh.vertexonboundary); md.mask.ice_levelset(pos)=0; 110 111 %ice front: 297 112 pos=find(md.mask.ice_levelset==0); md.mask.ocean_levelset(pos)=-1; md.mask.land_levelset(pos)=1; 298 113 299 114 %now, glaciers: 300 mesh_glacier_levelset=md.private.bamg.landmask; %see meshing for why301 115 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 end310 116 %}}} 311 117 %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; 118 if bas.isnameany('antarctica-grounded'), 119 120 md.mask.groundedice_levelset=ones(md.mesh.numberofvertices,1); 121 end 122 if bas.isnameany('ronne'), 123 124 md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1); 330 125 331 126 %correction to land and ocean levelset: ice shelf is not on land! … … 333 128 md.mask.ocean_levelset(pos)=1; 334 129 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 130 end 347 131 %}}} … … 349 133 [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326'); 350 134 %}}} 351 %flowequation: % {{{352 md=setflowequation(md,'SSA','all');353 %}}}354 135 %geometry: {{{ 355 if bas.iscontinentany('antarctica'), % {{{136 if bas.iscontinentany('antarctica'), 356 137 di=md.materials.rho_ice/md.materials.rho_water; 357 138 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 139 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; 140 md.geometry.bed=interpBedmap2(md.mesh.x,md.mesh.y,'bed'); 141 end % }}} 142 %Slr: {{{ 143 if bas.iscontinentany('antarctica'), 144 145 delH=textread([modeldatapath '/AdhikariSciAdvTrends/AIS_delH_trend.txt']); 146 longAIS=delH(:,1); 147 latAIS=delH(:,2); 148 delHAIS=delH(:,3); 149 index=delaunay(latAIS,longAIS); 150 151 lat=md.mesh.lat; 152 long=md.mesh.long+360; 153 pos=find(long>360);long(pos)=long(pos)-360; 154 155 delHAIS=InterpFromMesh2d(index,longAIS,latAIS,delHAIS,long,lat); 156 northpole=find_point(md.mesh.long,md.mesh.lat,0,90); delHAIS(northpole)=0; 157 158 md.slr.deltathickness=delHAIS(md.mesh.elements)*[1;1;1]/3/100; 159 160 md.slr.sealevel=zeros(md.mesh.numberofvertices,1); 161 md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1); 162 md.slr.Ngia=zeros(md.mesh.numberofvertices,1); 163 md.slr.Ugia=zeros(md.mesh.numberofvertices,1); 164 165 md.dsl.global_average_thermosteric_sea_level_change=[0;0]; 166 md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); 167 md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); 168 md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1); 423 169 424 170 end %}}} 171 % material properties: {{{ 172 md.materials=materials('hydro'); 173 %}}} 174 %diverse: {{{ 175 md.miscellaneous.name=bas.name; 425 176 % }}} 426 %velocities: %{{{427 if bas.iscontinentany('antarctica'),% {{{428 429 velnc =[jplsvn '/database/antarctica/velocity_ref_v28Apr2011bamber_900m_median.nc'];430 431 %read and process netcdf432 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.drag490 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 -> 2010501 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 end517 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)<=1978563 ncdata = [MARname '/MARv3.5-15km-monthly-ERA-40-' num2str(Years(ii)) '.nc'];564 elseif Years(ii)>1978565 ncdata = [MARname '/MARv3.5-15km-monthly-ERA-Interim-' num2str(Years(ii)) '.nc'];566 end567 568 st1 = ncread(ncdata,'STcorr'); % Surface temperature569 st1(st1<-0.9000e30)=nan;570 st1 = double(mean(st1,3)); % Annual mean surface temp571 572 st1s(:,ii)=st1(:);573 574 smb1 = ncread(ncdata,'SMBcorr'); % in mmWe/month575 smb1(smb1<-0.900e30)=nan;576 smb1 = double(sum(smb1,3)/1000); %mWe/year577 578 smb1s(:,ii)=smb1(:);579 580 end581 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 end648 649 %prognostic650 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 else768 %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 end781 md.thermal.spctemperature=[md.initialization.temperature;1]; %impose observed temperature on surface782 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/yr792 %Rignot and Steffen, 2008 - Petermann high rates up to 20 km out793 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,35798 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]; %5807 meltratesgl=[22,4]; %26808 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, 2003818 else819 iceshelves={};820 end821 %}}}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; end842 end843 844 md.basalforcings.floatingice_melting_rate(find(in))=meltrate;845 md.basalforcings.floatingice_melting_rate(find(in&glpt))=meltrategl;846 847 end848 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 end857 end858 end859 %}}}860 %Thermal spcs and heatflux: %{{{861 md.thermal.spctemperature=md.initialization.temperature; %impose observed temperature on surface862 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 else887 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 end901 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 177 918 178 sl.icecaps{ind}=md; … … 921 181 % ParameterizeContinents {{{ 922 182 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'}), 183 for ind=sl.basinindx('continent',{'hemisphereeast','hemispherewest'}), 930 184 disp(sprintf('Masks for basin %s\n', sl.icecaps{ind}.miscellaneous.name)); 931 185 md=sl.icecaps{ind}; bas=sl.basins{ind}; … … 934 188 [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326'); 935 189 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 190 %mask: %{{{ 943 191 %new type of mask: … … 990 238 md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); 991 239 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 240 %grounded ice: 1005 241 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; 242 1012 243 % }}} 1013 244 %}}} 1014 245 %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: 246 1020 247 md.slr.sealevel=zeros(md.mesh.numberofvertices,1); 248 md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1); 249 md.slr.Ngia=zeros(md.mesh.numberofvertices,1); 250 md.slr.Ugia=zeros(md.mesh.numberofvertices,1); 251 252 md.dsl.global_average_thermosteric_sea_level_change=[0;0]; 253 %md.slr.steric_rate=(1.1+.38)*ones(md.mesh.numberofvertices,1); %steric + water storage. 254 md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); 255 md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); 256 md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1); 1021 257 1022 258 %love numbers: 1023 md.slr.steric_rate=(1.1+.38)*ones(md.mesh.numberofvertices,1); %steric + water storage.1024 259 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 260 %}}} 1028 261 %geometry: {{{ 1029 262 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; 263 md.geometry.bed=-ones(md.mesh.numberofvertices,1); 1045 264 % }}} 1046 265 %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; 266 md.materials=materials('hydro'); 1055 267 % }}} 1056 268 sl.icecaps{ind}=md; 1057 269 end 1058 270 % }}} 271 1059 272 %Assemble Earth in 3D {{{ 1060 273 … … 1066 279 %create earth model by concatenating all the icecaps in 3d: 1067 280 sl.caticecaps('tolerance',tolerance,'loneedgesdetect',loneedgesdetect); 281 282 error; 1068 283 1069 284 %figure out how each icecap's mesh connects to the larger earth mesh: … … 1110 325 1111 326 % }}} 327 error; 1112 328 %Solve Sea-level equation on Earth only: {{{ 1113 329 md=sl.earth; %we don't do computations on ice sheets or land.
Note:
See TracChangeset
for help on using the changeset viewer.
