|
| 1 | +# MITgcm Velocity Rotation on LLC/Cubed-Sphere Grids |
| 2 | + |
| 3 | +## The Problem |
| 4 | + |
| 5 | +On MITgcm LLC (Lat-Lon-Cap) and cubed-sphere grids, UVEL and VVEL are stored in **local grid coordinates**, not geographic (eastward/northward). Each face of the cubed sphere has its own local coordinate system, so the i-direction and j-direction point in different geographic directions on different faces. |
| 6 | + |
| 7 | +When visualized on a lon-lat map without rotation, this causes visible discontinuities at face boundaries — velocities appear to jump abruptly where one face meets another, because the local axes change orientation. |
| 8 | + |
| 9 | +Tracer fields (THETA, SALT, ETAN, etc.) are scalar values and are unaffected by grid orientation — they display correctly as-is. |
| 10 | + |
| 11 | +### Grid staggering (secondary issue) |
| 12 | + |
| 13 | +UVEL lives on the XG/YC grid (western cell face), VVEL lives on the XC/YG grid (southern cell face), while tracers live on XC/YC (cell centers). This creates a small spatial offset but is NOT the primary cause of the visible discontinuities. The dominant effect is the face rotation. |
| 14 | + |
| 15 | +## The Solution: CS/SN Rotation |
| 16 | + |
| 17 | +MITgcm can output two grid files: |
| 18 | +- `AngleCS.data` / `CS.data` — cosine of the angle between grid-i and geographic-east |
| 19 | +- `AngleSN.data` / `SN.data` — sine of the angle between grid-i and geographic-east |
| 20 | + |
| 21 | +With these, geographic velocities are computed as: |
| 22 | + |
| 23 | +``` |
| 24 | +u_east = CS * UVEL - SN * VVEL |
| 25 | +v_north = SN * UVEL + CS * VVEL |
| 26 | +``` |
| 27 | + |
| 28 | +This transforms the local (i,j) velocity components to proper (east, north) components. |
| 29 | + |
| 30 | +## Implementation Plan |
| 31 | + |
| 32 | +### 1. Load CS/SN during `mitgcm_open()` |
| 33 | + |
| 34 | +In `MitgcmStore`, add: |
| 35 | +```c |
| 36 | +float *angle_cs; /* AngleCS or CS [ny*nx], NULL if unavailable */ |
| 37 | +float *angle_sn; /* AngleSN or SN [ny*nx], NULL if unavailable */ |
| 38 | +``` |
| 39 | + |
| 40 | +During `mitgcm_open()`, try loading in order: |
| 41 | +1. `AngleCS.data` / `AngleSN.data` |
| 42 | +2. `CS.data` / `SN.data` |
| 43 | + |
| 44 | +Both must be present; if either is missing, set both to NULL (no rotation). |
| 45 | + |
| 46 | +### 2. Mark velocity variables |
| 47 | + |
| 48 | +In `MitgcmVarData`, add: |
| 49 | +```c |
| 50 | +int velocity_component; /* 0=none, 1=U-component, 2=V-component */ |
| 51 | +``` |
| 52 | + |
| 53 | +During `mitgcm_scan_variables()`, detect velocity fields by name: |
| 54 | +- U-component: names starting with "U" or containing "UVEL" |
| 55 | +- V-component: names starting with "V" or containing "VVEL" |
| 56 | + |
| 57 | +This heuristic covers UVEL, VVEL, and diagnostic velocity fields. |
| 58 | + |
| 59 | +### 3. Pair U/V variables |
| 60 | + |
| 61 | +When scanning variables, if both a U-component and V-component exist in the same diagnostic group (same prefix, same dimensions), link them: |
| 62 | +```c |
| 63 | +USVar *paired_velocity; /* Pointer to the other component */ |
| 64 | +``` |
| 65 | + |
| 66 | +### 4. Apply rotation during `mitgcm_read_slice()` |
| 67 | + |
| 68 | +After reading the raw slab, if `angle_cs` is available and the variable is a velocity component: |
| 69 | + |
| 70 | +```c |
| 71 | +if (store->angle_cs && vd->velocity_component != 0 && vd->paired_velocity) { |
| 72 | + /* Need to read the paired component too */ |
| 73 | + float *other = malloc(slab_size * sizeof(float)); |
| 74 | + /* read paired component at same time/depth */ |
| 75 | + |
| 76 | + if (vd->velocity_component == 1) { |
| 77 | + /* This is U, other is V */ |
| 78 | + for (i = 0; i < slab_size; i++) { |
| 79 | + float u = data[i], v = other[i]; |
| 80 | + data[i] = store->angle_cs[i] * u - store->angle_sn[i] * v; |
| 81 | + } |
| 82 | + } else { |
| 83 | + /* This is V, other is U */ |
| 84 | + for (i = 0; i < slab_size; i++) { |
| 85 | + float u = other[i], v = data[i]; |
| 86 | + data[i] = store->angle_sn[i] * u + store->angle_cs[i] * v; |
| 87 | + } |
| 88 | + } |
| 89 | + free(other); |
| 90 | +} |
| 91 | +``` |
| 92 | + |
| 93 | +### 5. Handle missing CS/SN gracefully |
| 94 | + |
| 95 | +If CS/SN files are not present, velocities display in raw local coordinates (current behavior). No error, no warning — the data is still valid, just not geographically oriented. |
| 96 | + |
| 97 | +Optionally print an info message: |
| 98 | +``` |
| 99 | +MITgcm: AngleCS/AngleSN not found, velocity fields shown in local grid coordinates |
| 100 | +``` |
| 101 | + |
| 102 | +### Considerations |
| 103 | + |
| 104 | +- **Performance**: rotation requires reading two fields per slice instead of one. For interactive use this should be fine (binary reads are fast). |
| 105 | +- **Timeseries**: the `mitgcm_read_timeseries()` function would also need to read both components at each timestep for rotated output. This doubles the I/O for velocity time series. |
| 106 | +- **Grid staggering**: for full correctness, UVEL should be interpolated from XG to XC before rotation. This is a C-grid averaging: `u_at_center[i] = 0.5 * (uvel[i] + uvel[i-1])` in the i-direction (with wrapping within each face). This is a secondary refinement. |
| 107 | +- **Variable naming**: rotated variables could be renamed to indicate geographic orientation (e.g., "UVEL_east", "VVEL_north"), or the original names kept with a visual indicator. |
0 commit comments