diagnostics
-
- Posts: 14
- Joined: Wed Jun 30, 2010 2:58 pm
- Location: Texas A&M University
diagnostics
Hi,
I'm trying to do some momentum balance analysis using ROMS. I output the 2D momentum terms into the diagnostic output files, but I'm not quite sure how the momentum balance is completed. I chose a point in the model, and added up all the terms (ubar_sustr+ubar_bustr+ubar_prsgrd+ubar_cor+ubar_hadv+ubar_accel+ubar_hvisc) in the xi-momentum equation at this point, but the result is not zero. Could anyone give some explanations on this? Thanks!
I'm trying to do some momentum balance analysis using ROMS. I output the 2D momentum terms into the diagnostic output files, but I'm not quite sure how the momentum balance is completed. I chose a point in the model, and added up all the terms (ubar_sustr+ubar_bustr+ubar_prsgrd+ubar_cor+ubar_hadv+ubar_accel+ubar_hvisc) in the xi-momentum equation at this point, but the result is not zero. Could anyone give some explanations on this? Thanks!
Re: diagnostics
I have not worked extensively with the momentum diagnostics, and certainly not the depth averaged terms. But with the tracer diagnostics the sign convention can be interpreted as the time rate of change being on the left hand side, and all other terms on the right hand side.
So in the case of salt we expect:
salt_rate = salt_hadv + salt_vadv + salt_hdiff + salt_vdiff
So if you are used to thinking of the two advection terms as div(u.salt) - because we usually write it on the left hand side - this is:
salt_rate = -div_h(u.salt) + salt_hdiff + salt_vdiff
Therefore my expectation would be that you get this balance:
ubar_accel = ubar_sustr + ubar_bustr + ubar_prsgrd + ubar_cor + ubar_hadv + ubar_hvisc
You might want to check that the sign of the sustr and bustr terms is consistent with the sense of ubar.
So in the case of salt we expect:
salt_rate = salt_hadv + salt_vadv + salt_hdiff + salt_vdiff
So if you are used to thinking of the two advection terms as div(u.salt) - because we usually write it on the left hand side - this is:
salt_rate = -div_h(u.salt) + salt_hdiff + salt_vdiff
Therefore my expectation would be that you get this balance:
ubar_accel = ubar_sustr + ubar_bustr + ubar_prsgrd + ubar_cor + ubar_hadv + ubar_hvisc
You might want to check that the sign of the sustr and bustr terms is consistent with the sense of ubar.
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu
-
- Posts: 14
- Joined: Wed Jun 30, 2010 2:58 pm
- Location: Texas A&M University
Re: diagnostics
Hi Dr. Wilkins, yes your suggestion works! This time I computed 'ubar_sustr+ubar_bustr+ubar_prsgrd+ubar_cor+ubar_hadv+ubar_hvisc-ubar_accel' and the result is zero. Thank you!
Re: diagnostics
I'm using the diagnostic feature for some budgets. In the course of some error checking I compared the field from salt_rate to the difference between salt fields from the history files saved at the beginning and end of the diagnostic averaging period. The shapes and magnitudes, e.g. of a vertical profile of salt_rate compared to (salt2-salt1)/3600, look very similar, BUT there is an OFFSET that is as big as the signal. Any clues? The sum of all salt diagnostic terms balances perfectly, as expected.
Re: diagnostics
I did some more analysis on the problem of the salt_rate diagnostic. The plot shows the evolution of temperature and salt at a surface model location, both from history and average files. They look very similar. Then I calculate the time derivative of both and plot it on the bottom row, comparing it with the associated rate term from the diagnostics. As you can see, the salt-rate looks pretty bad. This pattern is similar at other model locations. Where is salt_rate calculated in the code?
![Image](http://pogo.ocean.washington.edu/parker/roms/rateTest.gif)
![Image](http://pogo.ocean.washington.edu/parker/roms/rateTest.gif)
Re: diagnostics
I wonder if the problem might arise in the two different versions of the calculation of DiaTwrk(i,j,k,itrc,iTrate) taken from step3d_t.F. I assume the first, active, version forces the salt_rate diagnostic to equal the sum of all the other terms. Then the second, commented out, version appears to calculate the salt_rate based on the difference of subsequent salt fields...
# ifdef DIAGNOSTICS_TS
!
! Compute time-rate-of-change diagnostic term.
!
DO k=1,N(ng)
DO j=JstrR,JendR
DO i=IstrR,IendR
DiaTwrk(i,j,k,itrc,iTrate)=t(i,j,k,nnew,itrc)- &
& DiaTwrk(i,j,k,itrc,iTrate)[/color]
!! DiaTwrk(i,j,k,itrc,iTrate)=t(i,j,k,nnew,itrc)- &
!! & t(i,j,k,nstp,itrc)
END DO
END DO
END DO
# endif
# ifdef DIAGNOSTICS_TS
!
! Compute time-rate-of-change diagnostic term.
!
DO k=1,N(ng)
DO j=JstrR,JendR
DO i=IstrR,IendR
DiaTwrk(i,j,k,itrc,iTrate)=t(i,j,k,nnew,itrc)- &
& DiaTwrk(i,j,k,itrc,iTrate)[/color]
!! DiaTwrk(i,j,k,itrc,iTrate)=t(i,j,k,nnew,itrc)- &
!! & t(i,j,k,nstp,itrc)
END DO
END DO
END DO
# endif
Re: diagnostics
I figured out the problem with my salt_rate test. salt_rate is NOT <ds/dt>, where <> denotes averaging over the time period for which diagnostics are collected. Instead (I think) salt_rate = <d(s*v)/dt>/<v> where v is the (changing) volume of a given grid cell. When I use this interpretation to check salt_rate using values from the avg and his fields it works perfectly.
- arango
- Site Admin
- Posts: 1356
- Joined: Wed Feb 26, 2003 4:41 pm
- Location: DMCS, Rutgers University
- Contact:
Re: diagnostics
Yes, ROMS is formulated in flux form. See the governing equations in
WikiROMS. All the diagnostics are also in flux form. Notice that the time-averaged file has also all the fluxes and quadratic terms. This becomes very handy when doing all these analysis. We just set the time averaging window for both the diagnostic (NDIA) and average (NAVG) files to be the same.
![Arrow :arrow:](./images/smilies/icon_arrow.gif)
- donoso
- Posts: 13
- Joined: Mon Jan 05, 2015 4:59 pm
- Location: Department of Geophysics. University of Concepcion
Re: diagnostics
Hello everyone,
I'm also interested in testing my tracer time rate of change. In particular, I'm trying to contrast the temperature time rate of change with the time derivative of temperature, after setting the time averaging window for both the diagnostic (NDIA) and average (NAVG) files to be the same. I followed the tracer equation flux form which would be:
temp_rate = <d(temp*v)/dt>/<v>
There, both "temp" and "v" are functions of time. So, I expanded them using the following Reynolds decomposition:
temp_rate = <d((<temp>+temp')*(<v>+v'))>/dt>/<v>
Which, if my math is correct, would lead to:
temp_rate = d(<temp>*<v>+<temp'*v'>)/dt/<v>
There, I think the quadratic term <temp'*v'> is missing to evaluate the time derivative because, currently, it is not possible to save that term into the average files. Am I right? If so, is there any other way to test the correspondence between temperature time rate of change and temperature? e.g. time integration of the first one would be easier?
Thanks in advance!
Cheers,
David.
I'm also interested in testing my tracer time rate of change. In particular, I'm trying to contrast the temperature time rate of change with the time derivative of temperature, after setting the time averaging window for both the diagnostic (NDIA) and average (NAVG) files to be the same. I followed the tracer equation flux form which would be:
temp_rate = <d(temp*v)/dt>/<v>
There, both "temp" and "v" are functions of time. So, I expanded them using the following Reynolds decomposition:
temp_rate = <d((<temp>+temp')*(<v>+v'))>/dt>/<v>
Which, if my math is correct, would lead to:
temp_rate = d(<temp>*<v>+<temp'*v'>)/dt/<v>
There, I think the quadratic term <temp'*v'> is missing to evaluate the time derivative because, currently, it is not possible to save that term into the average files. Am I right? If so, is there any other way to test the correspondence between temperature time rate of change and temperature? e.g. time integration of the first one would be easier?
Thanks in advance!
Cheers,
David.
Re: diagnostics
I think you can use: temp_rate = <d(temp*v)/dt>/<v> = <v>*<d(temp)/dt>/<v> + <temp>*<d(v)/dt>/<v> = <d(temp)/dt> + <temp>*<d(v)/dt>/<v>
and then rearrange to find:
<d(temp)/dt> = temp_rate - <temp>*<d(v)/dt>/<v>
You can calculate the LHS using the difference of temp from two history files (the ones surrounding the <> averaging window). And you can do the same for <d(v)/dt>. On the RHS you use <temp> from the AVERAGES, and I think you calculate <v> using <zeta> from the averages.
and then rearrange to find:
<d(temp)/dt> = temp_rate - <temp>*<d(v)/dt>/<v>
You can calculate the LHS using the difference of temp from two history files (the ones surrounding the <> averaging window). And you can do the same for <d(v)/dt>. On the RHS you use <temp> from the AVERAGES, and I think you calculate <v> using <zeta> from the averages.
- donoso
- Posts: 13
- Joined: Mon Jan 05, 2015 4:59 pm
- Location: Department of Geophysics. University of Concepcion
Re: diagnostics
Thank you Parker!
I will try your idea using some history files.
Just like you thought, I used <temp> and I calculated <v> (using <zeta>) both from the averages files. That is because I have a big domain, tidally forced and I don't have enough space storage to save a year of hourly history files. I wonder if, in general, it is not a good idea to derive averaged fields? I thought that an "average of a time derivative" was equal to a "time derivative of an average", which would lead your last equation be equal to:
d(<temp>)/dt = temp_rate - <temp>*d(<v>)/dt>/<v>
Where, both <temp> and <v> could be taken from the averages files. So, just to be sure, that equality is not correct math?
Cheers,
David.
![Smile :-)](./images/smilies/icon_smile.gif)
Just like you thought, I used <temp> and I calculated <v> (using <zeta>) both from the averages files. That is because I have a big domain, tidally forced and I don't have enough space storage to save a year of hourly history files. I wonder if, in general, it is not a good idea to derive averaged fields? I thought that an "average of a time derivative" was equal to a "time derivative of an average", which would lead your last equation be equal to:
d(<temp>)/dt = temp_rate - <temp>*d(<v>)/dt>/<v>
Where, both <temp> and <v> could be taken from the averages files. So, just to be sure, that equality is not correct math?
Cheers,
David.