[MITgcm-support] MITgcm density anomaly

[Martin Losch]

Hi John, 

I am returning this to the MITgcm-support list.

I don't normally use the UDIAG/SDIAG variables. What I do instead is define my diagnostics in diagnostics_main_init.F (there are plenty of examples of how to do that), so that I can give it a reasonable name, that I can understand. This does not answer your question, but provides a way to avoid it.

density: rho3d = rhoInSitu (it's just the local variable, passed as rhoInSitu)
I'd use the RHOAnom diagnostics, as that is the tested and approved version, and the diagnostic is written shortly after it has been computed. Having saig that, rhoInSitu does not change during the time step,  but diagnostics_state_fill is called from do_statevars_diags, and that is actually called at the beginning of the time forward_step, ie before do_oceanic_phys. That means that you get only the density of the previous time step, and that is inconsistent with T/S. So there is a time offset of one time step. Maybe Jean-Michel can confirm this.

Martin



On Sep 9, 2012, at 10:16 PM, John Pender wrote:

> Hi Martin-
> 
> I have a couple of quick questions regarding the diagnostic output referenced on this page:
> 
> 	http://mitgcm.org/public/r2_manual/latest/online_documents/node267.html
> 
> My first question has to do with the density anomaly.  I have been tasked with coding some diagnostic variables which include the cell-center density anomaly as part of the calculation.  The write-to-file call for most of the diagnostics (ETAN, UVEL, VVEL, WVEL, SDIAGx, UDIAGx, THETA, etc.) are made in
> 
> 	     MITgcm/pkg/diagnostics/diagnostics_fill_state.F
> 
> which I have copied into my local code directory.  When I look at the list of includes for this file
> 
> #include "SIZE.h"
> #include "EEPARAMS.h"
> #include "PARAMS.h"
> #include "GRID.h"
> #include "SURFACE.h"
> #include "DYNVARS.h"
> #include "NH_VARS.h"
> 
> most of which are in
> 	.../MITgcm/model/inc
> I see the obvious candidate for the density anomaly I should use lives in DYNVARS.h:     a variable called  rhoInSitu
> 
> pacman7 1090% grep rho DYN*
> C     rhoInSitu    :: In-Situ density anomaly [kg/m^3] at cell center level.
>      &                       rhoInSitu,
>       _RL  rhoInSitu(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
> 
> I can write this variable to file using UDIAG1 if I like.  All of this is fine, but I also have the option of writing the density anomaly to file as RHOAnoma.  The weird thing is that the write-to-file call for RHOAnom is not made from diagnostics_fill_state.F (as are ETAN, UVEL, VVEL, THETA, SDIAGx, UDIAGx, etc).  Instead, it's made from
> 	MITgcm/model/src/diags_rho.F
> which is itself called from
> 	MITgcm/model/src/so_oceanic_phys.F
> What's more the write-to-file calls don't use rhoInSitu, they use rho3d.  
> 
> None of this would bother me at all except that the actual numbers written to file in 
> 	RHOAnoma 
> 	using rho3d 
> 	from a CALL DIAGNOSTICS_FILL statement in MITgcm/model/src/diags_rho.F
> are not the same as the actual numbers written to file in
> 	UDIAG1    
> 	using rhoInSitu
> 	from a CALL DIAGNOSTICS_FILL statement in code/diagnostics_fill_state.F
> 
> As far as my diagnostics go I don't know as I have a choice as to what I must use for the cell-center density anomaly (rhoInSitu) but I admit RHOAnoma confuses me.  Is there a grid offset or is there a time offset or is something else going on?
> 
> 
> My second question concerns the suite of diagnostic outputs provided by MITgcm (SDIAGx and UDIAGx).  I was under the impression that these outputs were on the cell-center grid by default but there is a discrepancy in the netcdf files.  I have set up a small run with modifications to 
> 	code/diagnostics_fill_state.F
> An identical set of numbers (single layer) is written to each of
> 
> 	SDIAG1, SDIAG2, SDIAG3, SDIAG4, SDIAG5, SDIAG6, SDIAG7, SDIAG8, SDIAG9,  and SDIAG10
> 
> in a series of ten CALL DIAGNOSTICS_FILL( … ) commands, and another set of numbers (rhoInSitu, as it turns out) to each of
> 
> 	UDIAG1, UDIAG2, UDIAG3, UDIAG4, UDIAG5, UDIAG6, UDIAG7, UDIAG8, UDIAG9,  and UDIAG10
> 
> in another series of ten CALL DIAGNOSTICS_FILL( … ) commands.  If I look at, for instance, the UDIAG1 output I see (using ncdump for netdcf output) that the data is written to the center grid
> 
> 	...
> 	double UDIAG1(T, Zmd000100, Y, X) ;
> 	…
> 
> but if I look at the UDIAG5 output I see that the x coordinate is offset to the face grid.
> 
> 	...
> 	double UDIAG5(T, Zmd000100, Y, Xp1) ;
> 	…
> 
> The actual numbers written to file are the same for either case, so what is different is the presumed grid the data refer to.
> 
> I'm assuming that this is a bug and that I should presume the cell-center grid even if the netcdf files claim otherwise.
> 
> Thanks very much,
> John