[Dart-dev] [7223] DART/trunk/models/POP/dart_pop_mod.f90: the written documentation for the details of the POP grid is inconsistent with

[nancy at ucar.edu]

Revision: 7223
Author:   nancy
Date:     2014-10-28 13:52:18 -0600 (Tue, 28 Oct 2014)
Log Message:
-----------
the written documentation for the details of the POP grid is inconsistent with 
the current POP code.  the KMT array is read from a grid file which stores
the lowest valid level in each column (0 = land).  an offset KMU array
is computed from that for the velocity points.  the original code was
using an offset in the wrong direction when trying to index the 4 T points
that enclosed each U point.  it also did not handle the wrap-around in
longitude correctly, leaving those KMU values undefined.

the impacts were:  possibly failed forward operators for U and V currents
if they were along the 'seam' in the shifted-pole grid; and current values 
in the ocean along the 'seam' not being impacted by observations when they should
have been.  the seam in the shifted grid runs along south america, up the
atlantic ocean, and into greenland.  the impacts are observable, but minimal.

the fix uses the right 4 KMT points to compute each KMU; it handles the
wraparound in longitude correctly; and it warns if the pole is unshifted
since we believe that the model itself might not correctly support a grid
with a singular point at the north pole.  this warning can be commented out
since the assimilation will work correctly with an unshifted grid; it is
more of a concern that the model advance will not have correctly computed
values at the pole.

Modified Paths:
--------------
    DART/trunk/models/POP/dart_pop_mod.f90

-------------- next part --------------
Modified: DART/trunk/models/POP/dart_pop_mod.f90
===================================================================
--- DART/trunk/models/POP/dart_pop_mod.f90	2014-10-27 21:35:28 UTC (rev 7222)
+++ DART/trunk/models/POP/dart_pop_mod.f90	2014-10-28 19:52:18 UTC (rev 7223)
@@ -635,13 +635,57 @@
 read( topo_unit, rec=1) KMT
 close(topo_unit)
 
-KMU(1, 1) = 0
-do j=2,ny
-do i=2,nx
-   KMU(i,j) = min(KMT(i, j), KMT(i-1, j), KMT(i, j-1), KMT(i-1, j-1))
+! the equation numbered 3.2 on page 15 of this document:
+!  http://www.cesm.ucar.edu/models/cesm1.0/pop2/doc/sci/POPRefManual.pdf
+! is WRONG.  (WRONG == inconsistent with the current POP source code.)
+!
+! for any U(i,j), the T(i,j) point with the same index values is located
+! south and west. so the T points which surround any U(i,j) point are
+! in fact at indices i,i+1, and j,j+1 .
+!
+!  NO: KMU(i,j) = min(KMT(i, j), KMT(i-1, j), KMT(i, j-1), KMT(i-1, j-1)) 
+! YES: KMU(i,j) = min(KMT(i, j), KMT(i+1, j), KMT(i, j+1), KMT(i+1, j+1))
+!
+! the latter matches the POP source code, on yellowstone, lines 908 and 909 in:
+!  /glade/p/cesm/releases/cesm1_2_2/models/ocn/pop2/source/grid.F90
+!
+! wrap around longitude boundary at i == nx.  set the topmost (last) latitude
+! U row to the same value in all cases. in the shifted pole grid currently in 
+! use all these points are on land and so are 0.  in the original unshifted
+! lat/lon grid these last row U points are above the final T row and are believed
+! to be unused.  for completeness we set all values in the last U row to the 
+! minimum of the all T row values immediately below it, for all longitudes.
+
+do j=1,ny-1
+   do i=1,nx-1
+      KMU(i,j) = min(KMT(i, j), KMT(i+1, j), KMT(i, j+1), KMT(i+1, j+1))
+   enddo
+   KMU(nx,j) = min(KMT(nx, j), KMT(1, j), KMT(nx, j+1), KMT(1, j+1))
 enddo
-enddo
+KMU(:,ny) = minval(KMT(:,ny))
 
+! IF YOU KNOW WHAT YOU ARE DOING YOU CAN COMMENT THIS TEST AND ERROR
+! CODE OUT, but keep reading to be sure:
+! the dart interface code can handle interpolation in any supported
+! POP grid, including a standard latitude/longitude grid.  but it is
+! not clear that the POP model code can compute values correctly if
+! the top row of the grid is in fact a singular point at the north pole.
+! a shifted grid puts this row over land where singularities can be ignored.
+! if you are sure your version of POP can compute the top grid row values
+! correctly, comment out the code between START HERE and END HERE and
+! recompile.  otherwise, use a shifted-pole grid. 
+
+! START HERE
+if (any(KMT(:,ny) /= 0)) then
+   msgstring = 'north boundary of grid does not appear to be over land'
+   call error_handler(E_ERR,'read_topography', &
+          msgstring, source, revision, revdate, &
+          text2='to continue, comment out this error call in the source code', &
+          text3='at end of read_topography() in models/POP/dart_pop_mod.f90')
+   
+endif
+! END HERE
+
 end subroutine read_topography