dimensions:
  imax = 128;
  jmax = 64;
  nv = 4;
variables:
  //coordinates of Z
  float lat(imax,jmax);
    lat:long_name = "latitude";
    lat:units = "degrees_north";
    lat:bounds = "lat_bnds";
  float lon(imax,jmax);
    lon:long_name = "longitude";
    lon:units = "degrees_east";
    lon:bounds = "lon_bnds";
  //coordinates of the boundaries of D
  float lat_bnds(imax,jmax,nv);
  float lon_bnds(imax,jmax,nv);

dimensions:

//You could index by i and j 
zimax = 
zjmax = 7028 //number of cells in y direction: 7028
uimax = 1117 //number of u edges in x direction (i+1)
ujmax = 7028 // number of u edges in y direction (j)
vimax = 1116 //number of v edges in x direction (i)
vjmax = 7029 // number of v edges in y direction (j+1)
dimax = 1117 // number of d vertices in x direction (i+1)
djmax = 7029 // number of d vertices in y direction (j+1)


eta_rho = 60 ;
xi_rho = 160 ;
eta_u = 60 ;
xi_u = 159 ;
eta_v = 59 ;
xi_v = 160 ;
eta_psi = 59 ;
xi_psi = 1116 //number of cells in x direction: 1116 ;

//Or you could index by a vector
D = 7851393 //(i+1) * (j+1)
Z = 7843248 // i * j 
U = 7850276 // (i+1)j 
V = 7844364 // i(j+1)  

variables:
  //Z-grid
  //coordinates of centres
  double zx(zimax, zjmax);
    zx:standard_name = "longitude"; 
    zx:units = "degrees_east";
  double zy(zimax, zjmax);
    zy:standard_name = "latitude";
    zy:units = "degrees_north";

  //U-grid
  //coordinates of u edges
  double ux(uimax, ujmax);
    ux:standard_name = "longitude"; 
    ux:units = "degrees_east";
  double uy(uimax, ujmax);
    uy:standard_name = "latitude";
    uy:units = "degrees_north";

  //V-grid
  //coordinates of v edges
  double vx(vimax, vjmax);
    vx:standard_name = "longitude"; 
    vx:units = "degrees_east";
  double vy(vimax, vjmax);
    vy:standard_name = "latitude";
    vy:units = "degrees_north";


  //D-grid
  //coordinates of d vertices
  double dx(dimax, djmax);
    dx:standard_name = "longitude"; 
    dx:units = "degrees_east";
  double dy(dimax, djmax);
    dy:standard_name = "latitude";
    dy:units = "degrees_north";


  //Constituents are defined at the Z grid
  double c(zimax, zjmax);
    c:long_name = "rhodamine";
    c:units= "1e-3"
    
  //bedlevel
  double dp(dimax, djmax);
    dp:standard_name = "bedrock_altitude";
    db:unit="metre";
    db:positive="up";

  //bathymetry 
  double dp2(dimax, djmax);
    dp2:standard_name = "bedrock_altitude";
    db2:unit="metre";
    db2:positive="down";