Stream: python-dev

Topic: pop tools DZU, DZT

Deepak Cherian (Apr 03 2020 at 17:01):

cc @Matt Long @Frank Bryan

So I had to do this to get it to work:

dzt = ... (matches exactly)
dzu = 4 point min of dzt

# this is the funny part
dzu = xr.where(
        fortran_zindex >= KMU,   # fortran_zindex is integer in range [1, 62]
        dzt, # when condition is true
        dzu, # when condition is false
    )

Is this sensible? i.e. setting DZU at KMU and deeper to DZT

Deepak Cherian (Apr 03 2020 at 18:07):

Frank,

With

(
    expected[["DZU", "DZT"]]
    .isel(nlon=slice(-10, None), nlat=slice(70), z_t=0)
    .to_array(dim="variable")
    .plot(col="variable")
)

I get
pasted image

which seems totally wrong. These points are all land so it doesn't matter (does matter for testing though).
Once I set DZU = DZT in this region; everything is equal!

Deepak Cherian (Apr 03 2020 at 18:07):

expected is the grid file you shared

Frank Bryan (Apr 05 2020 at 20:48):

DZU should be the minimum of the 4 surrounding DZT no matter where you are. I'm not sure what the backstory is here - are you trying to avoid reading from the file and computing online?

Frank Bryan (Apr 05 2020 at 20:48):

If it is helpful, here is the IDL code that generated DZT & DZU in the netCDF grid file. IDL is zero-based indexing

;;; Create 3D cell thickness and cell volume arrays
for k=0,ndep-1 do begin
DWORK2D(,) = dz(k)
indx = where(k EQ KMT-1)
if ( max(indx) GE 0 ) then begin
DWORK2D(indx)=DZBC(indx)
endif
DZT(,,k) = DWORK2D

 ;;; DZU=min(DZT(i,j),DZT(i+1,j),DZT(i,j+1),DZT(i+1,j+1))
 DZU(*,*,k) = DWORK2D
 indx = where(shift(DWORK2D,-1,0) LT DZU(*,*,k))
 DZU(indx,k) = DWORK2D(indx)
 indx = where(shift(DWORK2D,0,-1) LT DZU(*,*,k))
 DZU(indx,k) = DWORK2D(indx)
 indx = where(shift(DWORK2D,-1,-1) LT DZU(*,*,k))
 DZU(indx,k) = DWORK2D(indx)

Deepak Cherian (Apr 06 2020 at 14:15):

The backstory is that I am trying to calculate DZU and DZT in python to avoid the need to generate "augmented" grid files. I am using your augmented grid file as a test case. See: https://github.com/NCAR/pop-tools/pull/44/files .

I get equal results (my python code compared to your grid file) once I fix DZU at the surface. In this image: user_uploads/2/IwMX9KsyFEYsBV_5OkvRvXQk/pasted_image.png
DZU in the augmented grid file has values that are not the 4-point minimum of DZT (which is all 1000). Once I set DZU in this region to be 1000 everything works. So there is a bug somewhere.

Frank Bryan (Apr 06 2020 at 14:36):

Your result is correct. All DZT and DZU in the first layer should be 10m.

I'll look at the code I used to generate the augmented file (a mess of IDL and netCDF operators) to see where the problem arises.
I know that there will be a problem along the tripole seam across the top of the grid. The neighbors across the seam are on the same row on the opposite side in the i-direction. I'll try to write out the proper operator in index notation . We might be able to work this into the differential operators themselves.

Frank Bryan (Apr 06 2020 at 19:48):

@Deepak Cherian I indeed see the spurious values in my netCDF file. There are 46 values in a cluster on the right edge of the grid (i=3599, j=16:61,k=0) that are bad. Somehow, they are getting munged when I write them out to netCDF. They are correct in the loop that computes them.

Looking at how I might fix this.

Deepak Cherian (Apr 07 2020 at 15:49):

@Frank Bryan thanks. FYI this error isn't blocking me right now so there's no urgency.

Frank Bryan (Apr 07 2020 at 15:52):

determined that something weird was going on in the IDL shift operations that was overwriting some values. I was able to get correct results with a fully unrolled loop and array index notation. I'll replace the augmented grid file later today.
[more questions about TPOS MITgcm coming your way]

Last updated: May 16 2025 at 17:14 UTC