bug in AngleCS and AngleSN in curvilinear grid

[stanislav-martyanov]

Hello!

I think there is a bug in AngleCS and AngleSN in MITgcm. Actually, it concerns the underlying 2D fields angleCosC and angleSinC calculated in s/r CALC_GRID_ANGLES. This s/r is called from ini_curvilinear_grid while the 2D arrays angleCosC and angleSinC are assigned to 1 and 0 in ini_grid, respectively (it corresponds to a non-rotated grid). At this point everything is OK.

But somewhere in CALC_GRID_ANGLES there is a problem, because for grids I used for MITgcm, the angleCosC and angleSinC values of the easternmost and northernmost edges of the grid are not correct. There should not be an abrupt change in corresponding values because the grid is rather smooth everywhere.

By the way, the erroneous edge values of angleCosC also appear when a simple geographical grid is used (with the curvilinear option enabled).

According to model's code, this problem mainly affects the wind field interpolation at those edge points and some other less frequently-used features.

I've saved the output files with AngleCS and AngleSN values so anyone can verify the problem. The matlab script to read them is:

fin = fopen('AngleCS.data','rb','ieee-be');
angleCS=fread(fin,[544 518],'float32');
fin2 = fopen('AngleSN.data','rb','ieee-be');
angleSN=fread(fin2,[544 518],'float32');

Files:
angles.zip

Kind regards,
Stanislav Martyanov

[mjlosch]

Thanks for submitting this bug report.
The computation of AngleCS and AngleSN in s/r CALC_GRID_ANGLES relies entirely on the grid information contained inYG, dxG, and dyG, so it would be extremely helpful for debugging this, if you could provide a (simple, small, maybe not with dimensions (544,518)) setup (i.e. code and input directory) to reproduce the problem. For that you could use your setup, coarse grain, and remove any bathymetry etc, to keep just as much to have the angles computed and saved to disk.

[stanislav-martyanov]

Thanks for submitting this bug report. The computation of AngleCS and AngleSN in s/r CALC_GRID_ANGLES relies entirely on the grid information contained inYG, dxG, and dyG, so it would be extremely helpful for debugging this, if you could provide a (simple, small, maybe not with dimensions (544,518)) setup (i.e. code and input directory) to reproduce the problem. For that you could use your setup, coarse grain, and remove any bathymetry etc, to keep just as much to have the angles computed and saved to disk.

Hello! Here are compressed files with CPP options and 'data' parameters needed to run the model in a simple configuration. I tried to attach other files together with the initial conditions for T and S, and also bathymetry but the file-size limitations here on Github did not allow this. Still, the initial conditions are not really necessary for this issue, and default MITgcm settings can be used for a simple run.

To run it you only need to turn off all the unnecessary packages in the data.pkg and packages.conf (actually, almost all of them). The model will abort shortly because it lacks the boundary files, forcings, etc, but it will (I hope) generate the output grid files.

I forgot to mention that the MITgcm versions I used were 65z (this config: a very old and currently obsolete configuration) and recently 67x. The above-mentioned problem with angles persists in both configurations so I believe that it is not a matter of a version.

Regards,
Stan

Files:
kar_bar_settings2.zip

[mjlosch]

Thanks, I can reproduce the problem: The last row (j=Ny) and the last column (i=Nx) of AngleCS/SN are not correctly computed in this regional domain.
I think that this is the reason: we specify all input grid fields from 1:Nx and 1:Ny. In a periodic domain (including cubed sphere / llc) all overlaps are filled correctly by the exchanges before computing the angles. In a non-periodic domain the coordinates (and dx/yG) are not correct for Nx+1 to Nx+Olx and Ny+1 to Ny+Oly, but the values in Nx+1 and Ny+1 are needed for computing the angles at Nx and Ny.

AngleCS/SN are never used in the code, but are only provided as a convenience output. Alternatively one can compute these angle (correctly) offline and provide them to the model (just not with OLD_GRID_IO). Further: a non periodic domain should either have land or open boundaries at Nx and Ny, which again reduces the problem.

I agree, it's a bug, but not a severe one. We could issue a warning at runtime, if a non-period domain is used. Or we could try to fix the values of YG at Nx/y+1 along the boundaries (not clear how).

[mjlosch]

suggestion: after l129 in ini_curvilinear_grid.F

      CALL EXCH_Z_3D_RS( yG, 1, myThid )

