Repository: MITgcm/xmitgcm

Issue:
How to find the lat/ lon values for LLC
Maybe this is still in the works, but is there some way to figure out the lat-lon values for the data? 

When I read in the data following the binder example (https://binder.pangeo.io/v2/gh/pangeo-data/pangeo-ecco-llc/master?filepath=visualize_SST.ipynb) there is no lat/lon information, just indices i,j. 

Maybe an example on how to do this would be useful. 

#158 #154 might be related.

Discussion:
rabernat: Tim thanks SO much for doing this! I plan to review this PR in detail ASAP. 

Sent from my iPhone

> On Apr 4, 2020, at 3:44 PM, Timothy Smith <notifications@github.com> wrote:
> 
> ﻿
> This PR is designed to give llcreader access to grid information on the ECCO portal (or elsewhere), to hopefully to close #161, close #166, and close #158 (possibly more).
> 
> Right now it gets the full grid for LLC 2160 and gets the available grid info on the ECCO portal for LLC 4320, which is some but not all. I was not able access the grid info mentioned by @rabernat in #161, which was here
> 
> vertical grid variables and coordinates
> this got a little tricky since everything is designed to deal with faceted data. I wrote some functions to read these in differently. I adapted the same infrastructure to the 1D variables, even though they could probably be dealt with eagerly (see _get_1d_chunk and _dask_array_vgrid in llcmodel.py).
> 
> k_p1 index
> I changed this implementation as follows. Currently, the code grabs k_levels for k_p1 where k_levels is the cell centered index. Even if the user does not select levels, this still results in len(k)==len(k_p1), which is technically not right, unless I'm missing something... It gets tricky though when the user selects individual, nonconsecutive levels. I wrote this PR so that, if the user selects cell centered k_levels=[0,10,11] then the output dataset has k_p1=[0,1,10,11,12] - so that k_p1 has the upper and lower indices for each vertical grid cell. I found this easier to do in order to correctly assign Zu,Zl, and Zp1. This implementation "seems right", but I'm interested if anyone has an idea for how to handle it more elegantly or if it's unnecessary/confusing. It also felt like a somewhat hacky solution, so again, feedback welcome. (see _get_kp1_levels in llcmodel.py). With this new implementation I had to change the expected size of k_p1 to be 51 and 91 for llc90 and llc2160/llc4320.
> 
> tests
> I just appended checks for this stuff in current tests rather than adding new ones. Tests for the k_p1 levels is in test_llc90_local_faces_load and 2160/4320 grids checked for in test_ecco_portal_faces/latlon. Perhaps could add more, testing the added read_grid and grid_vars_to_coords options added to get_dataset.
> 
> Some of the implementation started to feel a little hacky, and I'm open to anyone's suggestion on making it cleaner. Let me know what you think!
> 
> You can view, comment on, or merge this pull request online at:
> 
>   https://github.com/MITgcm/xmitgcm/pull/197
> 
> Commit Summary
> 
> access to llc2160 grid is a go, with a couple caveats for now
> use grid var names not filenames, get vertical grid vars working
> lil typo
> fixes loading grid dir issues
> it works for llc2160
> add directory to some LLC4320 grid vars
> separate grid metadata function, add shrunk_grid flag
> add grid_varnames to paramter list
> test that grid variables show up in portal datasets
> prefixes->variables
> try Zl...
> extra logic to create Zl, Zu, and Zp1
> make reading grid optional
> wow, not copying dicts really messes things up
> only build Zu, Zl if RF is there
> shouldn't kp1 have nz+1 levels
> only read grid if grid_varnames has entries
> also test selected k_p1 levels
> typo...
> ... more typos...
> rhoref, extra dx/dy metrics
> add new metrics to test
> typos, and add rhoRef to vgrid prefixes
> File Changes
> 
> M xmitgcm/llcreader/known_models.py (17)
> M xmitgcm/llcreader/llcmodel.py (201)
> M xmitgcm/llcreader/stores.py (41)
> M xmitgcm/test/test_llcreader.py (28)
> M xmitgcm/variables.py (28)
> Patch Links:
> 
> https://github.com/MITgcm/xmitgcm/pull/197.patch
> https://github.com/MITgcm/xmitgcm/pull/197.diff
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub, or unsubscribe.

rabernat: > I was not able access the grid info mentioned by @rabernat in #161, which was [here](https://pangeo-data.github.io/pangeo-datastore/master/ocean/llc4320.html#llc4320-grid)

FYI this has moved here: https://catalog.pangeo.io/browse/master/ocean/LLC4320/LLC4320_grid/
timothyas: Thanks for the review @rabernat! I will get to this in the next couple of nights. I agree, writing out chunking 1D variables felt ... awkward ... but it was a start. I'll change that around and address the comments as well.

Should I point the LLC4320 grid to the link you mentioned? https://catalog.pangeo.io/browse/master/ocean/LLC4320/LLC4320_grid/

I hope you're doing OK - staying safe, sane, & healthy!
rabernat: Thanks for the update Tim!

It looks like we are getting some test failures on python 3.6 only. I'm not sure what's behind this. They appear unrelated to your PR, but it would be good to try to resolve before merging.
timothyas: For sure! I'll start cracking into it tonight. I have been lazily ignoring it, but I agree it would be good to get to the bottom of before getting a new version on pip.
rabernat: Now that I've merged #200, could you rebase this? Hopefully the tests will all go green.
timothyas: Looks as good as green to me :) 

One last thing: should we keep the LLC4320 grid link to the ECCO portal, which has oddly an incomplete set of the coordinate files or switch to the Pangeo link?  The ECCO portal is currently missing e.g. Depth, but containing the seemingly necessary ones like dx/dy, xc/yc. Or maybe the ECCO portal could get updated ... 
rabernat: > should we keep the LLC4320 grid link to the ECCO portal, which has oddly an incomplete set of the coordinate files or switch to the Pangeo link?

Are you referring to this directory: https://data.nas.nasa.gov/ecco/data.php?dir=/eccodata/llc_4320/grid

Let's ping @menemenlis to see if we can get them to update the grid directory to include the complete set of files.
timothyas: Yep that's what I'm referring to and currently referencing in `known_models`. Compared to the 2160, that website is missing 

`AngleCS`, `AngleSN`, `Depth`, `rLowC/W/S`, `rSurfC/W/S` 

I think the `AngleCS` and `AngleSN` are probably the most critical out of that list. 
rabernat: I've followed up with Dimitris over email to try to get his input on this.
hongandyan: Are you looking for missing grid at ECCO NAS data portal?
FYI, ECCO Drive has complete grid
https://ecco.jpl.nasa.gov/drive/files/ECCO2/LLC4320/grid


rabernat: When we developed this tool, only ECCO Data Portal existed.

@hongandyan -- could you and / or @menemenlis explain the difference between "ECCO Drive" and "ECCO Data Portal"? Why would we use one vs. the other?

hongandyan: > When we developed this tool, only ECCO Data Portal existed.
> 
> @hongandyan -- could you and / or @menemenlis explain the difference between "ECCO Drive" and "ECCO Data Portal"? Why would we use one vs. the other?

sorry for the confusion.
To my understanding, ECCO Drive holds ECCO product, including LLC-xxx
while ECCO Data Portal hosts LLC-xxxx model output .
Your tools are great, forget about ECCO Drive.
Here I just want to provide some grid info as LLC4320-grid looks messed on ECCO Data Portal

