Bug in ana_grid.h for BENCHMARK

[jpringle]

There is a bug in the definition of lonp,latp,lonr,latr,lonu,latu and lonv,latv in ana_grid.h for the benchmark case. It does not affect the solution, since the metric definitions are correct. However, it can make interpreting plots difficult, and will cause errors if this analytical grid is then processed with any of the nesting tools.

The orignal code is

Code: Select all

#if defined BENCHMARK
!
!  Spherical coordinates set-up.
!
      dx=Xsize/REAL(Lm(ng),r8)
      dy=Esize/REAL(Mm(ng),r8)
      spherical=.TRUE.
      DO j=Jmin,Jmax
        val1=-70.0_r8+dy*(REAL(j,r8)-0.5_r8)
        val2=-70.0_r8+dy*REAL(j,r8)
        DO i=Imin,Imax
          lonr(i,j)=dx*(REAL(i,r8)-0.5_r8)
          latr(i,j)=val1
          lonu(i,j)=dx*REAL(i,r8)
          lonp(i,j)=lonu(i,j)
          latu(i,j)=latr(i,j)
          lonv(i,j)=lonr(i,j)
          latv(i,j)=val2
          latp(i,j)=latv(i,j)
        END DO
      END DO

To confirm that this is wrong, look at the lat/lon information written to the history file. If the grid structure is G, run the following code:

Code: Select all

markersize=12;
color='r';
linewidth=1;plot(G.lon_psi,G.lat_psi,['o',color],'markersize',markersize)
hold on
plot(G.lon_psi,G.lat_psi,['-',color],'markersize',markersize)
plot(G.lon_psi',G.lat_psi',['-',color],'markersize',markersize)
plot(G.lon_rho,G.lat_rho,['s',color],'markersize',markersize)
plot(G.lon_v,G.lat_v,['^',color],'markersize',markersize)
plot(G.lon_u,G.lat_u,['>',color],'markersize',markersize)

Compare grid to that shown in https://www.myroms.org/wiki/index.php/N ... _Technique; you will see that the psi and velocity grids are shifted relative to rho grid. A correct version of the code is

Code: Select all

#if defined BENCHMARK
!
!  Spherical coordinates set-up.
!
      dx=Xsize/REAL(Lm(ng),r8)
      dy=Esize/REAL(Mm(ng),r8)
      spherical=.TRUE.
      !changed by JMP to define grids by Psi grid
      DO j=Jmin,Jmax
        DO i=Imin,Imax
          lonp(i,j)=dx*REAL(i,r8)+ 0.0_r8-dx      !start psi grid at 0 lonitude
          latp(i,j)=dy*REAL(j,r8)-70.0_r8-dy      !and at -70 latitude. The minus dx, dy compensates for i,j starting at 0 in this loop, where psi starts at 1,1
          lonr(i,j)=lonp(i,j) + 0.5_r8*dx         !these offsets are confusing as heck
          latr(i,j)=latp(i,j) + 0.5_r8*dy         !to make sense of them, consider what
          lonu(i,j)=lonp(i,j)                     !lat,lon shift exist for a constant i,j
          latu(i,j)=latp(i,j) + 0.5_r8*dy         !between variables...
          lonv(i,j)=lonp(i,j) + 0.5_r8*dx         !e.g. what is the spatial offset between
          latv(i,j)=latp(i,j)                     !rho(1,1) and psi(1,1)?
        END DO
      END DO

This error is not all that important since it does not affect the solution. But since people might use the Benchmark case to modify to get analytical spherical grids for their own use, it is worth fixing.

Cheers,
Jamie Pringle