insert

      IF (.not.useCubedSphereExchange) THEN
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         j = sNy
         DO i = 1-Olx, sNx+Olx
          IF ( ABS(yG(i,j+1,bi,bj) - yG(i,j,bi,bj)) .GT.
     &         2.*ABS(yG(i,j,bi,bj) - yG(i,j-1,bi,bj)) ) THEN
           yG(i,j+1,bi,bj) = 2.*yG(i,j,bi,bj) - yG(i,j-1,bi,bj)
          ENDIF
         ENDDO
         i = sNx
         DO j = 1-Oly, sNy+Oly
          IF ( ABS(yG(i+1,j,bi,bj) - yG(i,j,bi,bj)) .GT.
     &         2.*ABS(yG(i,j,bi,bj) - yG(i-1,j,bi,bj)) ) THEN
           yG(i+1,j,bi,bj) = 2.*yG(i,j,bi,bj) - yG(i-1,j,bi,bj)
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDIF

to extrapolate YG if there is a jump between Nx/y and Nx/y+1. Similarly in ini_spherical_polar_grid.F

@stanislav-martyanov maybe you would like to try that.

[jm-c]

With cartesian or SphericalPolarGrid, the grid setting is simple, only few user option and it's easy to set the grid correctly regarding the connectivity (doublely periodid by default). But with curvilinear-grid, the user need to provide all the grid metrics and it's less easy to get a consistent grid-setting regarding the connectivity. And the code does not check for any inconsistency (that depends on the use of OBC and also land barriers).

@stanislav-martyanov in your case, as I understand, the grid-metrics are not supposed to be periodic in X and Y, and what you found is that the grid-angles are not right at the last row/column but likely other quantities are not right neither and the last row and column should not be part of the domain (either land or outside Open-Boundary).
So the main question would be, does the model needs these angles coeff. at the last row and column of the domain when this part of the domain is not computed (either land or outside OB).

[jm-c]

Also some discussion relevant to similar questions in issue #362

[stanislav-martyanov]

AngleCS/SN are never used in the code, but are only provided as a convenience output.

Yes, AngleCS/SN are not used anywhere in the model (these are just the output names used in WRITE_GRID), but, as I understand, the real fields they represent (angleCosC and angleSinC) are (sometimes). As the search tells, these 2D fields can be used in:

1. rotate_uv2en
2. GMREDI_K3D, GMREDI_INIT_VARIA
3. EXF_SET_UV
4. MOM_W_CORIOLIS_NH, MOM_QUASIHYDROSTATIC, MOM_U_CORIOLIS_NH
5. cost_gencost_customize
6. DIAGNOSTICS_FILL_STATE

Though I must admit I did not deeply dig all these subroutines. It looks like the most important is EXF_SET_UV because it interpolates and rotates the input wind field on-the-fly and is used almost always in realistic configurations. But I did not investigate yet if these interpolated wind values at the edges are actually used in calculations.

So the main question would be, does the model needs these angles coeff. at the last row and column of the domain when this part of the domain is not computed (either land or outside OB).

I believe it does not but I'm not absolutely sure because of the aforementioned routines I'm not fully familiar with. Anyway, I agree with mjlosch that the issue is not severe.

To solve this issue we've just slightly modified our post-processing tools without any modifications to the MITgcm code itself, so now it is not a problem.

--
This issue appeared absolutely accidentally when we compared the computed velocity field at the vicinity of domain's edges with the imposed external velocity data (OBCS). For that task we were to rotate the vectors and decided to use AngleCS/SN generated by the model, so the problem with AngleCS/SN values revealed itself.

So, if anybody decides to do anything similar in the future, it would be useful for him/her to know about this issue to avoid a surprise and a headache :-)

Anyway, it is up to the team to decide if any correction/warning/doc-info should be added.

Regards,
Stan

[mjlosch]

@stanislav-martyanov you are right, the angles are indeed used (sorry about my previous, wrong statement) and they have to be correct on the computational domain.

I have the impression that in your configuration this is not a problem because the row/column in question are an open boundary where velocity and tracer values are prescribed and nothing is computed, so that the problem only appears in post processing if you are trying access values near the boundary. Most likely most configurations will be like this (either open boundary conditions or wall), but I can imagine configurations, either with a curvilinear grid, or an idealized sector model with walls, etc., where the last row/columns are not open boundaries and the wall has been put at i=1 or j =1.

So I think some sort of check is in place if the domain is periodic. I am not sure if it is worth thinking about an extrapolation scheme for YG as outlined above (I am sure that I can be done better than my suggestion).

[jm-c]

Regarding this comment (from @mjlosch):

I am not sure if it is worth thinking about an extrapolation scheme for YG as outlined above (I am sure that I can be done better than my suggestion).

I think I answered this point on Apr 22:

as I understand, the grid-metrics are not supposed to be periodic in X and Y, and what you found is that the grid-angles are not right at the last row/column but likely other quantities are not right neither and the last row and column should not be part of the domain (either land or outside Open-Boundary).

I suggest to close this issue, because (1) it has not been active for few weeks ; (2) not clear it the reported problem is when using OBC (which makes this issue less useful to keep it open) ; (3) if the set-up does not use OBC then this discussion could continue in issue #362.