rabernat: I just got the following info from @menemenlis:

> I have added all the grid files that I could locate to:
/nobackupp2/dmenemen/llc_2160/grid
and
/nobackupp2/dmenemen/llc_4320/grid
> Some from Hong’s ECCO drive and some that I recreated manually from the tile* files.
>
> I attach below the list of available grid files.  Are there additional grid files needed?
> In particular, Ryan, did you want me to put the compressed binary mask files there also?
> If there are more grid files needed and you have them available somewhere, please let me know.

<details>

```
pfe25 ~.$ cd /nobackupp2/dmenemen/llc_2160/grid
pfe25 grid.$ ls -l
total 77953672
-rw-r--r-- 1 dmenemen g26209   242611200 May 22  2014 AngleCS.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 AngleCS.meta
-rw-r--r-- 1 dmenemen g26209   242611200 May 22  2014 AngleSN.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 AngleSN.meta
-rw-r--r-- 1 dmenemen g26209         364 May 22  2014 DRC.data
-rw-r--r-- 1 dmenemen g26209         144 May 22  2014 DRC.meta
-rw-r--r-- 1 dmenemen g26209         360 May 22  2014 DRF.data
-rw-r--r-- 1 dmenemen g26209         144 May 22  2014 DRF.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 DXC.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 DXC.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 DXF.data
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 DXG.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 DXG.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 DXV.data
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 DYC.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 DYC.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 DYF.data
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 DYG.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 DYG.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 DYU.data
-rw-r--r-- 1 dmenemen g26209   242611200 May 22  2014 Depth.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 Depth.meta
-rw-r--r-- 1 dataset  n1856           15 May 21  2019 FileSystemUp.txt
-rw-r--r-- 1 dmenemen g26209         360 May 22  2014 PHrefC.data
-rw-r--r-- 1 dmenemen g26209         144 May 22  2014 PHrefC.meta
-rw-r--r-- 1 dmenemen g26209         364 May 22  2014 PHrefF.data
-rw-r--r-- 1 dmenemen g26209         144 May 22  2014 PHrefF.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 RAC.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 RAC.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 RAS.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 RAS.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 RAW.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 RAW.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 RAZ.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 RAZ.meta
-rw-r--r-- 1 dmenemen g26209         360 May 22  2014 RC.data
-rw-r--r-- 1 dmenemen g26209         144 May 22  2014 RC.meta
-rw-r--r-- 1 dmenemen g26209         364 May 22  2014 RF.data
-rw-r--r-- 1 dmenemen g26209         144 May 22  2014 RF.meta
-rw-r--r-- 1 dmenemen g26209         360 May 22  2014 RhoRef.data
-rw-r--r-- 1 dmenemen g26209         144 May 22  2014 RhoRef.meta
-rw-r--r-- 1 dmenemen g26209   242611200 May 22  2014 U2zonDir.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 U2zonDir.meta
-rw-r--r-- 1 dmenemen g26209   242611200 May 22  2014 V2zonDir.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 V2zonDir.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 XC.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 14:05 XC.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 XG.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 14:05 XG.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 YC.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 14:56 YC.meta
-rw-r--r-- 1 dmenemen g26209   242611200 Nov  6  2019 YG.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 14:07 YG.meta
-rw-r--r-- 1 dmenemen g26209 21835008000 May 22  2014 hFacC.data
-rw-r--r-- 1 dmenemen g26209         223 May 22  2014 hFacC.meta
-rw-r--r-- 1 dmenemen g26209 21835008000 May 22  2014 hFacS.data
-rw-r--r-- 1 dmenemen g26209         223 May 22  2014 hFacS.meta
-rw-r--r-- 1 dmenemen g26209 21835008000 May 22  2014 hFacW.data
-rw-r--r-- 1 dmenemen g26209         223 May 22  2014 hFacW.meta
-rw-r--r-- 1 dmenemen g26209   242611200 May 22  2014 rLowC.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 rLowC.meta
-rw-r--r-- 1 dmenemen g26209   242611200 May 22  2014 rLowS.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 rLowS.meta
-rw-r--r-- 1 dmenemen g26209   242611200 May 22  2014 rLowW.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 rLowW.meta
-rw-r--r-- 1 dmenemen g26209   242611200 May 22  2014 rSurfC.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 rSurfC.meta
-rw-r--r-- 1 dmenemen g26209   242611200 May 22  2014 rSurfS.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 rSurfS.meta
-rw-r--r-- 1 dmenemen g26209   242611200 May 22  2014 rSurfW.data
-rw-r--r-- 1 dmenemen g26209         154 May 22  2014 rSurfW.meta
-rw-r--r-- 1 dmenemen g26209        1167 Nov  6  2019 read_tile_mitgrid.m
-rw-r--r-- 1 dmenemen g26209          43 Sep  9  2019 readme.txt
-rw-r--r-- 1 dmenemen g26209        1143 Jan  6 12:37 thk90.mat
-rw-r--r-- 1 dmenemen g26209  1792696448 Jun  1 14:57 tile001.mitgrid
-rw-r--r-- 1 dmenemen g26209  1792696448 Jun  1 15:16 tile002.mitgrid
-rw-r--r-- 1 dmenemen g26209   597749888 Jun  1 15:13 tile003.mitgrid
-rw-r--r-- 1 dmenemen g26209  1792696448 Jun  1 14:26 tile004.mitgrid
-rw-r--r-- 1 dmenemen g26209  1792696448 Jun  1 14:58 tile005.mitgrid
pfe25 grid.$ cd /nobackupp2/dmenemen/llc_4320/grid
pfe25 grid.$ ls -l
total 275781536
-rw-r--r-- 1 dmenemen g26209   970444800 Jun  2 02:32 AngleCS.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  2 02:31 AngleCS.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Jun  1 21:49 AngleSN.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 21:35 AngleSN.meta
-rw-r--r-- 1 dmenemen g26209         364 Jan  6 11:39 DRC.data
-rw-r--r-- 1 dmenemen g26209         144 Jan  6 11:39 DRC.meta
-rw-r--r-- 1 dmenemen g26209         360 Jan  6 11:39 DRF.data
-rw-r--r-- 1 dmenemen g26209         144 Jan  6 11:39 DRF.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 DXC.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 21:25 DXC.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 DXF.data
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 DXG.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 21:39 DXG.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 DXV.data
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 DYC.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 14:27 DYC.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 DYF.data
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 DYG.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 17:48 DYG.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 DYU.data
-rw-r--r-- 1 dmenemen g26209   970444800 Mar  9 16:28 Depth.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  2 02:26 Depth.meta
-rw-r--r-- 1 dataset  n1856           15 May 21  2019 FileSystemUp.txt
-rw-r--r-- 1 dmenemen g26209         360 Jan  6 11:40 PHrefC.data
-rw-r--r-- 1 dmenemen g26209         144 Jan  6 11:40 PHrefC.meta
-rw-r--r-- 1 dmenemen g26209         364 Jan  6 11:40 PHrefF.data
-rw-r--r-- 1 dmenemen g26209         144 Jan  6 11:40 PHrefF.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 RAC.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  2 02:34 RAC.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 RAS.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  2 02:26 RAS.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 RAW.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  2 02:31 RAW.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 RAZ.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  2 02:32 RAZ.meta
-rw-r--r-- 1 dmenemen g26209         360 Jan  6 11:40 RC.data
-rw-r--r-- 1 dmenemen g26209         144 Jan  6 11:40 RC.meta
-rw-r--r-- 1 dmenemen g26209         364 Jan  6 11:40 RF.data
-rw-r--r-- 1 dmenemen g26209         144 Jan  6 11:40 RF.meta
-rw-r--r-- 1 dmenemen g26209         360 Jan  6 11:40 RhoRef.data
-rw-r--r-- 1 dmenemen g26209         144 Jan  6 11:40 RhoRef.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Jun  2 02:31 U2zonDir.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  2 02:31 U2zonDir.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Jun  2 02:12 V2zonDir.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  2 02:17 V2zonDir.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 XC.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 21:25 XC.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 XG.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 17:43 XG.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 YC.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 21:44 YC.meta
-rw-r--r-- 1 dmenemen g26209   970444800 Sep 22  2019 YG.data
-rw-r--r-- 1 dmenemen g26209         154 Jun  1 21:30 YG.meta
-r--r--r-- 1 dmenemen g26209 87340032000 Feb 25 16:22 hFacC.data
-rw-r--r-- 1 dmenemen g26209         223 Jun  2 02:03 hFacC.meta
-r--r--r-- 1 dmenemen g26209 87340032000 Feb 25 17:36 hFacS.data
-rw-r--r-- 1 dmenemen g26209         223 Jun  1 21:25 hFacS.meta
-r--r--r-- 1 dmenemen g26209 87340032000 Feb 25 17:45 hFacW.data
-rw-r--r-- 1 dmenemen g26209         223 Jun  1 17:43 hFacW.meta
-rw-r--r-- 1 dmenemen g26209         307 Mar  9 16:16 mk_Depth.data.m
-rw-r--r-- 1 dmenemen g26209        1138 Sep 22  2019 read_tile_mitgrid.m
-rw-r--r-- 1 dmenemen g26209        4254 Jun  1 21:30 readme.txt
-rw-r--r-- 1 dmenemen g26209        1143 Jan  6 11:43 thk90.mat
```

