[MITgcm-support] MITgcm density anomaly

[John Pender]

Hi Jean-Michel-

I am only just becoming aware of parameters like staggerTimeStep so, no, I have not yet included this in my code.

Originally, the custom diagnostics were evaluated as a postprocessing operation using netcdf files imported into matlab.  One of them, for instance, uses
	UVEL
	RHOAnom
	grid
as written to file, where every file uses the same time grid (output on the half hour).  I would average UVEL(x,y,t) and UVEL(x+1,y,t) to get UVEL at the cell center and then calculate Uvel - Ubar for the cell centers, where Ubar is the weighted average velocity for the water column.  This seems to work fine, though I might need to take another look at the interface between seawater and the ocean floor.  This quantity is then multiplied times a pressure-related quantity which I would calculate from the density anomaly to my boss' satisfaction.

It got tricker when I was asked to calculate this diagnostic from within MITgcm because that's when I noticed that RHOAnoma wasn't written to file (when asked by input/data.diagnostics) at the same time as most of the other diagnostic output (UVEL, ETAN, etc).  And this appears to be the case because the new values for RHOAnoma are not calculated at the same time as the new values for UVEL and the other dynamical variables.  If I calculate my diagnostic in 
	code/diagnostics_fill_state.F
the numbers for UVEL will be up to date but the numbers for rhoInSitu will correspond to an earlier time.

I have written UVEL to file from 
	code/diagnostics_fill_state.F
and I have written UVEL to file from
	code/diags_rho.F
and the numbers are the same, so I could put my diagnostic code in diags_rho.f and write to file from there but it would seem that the numbers for rhoInSitu would then refer to a somewhat later time compared to the numbers for UVEL. 

So I guess what my boss would like me to figure out is this:  if I were to take the average of 
	the numbers for rhoInSitu that are in diagnostics_fill_state.F 
and 
	the numbers for rhoInSitu that are in diags_rho.F
would I get numbers for rhoInSitu that correspond to the same point in time as the numbers for UVEL? 

Thanks very much for your help on this.  
John



On Sep 14, 2012, at 10:43 AM, Jean-Michel Campin wrote:

> Hi John,
> 
> Sorry for not beeing clear (and generating confusion).
> The half time-step I mentionned was refering to the section 2.8 of 
> the mamual (e.g., figure 2.7), when a stagger timestepping is used.
> And in addition, I type this wrong:
>>>     or with UTHMASS,USLTMASS (computed using uVel^(n-1/2) if staggerTimeStep)
> since it should have been: (computed using uVel^(n+1/2) if staggerTimeStep)
> 
> And yes, your interpretation (a,b,c) is right when staggerTimeStep=T.
> But regarding the output file (NetCDF or *.meta files), the model 
> does not make this distinction (iteration number is just a plain
> integer, corresponding to when the file is written).
> 
> I am not sure about the new/modified pieces of code you added in
> code/diagnostics_fill_state.F
> Are you using staggerTimeStep ?
> 
> And if, when using staggerTimeStep=T., we really want to make a clean
> comparison between URHOMASS and  UTHMASS,USLTMASS, the best option would be
> to call DIAGS_RHO_G (within bi,bj loops) from diagnostics_fill_state.F,
> in the section 
>     IF ( selectVars.EQ.1 .OR. selectVars.EQ.3 ) THEN
> but only if staggerTimeStep=T., and put the call to DIAGS_RHO_G
> in do_oceanic_phys.F within IF ( .NOT. staggerTimeStep ) THEN / ENDIF
> 
> If you want to make this change and need help, let us know.
> 
> And just to finish: There are 2 reasons why it has not yet been done:
> a) historical, since rhoInSitu  used to be just a local array 
>   (not in common block)
> b) this would make the comparison between WRHOMASS and 
>  WdRHO_P,WdRHOdP less straitforward (and more tricky to fix since
>  these later 2 diags use rhoKm1 which is just a 2-D local array).
> 
> Cheers,
> Jean-Michel
> 
>>> URHOMASS with either UVELTH,UVELSLT (no hFac) 
>>>     or with UTHMASS,USLTMASS (computed using uVel^(n-1/2) if staggerTimeStep)
> On Fri, Sep 14, 2012 at 08:52:23AM -0800, John Pender wrote:
>> Hi Jean-Michel-
>> 
>> Thanks very much for getting back to me.  I'm a bit confused by the n+1/2 vs n-1/2 notation.  Does this suggest that some variables are calculated on the timestep and some variables are calculated on the half time step?  I can see how this would work:  
>> 
>> 	a)  use density and other parameters to calculate pressure gradients, which accelerate matter
>> 			gives new values to uvel, vvel, etc	
>> 			* on time step
>> 	b)  use the velocities to update density, etc (perhaps with a flux divergence)
>> 			* on half time step
>> 	c)  repeat
>> 
>> but when I do an ncdump on, for instance, UVEL and on RHOAnoma (which Martin says is well vetted) I see that the data appear to have been written for the same sets of times.
>> 
>> I have been tasked with calculating and writing diagnostic variables that use velocities (uvel and vvel) and the density anomaly.  I have already done this in
>> 	code/diagnostics_fill_state.F
>> but I am getting rumblings from above as to whether I shouldn't take more care.  If the rhoInSitu numbers in diagnostics_fill_state.F are a half step out of date then would it make sense for me to bin these numbers in something like rhoInSituOld, and then move the diagnostic code to DIAGS_RHO_G.  rhoInSitu would have been updated (picture me waving my arms around a bit here) to a time a little later than uvel and vvel, but if I take the average of rhoInSitu and rhoInSituOld then (hopefully) I would have a quantity that sits on the same timestep as the other parameters.
>> 
>> Thanks very much,
>> John
>> 
>> 
>> 
>> 
>> 
>> On Sep 13, 2012, at 2:41 PM, Jean-Michel Campin wrote:
>> 
>>> Hi John,
>>> 
>>> Just for clarification regarding this:
>>>> 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.
>>> 
>>> There should not be a 1 time-step offset between T/S and RHOAnom, since RHOAnom 
>>> is filled up from subroutine DIAGS_RHO_G, called from do_oceanic_phys.F just 
>>> after computing density from identical T/S fields than the one diagnosed in 
>>> diagnostics_fill_state.F. 
>>> 
>>> However, things are less clear regarding diagnostics of density transport
>>> [U,V,W]RHOMASS, because you can argue that when staggerTimeStep is used,
>>> it would be better to diagnose Rho^(n)*uVel^(n+1/2) (to be more consistent
>>> with T/S advection) instead of Rho^(n)*uVel^(n-1/2) as we currently do.
>>> And you should be careful when comparing
>>> URHOMASS with either UVELTH,UVELSLT (no hFac) 
>>>     or with UTHMASS,USLTMASS (computed using uVel^(n-1/2) if staggerTimeStep)
>>> 
>>> Cheers,
>>> Jean-Michel
>> 
>> 
>> _______________________________________________
>> MITgcm-support mailing list
>> MITgcm-support at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-support
> 
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support