</details>

@timothyas - can you check that things are working on your branch with these new additions?

Patch:
diff --git a/doc/llcreader.rst b/doc/llcreader.rst
index b7e02d0..372a8f0 100644
--- a/doc/llcreader.rst
+++ b/doc/llcreader.rst
@@ -98,7 +98,7 @@ By default, all variables and all timesteps are loaded::
   >>> ds = model.get_dataset(k_chunksize=90)
   >>> print(ds)
   <xarray.Dataset>
-  Dimensions:   (face: 13, i: 4320, i_g: 4320, j: 4320, j_g: 4320, k: 90, k_l: 90, k_p1: 90, k_u: 90, time: 9030)
+  Dimensions:   (face: 13, i: 4320, i_g: 4320, j: 4320, j_g: 4320, k: 90, k_l: 90, k_p1: 91, k_u: 90, time: 9030)
   Coordinates:
     * face      (face) int64 0 1 2 3 4 5 6 7 8 9 10 11 12
     * i         (i) int64 0 1 2 3 4 5 6 7 ... 4313 4314 4315 4316 4317 4318 4319
@@ -108,30 +108,54 @@ By default, all variables and all timesteps are loaded::
     * k         (k) int64 0 1 2 3 4 5 6 7 8 9 10 ... 80 81 82 83 84 85 86 87 88 89
     * k_u       (k_u) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
     * k_l       (k_l) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
-    * k_p1      (k_p1) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
+    * k_p1      (k_p1) int64 0 1 2 3 4 5 6 7 8 9 ... 81 82 83 84 85 86 87 88 89 90
       niter     (time) int64 ...
     * time      (time) datetime64[ns] 2011-09-13 ... 2012-09-23T05:00:00
+      drC       (k_p1) >f4 dask.array<chunksize=(91,), meta=np.ndarray>
+      drF       (k) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
+      dxC       (face, j, i_g) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dxF       (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dxG       (face, j_g, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dyC       (face, j_g, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dyF       (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dyG       (face, j, i_g) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      hFacC     (k, face, j, i) float32 dask.array<chunksize=(90, 3, 4320, 4320), meta=np.ndarray>
+      hFacS     (k, face, j_g, i) float32 dask.array<chunksize=(90, 3, 4320, 4320), meta=np.ndarray>
+      hFacW     (k, face, j, i_g) float32 dask.array<chunksize=(90, 3, 4320, 4320), meta=np.ndarray>
+      PHrefC    (k) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
+      PHrefF    (k_p1) >f4 dask.array<chunksize=(91,), meta=np.ndarray>
+      rA        (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      rAs       (face, j_g, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      rAw       (face, j, i_g) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      Z         (k) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
+      Zp1       (k_p1) >f4 dask.array<chunksize=(91,), meta=np.ndarray>
+      rhoRef    (k) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
+      XC        (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      YC        (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      Zl        (k_l) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
+      Zu        (k_u) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
   Data variables:
-      Eta       (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      KPPhbl    (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      oceFWflx  (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      oceQnet   (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      oceQsw    (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      oceSflux  (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      oceTAUX   (time, face, j, i_g) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      oceTAUY   (time, face, j_g, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      PhiBot    (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      Salt      (time, k, face, j, i) >f4 dask.array<shape=(9030, 90, 13, 4320, 4320), chunksize=(1, 90, 3, 4320, 4320)>
-      SIarea    (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      SIheff    (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      SIhsalt   (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      SIhsnow   (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      SIuice    (time, face, j, i_g) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      SIvice    (time, face, j_g, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
-      Theta     (time, k, face, j, i) >f4 dask.array<shape=(9030, 90, 13, 4320, 4320), chunksize=(1, 90, 3, 4320, 4320)>
-      U         (time, k, face, j, i_g) >f4 dask.array<shape=(9030, 90, 13, 4320, 4320), chunksize=(1, 90, 3, 4320, 4320)>
-      V         (time, k, face, j_g, i) >f4 dask.array<shape=(9030, 90, 13, 4320, 4320), chunksize=(1, 90, 3, 4320, 4320)>
-      W         (time, k_l, face, j, i) >f4 dask.array<shape=(9030, 90, 13, 4320, 4320), chunksize=(1, 90, 3, 4320, 4320)>
+      Eta       (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      KPPhbl    (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      oceFWflx  (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      oceQnet   (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      oceQsw    (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      oceSflux  (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      oceTAUX   (time, face, j, i_g) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      oceTAUY   (time, face, j_g, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      PhiBot    (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      Salt      (time, k, face, j, i) float32 dask.array<chunksize=(1, 90, 3, 4320, 4320), meta=np.ndarray>
+      SIarea    (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      SIheff    (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      SIhsalt   (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      SIhsnow   (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      SIuice    (time, face, j, i_g) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      SIvice    (time, face, j_g, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      Theta     (time, k, face, j, i) float32 dask.array<chunksize=(1, 90, 3, 4320, 4320), meta=np.ndarray>
+      U         (time, k, face, j, i_g) float32 dask.array<chunksize=(1, 90, 3, 4320, 4320), meta=np.ndarray>
+      V         (time, k, face, j_g, i) float32 dask.array<chunksize=(1, 90, 3, 4320, 4320), meta=np.ndarray>
+      W         (time, k_l, face, j, i) float32 dask.array<chunksize=(1, 90, 3, 4320, 4320), meta=np.ndarray>
+
 
 This dataset is useless for computations on a laptop, because the individual
 chunks require nearly 20 GB of memory.
@@ -141,8 +165,9 @@ Get a single 2D variable::
 
   >>> ds = model.get_dataset(varnames=['Eta'])
   >>> print(ds)
+
   <xarray.Dataset>
-  Dimensions:  (face: 13, i: 4320, i_g: 4320, j: 4320, j_g: 4320, k: 90, k_l: 90, k_p1: 90, k_u: 90, time: 9030)
+  Dimensions:  (face: 13, i: 4320, i_g: 4320, j: 4320, j_g: 4320, k: 90, k_l: 90, k_p1: 91, k_u: 90, time: 9030)
   Coordinates:
     * face     (face) int64 0 1 2 3 4 5 6 7 8 9 10 11 12
     * i        (i) int64 0 1 2 3 4 5 6 7 ... 4313 4314 4315 4316 4317 4318 4319
@@ -152,18 +177,42 @@ Get a single 2D variable::
     * k        (k) int64 0 1 2 3 4 5 6 7 8 9 10 ... 80 81 82 83 84 85 86 87 88 89
     * k_u      (k_u) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
     * k_l      (k_l) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
-    * k_p1     (k_p1) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
+    * k_p1     (k_p1) int64 0 1 2 3 4 5 6 7 8 9 ... 81 82 83 84 85 86 87 88 89 90
       niter    (time) int64 ...
     * time     (time) datetime64[ns] 2011-09-13 ... 2012-09-23T05:00:00
+      drC      (k_p1) >f4 dask.array<chunksize=(91,), meta=np.ndarray>
+      drF      (k) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
+      dxC      (face, j, i_g) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dxF      (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dxG      (face, j_g, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dyC      (face, j_g, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dyF      (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dyG      (face, j, i_g) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      hFacC    (k, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      hFacS    (k, face, j_g, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      hFacW    (k, face, j, i_g) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      PHrefC   (k) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
+      PHrefF   (k_p1) >f4 dask.array<chunksize=(91,), meta=np.ndarray>
+      rA       (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      rAs      (face, j_g, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      rAw      (face, j, i_g) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      Z        (k) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
+      Zp1      (k_p1) >f4 dask.array<chunksize=(91,), meta=np.ndarray>
+      rhoRef   (k) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
+      XC       (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      YC       (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      Zl       (k_l) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
+      Zu       (k_u) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
   Data variables:
-      Eta      (time, face, j, i) >f4 dask.array<shape=(9030, 13, 4320, 4320), chunksize=(1, 3, 4320, 4320)>
+      Eta      (time, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+
 
 Get a few vertical levels from some 3D variables::
 
   >>> ds = model.get_dataset(varnames=['Salt', 'Theta'], k_levels=[1, 10, 40])
   >>> print(ds)
   <xarray.Dataset>
-  Dimensions:  (face: 13, i: 4320, i_g: 4320, j: 4320, j_g: 4320, k: 3, k_l: 3, k_p1: 3, k_u: 3, time: 9030)
+  Dimensions:  (face: 13, i: 4320, i_g: 4320, j: 4320, j_g: 4320, k: 3, k_l: 3, k_p1: 6, k_u: 3, time: 9030)
   Coordinates:
     * face     (face) int64 0 1 2 3 4 5 6 7 8 9 10 11 12
     * i        (i) int64 0 1 2 3 4 5 6 7 ... 4313 4314 4315 4316 4317 4318 4319
@@ -173,13 +222,39 @@ Get a few vertical levels from some 3D variables::
     * k        (k) int64 1 10 40
     * k_u      (k_u) int64 1 10 40
     * k_l      (k_l) int64 1 10 40
-    * k_p1     (k_p1) int64 1 10 40
+    * k_p1     (k_p1) int64 1 2 10 11 40 41
       niter    (time) int64 ...
     * time     (time) datetime64[ns] 2011-09-13 ... 2012-09-23T05:00:00
+      drC      (k_p1) >f4 dask.array<chunksize=(6,), meta=np.ndarray>
+      drF      (k) >f4 dask.array<chunksize=(3,), meta=np.ndarray>
+      dxC      (face, j, i_g) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dxF      (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dxG      (face, j_g, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dyC      (face, j_g, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dyF      (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      dyG      (face, j, i_g) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      hFacC    (k, face, j, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      hFacS    (k, face, j_g, i) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      hFacW    (k, face, j, i_g) float32 dask.array<chunksize=(1, 3, 4320, 4320), meta=np.ndarray>
+      PHrefC   (k) >f4 dask.array<chunksize=(3,), meta=np.ndarray>
+      PHrefF   (k_p1) >f4 dask.array<chunksize=(6,), meta=np.ndarray>
+      rA       (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      rAs      (face, j_g, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      rAw      (face, j, i_g) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      Z        (k) >f4 dask.array<chunksize=(3,), meta=np.ndarray>
+      Zp1      (k_p1) >f4 dask.array<chunksize=(6,), meta=np.ndarray>
+      rhoRef   (k) >f4 dask.array<chunksize=(3,), meta=np.ndarray>
+      XC       (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      YC       (face, j, i) float32 dask.array<chunksize=(3, 4320, 4320), meta=np.ndarray>
+      Zl       (k_l) >f4 dask.array<chunksize=(3,), meta=np.ndarray>
+      Zu       (k_u) >f4 dask.array<chunksize=(3,), meta=np.ndarray>
   Data variables:
-      Salt     (time, k, face, j, i) >f4 dask.array<shape=(9030, 3, 13, 4320, 4320), chunksize=(1, 1, 3, 4320, 4320)>
-      Theta    (time, k, face, j, i) >f4 dask.array<shape=(9030, 3, 13, 4320, 4320), chunksize=(1, 1, 3, 4320, 4320)>
+      Salt     (time, k, face, j, i) float32 dask.array<chunksize=(1, 1, 3, 4320, 4320), meta=np.ndarray>
+      Theta    (time, k, face, j, i) float32 dask.array<chunksize=(1, 1, 3, 4320, 4320), meta=np.ndarray>
 
+Note that when vertical levels are subset like this, any vertical coordinate
+associated with dimension `k_p1` will have levels above and below the selected
+`k_levels`, which are at cell center.
 
 A list of all available variables can be seen as follows::
 
diff --git a/xmitgcm/llcreader/known_models.py b/xmitgcm/llcreader/known_models.py
index 3b2ea59..3b6f4ac 100644
--- a/xmitgcm/llcreader/known_models.py
+++ b/xmitgcm/llcreader/known_models.py
@@ -38,6 +38,15 @@ class LLC2160Model(BaseLLCModel):
     varnames = ['Eta', 'KPPhbl', 'oceFWflx', 'oceQnet', 'oceQsw', 'oceSflux',
                 'oceTAUX', 'oceTAUY', 'PhiBot', 'Salt', 'SIarea', 'SIheff',
                 'SIhsalt', 'SIhsnow', 'SIuice', 'SIvice', 'Theta', 'U', 'V', 'W']
+    grid_varnames = ['AngleCS','AngleSN','Depth',
+                     'DRC','DRF',
+                     'DXC','DXF','DXG',
+                     'DYC','DYF','DYG',
+                     'hFacC','hFacS','hFacW','PHrefC','PHrefF',
+                     'RAC','RAS','RAW',
+                     'RC','RF','RhoRef','rLowC','rLowS','rLowW',
+                     'rSurfC','rSurfS','rSurfW','XC','YC']
+    # corner point problems: 'RAZ','XG','YG','DXV','DYU'
     mask_override = {'oceTAUX': 'c', 'oceTAUY': 'c'}
 
 
@@ -46,13 +55,21 @@ class LLC4320Model(BaseLLCModel):
     nz = 90
     delta_t = 25
     iter_start = 10368
-    iter_stop = 1310544 + 1
+    iter_stop = 1495152 + 1
     iter_step = 144
     time_units='seconds since 2011-09-10'
     calendar = 'gregorian'
     varnames = ['Eta', 'KPPhbl', 'oceFWflx', 'oceQnet', 'oceQsw', 'oceSflux',
                 'oceTAUX', 'oceTAUY', 'PhiBot', 'Salt', 'SIarea', 'SIheff',
                 'SIhsalt', 'SIhsnow', 'SIuice', 'SIvice', 'Theta', 'U', 'V', 'W']
+    grid_varnames = ['AngleCS','AngleSN','Depth',
+                     'DRC','DRF',
+                     'DXC','DXF','DXG',
+                     'DYC','DYF','DYG',
+                     'hFacC','hFacS','hFacW','PHrefC','PHrefF',
+                     'RAC','RAS','RAW','RC','RF',
+                     'RhoRef','XC','YC']
+    # corner point problems: 'RAZ','XG','YG','DXV','DYU'
     mask_override = {'oceTAUX': 'c', 'oceTAUY': 'c'}
 
 
@@ -61,9 +78,10 @@ class ECCOPortalLLC2160Model(LLC2160Model):
     def __init__(self):
         fs = _make_http_filesystem()
         base_path = 'https://data.nas.nasa.gov/ecco/download_data.php?file=/eccodata/llc_2160/compressed'
+        grid_path = 'https://data.nas.nasa.gov/ecco/download_data.php?file=/eccodata/llc_2160/grid'
         mask_path = 'https://storage.googleapis.com/pangeo-ecco/llc/masks/llc_2160_masks.zarr/'
         store = stores.NestedStore(fs, base_path=base_path, mask_path=mask_path,
-                                   shrunk=True, join_char='/')
+                                   grid_path=grid_path, shrunk=True, join_char='/')
         super(ECCOPortalLLC2160Model, self).__init__(store)
 
 
@@ -72,9 +90,10 @@ class ECCOPortalLLC4320Model(LLC4320Model):
     def __init__(self):
         fs = _make_http_filesystem()
         base_path = 'https://data.nas.nasa.gov/ecco/download_data.php?file=/eccodata/llc_4320/compressed'
+        grid_path = 'https://data.nas.nasa.gov/ecco/download_data.php?file=/eccodata/llc_4320/grid'
         mask_path = 'https://storage.googleapis.com/pangeo-ecco/llc/masks/llc_4320_masks.zarr/'
         store = stores.NestedStore(fs, base_path=base_path, mask_path=mask_path,
-                                   shrunk=True,  join_char='/')
+                                   grid_path=grid_path, shrunk=True, join_char='/')
         super(ECCOPortalLLC4320Model, self).__init__(store)
 
 
diff --git a/xmitgcm/llcreader/llcmodel.py b/xmitgcm/llcreader/llcmodel.py
index 7b30274..0f4a093 100644
--- a/xmitgcm/llcreader/llcmodel.py
+++ b/xmitgcm/llcreader/llcmodel.py
@@ -8,6 +8,34 @@ import warnings
 from .duck_array_ops import concatenate
 from .shrunk_index import all_index_data
 
+def _get_grid_metadata():
+    # keep this separate from get_var_metadata
+    # because grid stuff is weird
+    from ..mds_store import _get_all_grid_variables
+    from ..variables import extra_grid_variables, vertical_coordinates
+
+    # get grid info
+    grid_vars = _get_all_grid_variables('llc')
+    grid_vars.update(extra_grid_variables)
+
+    # make dictionary with keys as filenames
+    grid_metadata = {}
+    for key,val in grid_vars.items():
+        # masks use hFac filename to be computed in mds_store
+        if 'filename' in val and key[:4]!='mask':
+            val.update({'real_name':key})
+            grid_metadata[val['filename']] = val
+        else:
+            grid_metadata[key] = val
+
+    # force RF to point to Zp1, deal with this manually..
+    grid_metadata['RF']=vertical_coordinates['Zp1']
+    grid_metadata['RF']['real_name'] = 'Zp1'
+    for zv in ['Zu','Zl']:
+        grid_metadata[zv] = vertical_coordinates[zv]
+
+    return grid_metadata
+
 def _get_var_metadata():
     # The LLC run data comes with zero metadata. So we import metadata from
     # the xmitgcm package.
@@ -18,7 +46,7 @@ def _get_var_metadata():
 
     diag_file = StringIO(diagnostics)
     available_diags = parse_available_diagnostics(diag_file)
-    var_metadata = state_variables
+    var_metadata = state_variables.copy()
     var_metadata.update(package_state_variables)
     var_metadata.update(available_diags)
 
@@ -28,10 +56,18 @@ def _get_var_metadata():
     for a, b in aliases.items():
         var_metadata[a] = var_metadata[b]
 
+    # add grid metadata
+    var_metadata.update(_get_grid_metadata())
+
     return var_metadata
 
 _VAR_METADATA = _get_var_metadata()
 
+def _is_vgrid(vname):
+    # check for 1d, vertical grid variables
+    dims = _VAR_METADATA[vname]['dims']
+    return len(dims)==1 and dims[0][0]=='k'
+
 def _get_variable_point(vname, mask_override):
     # fix for https://github.com/MITgcm/xmitgcm/issues/191
     if vname in mask_override:
@@ -350,14 +386,17 @@ def _get_facet_chunk(store, varname, iternum, nfacet, klevels, nx, nz, dtype,
 
     file = fs.open(path)
 
-    # insert singleton axis for time and k level
-    facet_shape = (1, 1,) + _facet_shape(nfacet, nx)
+    # insert singleton axis for time (if not grid var) and k level
+    facet_shape = (1,) + _facet_shape(nfacet, nx)
+    facet_shape = (1,) + facet_shape if iternum is not None else facet_shape
 
     level_data = []
 
     # the store tells us whether we need a mask or not
     point = _get_variable_point(varname, mask_override)
-    if store.shrunk:
+
+    if (store.shrunk and iternum is not None) or \
+       (store.shrunk_grid and iternum is None):
         index = all_index_data[nx][point]
         zgroup = store.open_mask_group()
         mask = zgroup['mask_' + point].astype('bool')
@@ -399,6 +438,20 @@ def _get_facet_chunk(store, varname, iternum, nfacet, klevels, nx, nz, dtype,
 
     return np.concatenate(level_data, axis=1)
 
+def _get_1d_chunk(store, varname, klevels, nz, dtype):
+    """for 1D vertical grid variables"""
+
+    fs, path = store.get_fs_and_full_path(varname, None)
+
+    file = fs.open(path)
+
+    # read all levels for 1D variables
+    read_length = nz*dtype.itemsize # all levels in bytes
+    buffer = file.read(read_length)
+    data = np.frombuffer(buffer,dtype=dtype)
+
+    # now subset
+    return data[klevels]
 
 class BaseLLCModel:
     """Class representing an LLC Model Dataset.
@@ -428,8 +481,8 @@ class BaseLLCModel:
         Final model iteration number (exclusive; follows python range conventions)
     iter_step : int
         Spacing between iterations
-    varnames : list
-        List of variable names contained in the dataset
+    varnames, grid_varnames : list
+        List of data variable and grid variable names contained in the dataset
     mask_override : dict
         Override inference of masking variable, e.g. ``{'oceTAUX': 'c'}``
     """
@@ -445,6 +498,7 @@ class BaseLLCModel:
     iter_stop = None
     iter_step = None
     varnames = []
+    grid_varnames = []
     mask_override = {}
 
     def __init__(self, store):
@@ -454,7 +508,7 @@ class BaseLLCModel:
         ----------
         store : llcreader.BaseStore
         mask_ds : zarr.Group
-            Must contain variables `mask_c`, `masc_w`, `mask_s`
+            Must contain variables `mask_c`, `mask_w`, `mask_s`
         """
         self.store = store
         self.shape = (self.nz, self.nface, self.nx, self.nx)
@@ -475,6 +529,21 @@ class BaseLLCModel:
             masks[point] = _faces_to_facets(mask_faces)
         return masks
 
+    def _get_kp1_levels(self,k_levels):
+        # determine kp1 levels
+        # get borders to all k (center) levels
+        # ki used to get Zu, Zl later
+        ku = np.concatenate([k_levels[1:],[k_levels[-1]+1]])
+        kp1 = []
+        ki=[]
+        for i,(x,y) in enumerate(zip(k_levels,ku)):
+            kp1+= [x] if x not in kp1 else []
+            kp1+= [y] if y-x==1 else [x+1]
+
+
+        kp1=np.array(kp1)
+
+        return kp1
 
     def _make_coords_faces(self, all_iters):
         time = self.delta_t * all_iters
@@ -525,7 +594,7 @@ class BaseLLCModel:
     def _dask_array(self, nfacet, varname, iters, klevels, k_chunksize):
         # return a dask array for a single facet
         facet_shape =  _facet_shape(nfacet, self.nx)
-        time_chunks = (len(iters) * (1,),)
+        time_chunks = (len(iters) * (1,),) if iters is not None else ()
         k_chunks = (tuple([len(c)
                           for c in _chunks(klevels, k_chunksize)]),)
         chunks = time_chunks + k_chunks + tuple([(s,) for s in facet_shape])
@@ -534,16 +603,45 @@ class BaseLLCModel:
         dsk = {}
         token = tokenize(varname, self.store, nfacet)
         name = '-'.join([varname, token])
-        for n_iter, iternum in enumerate(iters):
+
+        # iters == None for grid variables
+        if iters is not None:
+            for n_iter, iternum in enumerate(iters):
+                for n_k, these_klevels in enumerate(_chunks(klevels, k_chunksize)):
+                    key = name, n_iter, n_k, 0, 0, 0
+                    task = (_get_facet_chunk, self.store, varname, iternum,
+                             nfacet, these_klevels, self.nx, self.nz, self.dtype,
+                             self.mask_override)
+                    dsk[key] = task
+        else:
             for n_k, these_klevels in enumerate(_chunks(klevels, k_chunksize)):
-                key = name, n_iter, n_k, 0, 0, 0
-                task = (_get_facet_chunk, self.store, varname, iternum,
+                key = name, n_k, 0, 0, 0
+                task = (_get_facet_chunk, self.store, varname, None,
                          nfacet, these_klevels, self.nx, self.nz, self.dtype,
                          self.mask_override)
                 dsk[key] = task
 
         return dsa.Array(dsk, name, chunks, self.dtype)
 
+    def _dask_array_vgrid(self, varname, klevels, k_chunksize):
+        # return a dask array for a 1D vertical grid var
+
+        # single chunk for 1D variables
+        chunks = ((len(klevels),),)
+
+        # manually build dask graph
+        dsk = {}
+        token = tokenize(varname, self.store)
+        name = '-'.join([varname, token])
+
+        nz = self.nz if _VAR_METADATA[varname]['dims'] != ['k_p1'] else self.nz+1
+        task = (_get_1d_chunk, self.store, varname,
+                list(klevels), nz, self.dtype)
+
+        key = name, 0
+        dsk[key] = task
+
+        return dsa.Array(dsk, name, chunks, self.dtype)
 
     def _get_facet_data(self, varname, iters, klevels, k_chunksize):
         mask, index = self._get_mask_and_index_for_variable(varname)
@@ -553,7 +651,10 @@ class BaseLLCModel:
         if len(dims)==2:
             klevels = [0,]
 
-        data_facets = [self._dask_array(nfacet, varname, iters, klevels, k_chunksize)
+        if _is_vgrid(varname):
+            data_facets = self._dask_array_vgrid(varname,klevels,k_chunksize)
+        else:
+            data_facets = [self._dask_array(nfacet, varname, iters, klevels, k_chunksize)
                        for nfacet in range(5)]
 
         if len(dims)==2:
@@ -565,7 +666,7 @@ class BaseLLCModel:
 
     def get_dataset(self, varnames=None, iter_start=None, iter_stop=None,
                     iter_step=None, k_levels=None, k_chunksize=1,
-                    type='faces'):
+                    type='faces', read_grid=True, grid_vars_to_coords=True):
         """
         Create an xarray Dataset object for this model.
 
@@ -588,6 +689,10 @@ class BaseLLCModel:
             How many vertical levels per Dask chunk.
         type : {'faces', 'latlon'}, optional
             What type of dataset to create
+        read_grid : bool, optional
+            Whether to read the grid info
+        grid_vars_to_coords : bool, optional
+            Whether to promote grid variables to coordinate status
 
         Returns
         -------
@@ -613,18 +718,32 @@ class BaseLLCModel:
 
         varnames = varnames or self.varnames
 
+        # grid stuff
+        read_grid = read_grid and len(self.grid_varnames)!=0
+        grid_vars_to_coords = grid_vars_to_coords and read_grid
+        grid_varnames = self.grid_varnames if read_grid else []
+
         ds = self._make_coords_faces(iters)
         if type=='latlon':
             ds = _faces_coords_to_latlon(ds)
 
         k_levels = k_levels or np.arange(self.nz)
-        ds = ds.sel(k=k_levels, k_l=k_levels, k_u=k_levels, k_p1=k_levels)
+        kp1_levels = self._get_kp1_levels(k_levels)
+
+        ds = ds.sel(k=k_levels, k_l=k_levels, k_u=k_levels, k_p1=kp1_levels)
 
         # get the data in facet form
         data_facets = {vname:
                        self._get_facet_data(vname, iters, k_levels, k_chunksize)
                        for vname in varnames}
 
+        # get the grid in facet form
+        # do separately for vertical coords on kp1_levels
+        grid_facets = {}
+        for vname in grid_varnames:
+            my_k_levels = k_levels if _VAR_METADATA[vname]['dims'] !=['k_p1'] else kp1_levels
+            grid_facets[vname] = self._get_facet_data(vname, None, my_k_levels, k_chunksize)
+
         # transform it into faces or latlon
         data_transformers = {'faces': _all_facets_to_faces,
                              'latlon': _all_facets_to_latlon}
@@ -632,15 +751,45 @@ class BaseLLCModel:
         transformer = data_transformers[type]
         data = transformer(data_facets, _VAR_METADATA)
 
+        # separate horizontal and vertical grid variables
+        hgrid_facets = {key: grid_facets[key]
+                for key in grid_varnames if not _is_vgrid(key)}
+        vgrid_facets = {key: grid_facets[key]
+                for key in grid_varnames if _is_vgrid(key)}
+
+        # do not transform vertical grid variables
+        data.update(transformer(hgrid_facets, _VAR_METADATA))
+        data.update(vgrid_facets)
+
         variables = {}
-        for vname in varnames:
+        gridlist = ['Zl','Zu'] if read_grid else []
+        for vname in varnames+grid_varnames:
             meta = _VAR_METADATA[vname]
             dims = meta['dims']
             if type=='faces':
                 dims = _add_face_to_dims(dims)
-            dims = ['time',] + dims
+            dims = ['time',] + dims if vname not in grid_varnames else dims
             attrs = meta['attrs']
-            variables[vname] = xr.Variable(dims, data[vname], attrs)
+
+            # Handle grid names different from filenames
+            fname = vname
+            vname = meta['real_name'] if 'real_name' in meta else vname
+            if fname in grid_varnames:
+                gridlist.append(vname)
+
+            variables[vname] = xr.Variable(dims, data[fname], attrs)
+
+        # handle vertical coordinate after the fact
+        if read_grid and 'RF' in grid_varnames:
+            ki = np.array([list(kp1_levels).index(x) for x in k_levels])
+            for zv,sl in zip(['Zl','Zu'],[ki,ki+1]):
+                variables[zv] = xr.Variable(_VAR_METADATA[zv]['dims'],
+                                            data['RF'][sl],
+                                            _VAR_METADATA[zv]['attrs'])
 
         ds = ds.update(variables)
+
+        if grid_vars_to_coords:
+            ds = ds.set_coords(gridlist)
+
         return ds
diff --git a/xmitgcm/llcreader/stores.py b/xmitgcm/llcreader/stores.py
index 6748cce..a8bff12 100644
--- a/xmitgcm/llcreader/stores.py
+++ b/xmitgcm/llcreader/stores.py
@@ -13,33 +13,51 @@ class BaseStore:
         Where to find the data within the filesystem
     shrunk : bool, optional
         Whether the data files have been tagged with `.shrunk`
-    mask_fs : fsspec.AbstractFileSystem, optional
-        Where to find the mask datasets to decode the compression
-    mask_path : str, optional
-        Path the the mask datasets on the ``mask_fs`` filesystem
+    mask_fs, grid_fs : fsspec.AbstractFileSystem, optional
+        Where to find the mask or grid datasets to decode the compression
+    mask_path, grid_path : str, optional
+        Path to the mask or grid datasets on the ``mask_fs`` or ``grid_fs`` filesystem
+    shrunk_grid : bool, optional
+        Whether the grid files have been tagged with `.shrunk`
+        not always the same as for data variables
     join_char : str or None
         Character to use to join paths. Falls back on os.path.join if None.
     """
 
     def __init__(self, fs, base_path='/', shrunk=False,
-                 mask_fs=None, mask_path=None, join_char=None):
+                 mask_fs=None, mask_path=None, 
+                 grid_fs=None, grid_path=None,
+                 shrunk_grid=False, join_char=None):
         self.base_path = base_path
         self.fs = fs
         self.shrunk = shrunk
         self.mask_fs = mask_fs or self.fs
         self.mask_path = mask_path
+        self.grid_fs = grid_fs or self.fs
+        self.grid_path = grid_path
+        self.shrunk_grid = shrunk_grid
         self.join_char = join_char
         if shrunk and (mask_path is None):
             raise ValueError("`mask_path` can't be None if `shrunk` is True")
 
 
     def _directory(self, varname, iternum):
-        return self.base_path
+        if iternum is not None:
+            return self.base_path
+        else:
+            return self.grid_path
 
     def _fname(self, varname, iternum):
-        fname = varname + '.%010d.data' % iternum
-        if self.shrunk:
-            fname += '.shrunk'
+
+        if iternum is not None:
+            fname = varname + '.%010d.data' % iternum
+            if self.shrunk:
+                fname += '.shrunk'
+        else:
+            fname = varname + '.data'
+            if self.shrunk_grid:
+                fname = varname + '.shrunk'
+
         return fname
 
     def _join(self, *args):
@@ -103,4 +121,7 @@ class NestedStore(BaseStore):
     iteration number."""
 
     def _directory(self, varname, iternum):
-        return self._join(self.base_path, '%010d' % iternum)
+        if iternum is not None:
+            return self._join(self.base_path, '%010d' % iternum)
+        else:
+            return self.grid_path
diff --git a/xmitgcm/mds_store.py b/xmitgcm/mds_store.py
index 8b6fa1b..bf4bb3e 100644
--- a/xmitgcm/mds_store.py
+++ b/xmitgcm/mds_store.py
@@ -763,7 +763,7 @@ def _guess_layers(data_dir):
     return all_layers
 
 
-def _get_all_grid_variables(geometry, grid_dir, layers={}):
+def _get_all_grid_variables(geometry, grid_dir=None, layers={}):
     """"Put all the relevant grid metadata into one big dictionary."""
     possible_hcoords = {'cartesian': horizontal_coordinates_cartesian,
                         'llc': horizontal_coordinates_llc,
@@ -772,7 +772,7 @@ def _get_all_grid_variables(geometry, grid_dir, layers={}):
     hcoords = possible_hcoords[geometry]
 
     # look for extra variables, if they exist in grid_dir
-    extravars = _get_extra_grid_variables(grid_dir)
+    extravars = _get_extra_grid_variables(grid_dir) if grid_dir is not None else {}
 
     allvars = [hcoords, vertical_coordinates, horizontal_grid_variables,
                vertical_grid_variables, volume_grid_variables, mask_variables,
diff --git a/xmitgcm/variables.py b/xmitgcm/variables.py
index 24f1691..56aa0b6 100644
--- a/xmitgcm/variables.py
+++ b/xmitgcm/variables.py
@@ -187,7 +187,7 @@ horizontal_grid_variables = OrderedDict(
     rAs=dict(dims=["j_g", "i"], attrs=dict(
         standard_name="cell_area_at_v_location",
         long_name="cell area", units="m2", coordinates="YG XC"),
-        filename='RAZ'),
+        filename='RAS'),
 )
 
 vertical_grid_variables = OrderedDict(
@@ -651,7 +651,31 @@ extra_grid_variables = OrderedDict(
     maskCtrlS=dict(dims=['k', 'j_g', 'i'], attrs=dict(
         standard_name="ctrl_vector_3d_mask_at_v_location",
         long_name='CTRL 3D mask where ctrl vector is active at v location',
-        units=''))
+        units='')),
+    # Reference density profile
+    rhoRef=dict(dims=['k'],attrs=dict(
+        standard_name="reference_density_profile",
+        long_name="1D, vertical reference density profile",
+        coordinate="Z",
+        units='kg m-3'),
+        filename='RhoRef'),
+    # Additional grid metrics
+    dxF=dict(dims=['j','i'],attrs=dict(
+        standard_name="cell_x_size_at_t_location",
+        long_name="cell x size", units="m", coordinate="YC XC"),
+        filename='DXF'),
+    dyF=dict(dims=['j','i'],attrs=dict(
+        standard_name="cell_y_size_at_t_location",
+        long_name="cell y size", units="m", coordinate="YC XC"),
+        filename='DYF'),
+    dxV=dict(dims=['j_g','i_g'],attrs=dict(
+        standard_name="cell_x_size_at_f_location",
+        long_name="cell x size", units="m", coordinate="YG XG"),
+        filename='DXV'),
+    dyU=dict(dims=['j_g','i_g'],attrs=dict(
+        standard_name="cell_y_size_at_f_location",
+        long_name="cell y size", units="m", coordinate="YG XG"),
+        filename='DYU')
     # Printed from write_grid when sigma coordinates are used
     # AHybSigmF, BHybSigF, ...
     # Unclear where on the grid these variables exist,


Test patch:
diff --git a/xmitgcm/test/test_llcreader.py b/xmitgcm/test/test_llcreader.py
index f2c1d16..c41e16b 100644
--- a/xmitgcm/test/test_llcreader.py
+++ b/xmitgcm/test/test_llcreader.py
@@ -1,4 +1,5 @@
 import pytest
+from dask.array.core import Array as dsa
 
 llcreader = pytest.importorskip("xmitgcm.llcreader")
 
@@ -8,6 +9,16 @@ EXPECTED_VARS = ['Eta', 'KPPhbl', 'oceFWflx', 'oceQnet', 'oceQsw', 'oceSflux',
             'oceTAUX', 'oceTAUY', 'PhiBot', 'Salt', 'SIarea', 'SIheff',
             'SIhsalt', 'SIhsnow', 'SIuice', 'SIvice', 'Theta', 'U', 'V', 'W']
 
+EXPECTED_COORDS = {2160: ['CS','SN','Depth',
+                          'drC','drF','dxC','dxF','dxG','dyC','dyF','dyG',
+                          'hFacC','hFacS','hFacW','PHrefC','PHrefF','rA','rAs','rAw',
+                          'Z','Zp1','Zl','Zu','rhoRef','rLowC','rLowS','rLowW',
+                          'rSurfC','rSurfS','rSurfW','XC','YC'],
+                   4320: ['CS','SN','Depth',
+                          'drC','drF','dxC','dxF','dxG','dyC','dyF','dyG',
+                          'hFacC','hFacS','hFacW','PHrefC','PHrefF',
+                          'rA','rAs','rAw','rhoRef','Z','Zp1','Zl','Zu','XC','YC']}
+
 ########### Generic llcreader tests on local data ##############################
 
 @pytest.fixture(scope='module')
@@ -28,30 +39,35 @@ def test_llc90_local_faces(local_llc90_store, llc90_kwargs):
     ds_faces = model.get_dataset(**llc90_kwargs)
     assert set(llc90_kwargs['varnames']) == set(ds_faces.data_vars)
     assert ds_faces.dims == {'face': 13, 'i': 90, 'i_g': 90, 'j': 90, 'j_g': 90,
-                             'k': 50, 'k_u': 50, 'k_l': 50, 'k_p1': 50, 'time': 2}
+                             'k': 50, 'k_u': 50, 'k_l': 50, 'k_p1': 51, 'time': 2}
 
 def test_llc90_local_latlon(local_llc90_store, llc90_kwargs):
     store = local_llc90_store
     model = llcreader.LLC90Model(store)
     ds_latlon = model.get_dataset(type='latlon', **llc90_kwargs)
     assert set(llc90_kwargs['varnames']) == set(ds_latlon.data_vars)
-    assert ds_latlon.dims == {'i': 360, 'time': 2, 'k_p1': 50, 'face': 13,
+    assert ds_latlon.dims == {'i': 360, 'time': 2, 'k_p1': 51, 'face': 13,
                               'i_g': 360, 'k_u': 50, 'k': 50, 'k_l': 50,
                               'j_g': 270, 'j': 270}
 
 @pytest.mark.parametrize('rettype', ['faces', 'latlon'])
-@pytest.mark.parametrize('k_levels', [None, [0, 2, 7, 9, 10, 20]])
+@pytest.mark.parametrize('k_levels, kp1_levels', 
+        [(None,None),
+         ([0, 2, 7, 9, 10, 20],
+          [0,1,2,3,7,8,9,10,11,20,21])])
 @pytest.mark.parametrize('k_chunksize', [1, 2])
 def test_llc90_local_faces_load(local_llc90_store, llc90_kwargs, rettype, k_levels,
-                                k_chunksize):
+                                kp1_levels, k_chunksize):
     store = local_llc90_store
     model = llcreader.LLC90Model(store)
     ds = model.get_dataset(k_levels=k_levels, k_chunksize=k_chunksize,
                            type=rettype, **llc90_kwargs)
     if k_levels is None:
         assert list(ds.k.values) == list(range(50))
+        assert list(ds.k_p1.values) == list(range(51))
     else:
         assert list(ds.k.values) == k_levels
+        assert list(ds.k_p1.values) == kp1_levels
     assert all([cs==k_chunksize for cs in ds['T'].data.chunks[1]])
 
     ds.load()
@@ -72,8 +88,15 @@ def test_ecco_portal_faces(ecco_portal_model):
     nx = ecco_portal_model.nx
     assert ds_faces.dims == {'face': 13, 'i': nx, 'i_g': nx, 'j': nx,
                               'j_g': nx, 'k': 90, 'k_u': 90, 'k_l': 90,
-                              'k_p1': 90, 'time': 3}
+                              'k_p1': 91, 'time': 3}
     assert set(EXPECTED_VARS) == set(ds_faces.data_vars)
+    assert set(EXPECTED_COORDS[nx]).issubset(set(ds_faces.coords))
+
+    # make sure vertical coordinates are in one single chunk
+    for fld in ds_faces[['Z','Zl','Zu','Zp1']].coords:
+        if isinstance(ds_faces[fld].data,dsa):
+            assert len(ds_faces[fld].data.chunks)==1
+            assert (len(ds_faces[fld]),)==ds_faces[fld].data.chunks[0]
 
 def test_ecco_portal_load(ecco_portal_model):
     # an expensive test because it actually loads data
@@ -88,6 +111,13 @@ def test_ecco_portal_latlon(ecco_portal_model):
     ds_ll = ecco_portal_model.get_dataset(iter_stop=iter_stop, type='latlon')
     nx = ecco_portal_model.nx
     assert ds_ll.dims == {'i': 4*nx, 'k_u': 90, 'k_l': 90, 'time': 3,
-                             'k': 90, 'j_g': 3*nx, 'i_g': 4*nx, 'k_p1': 90,
+                             'k': 90, 'j_g': 3*nx, 'i_g': 4*nx, 'k_p1': 91,
                              'j': 3*nx, 'face': 13}
     assert set(EXPECTED_VARS) == set(ds_ll.data_vars)
+    assert set(EXPECTED_COORDS[nx]).issubset(set(ds_ll.coords))
+
+    # make sure vertical coordinates are in one single chunk
+    for fld in ds_ll[['Z','Zl','Zu','Zp1']].coords:
+        if isinstance(ds_ll[fld].data,dsa):
+            assert len(ds_ll[fld].data.chunks)==1
+            assert (len(ds_ll[fld]),)==ds_ll[fld].data.chunks[0]