Offline Heat Budget calculations from ROMS model output

[aneeshcs]

Hi,

Has anyone computed the individual terms in the heat budget equation (as shown below) using the ROMS model output, computing the advection, temperature tendency term and the remaining terms? If so, is it possible to completely close this budget to good accuracy or what is the order of magnitude of the residual terms in the budget?

heateqn.jpg (43.49 KiB) Viewed 25133 times

Thanks

[wilkin]

The DIAGNOSTICS_TS and AVERAGES options allow you to save everything you need, but it will take some care to get the balance you seek.

#define DIAGNOSTICS_TS will save (to the diagnostics file) the vertical advection, vertical mixing, lateral mixing, and time rate of changes terms provided you also activate the logical switches in ocean.in:

Dout(iTrate) == T F ! temp_rate, ... time rate of change
Dout(iThadv) == T F ! temp_hadv, ... horizontal total advection
Dout(iTvadv) == T F ! temp_vadv, ... vertical advection
Dout(iThdif) == T F ! temp_hdiff, ... horizontal total diffusion
Dout(iTvdif) == T F ! temp_vdiff, ... vertical diffusion

#define AVERAGES will save (to the averages file) the averages of heat flux, velocity <U>, temperature <T>, and quadratic product <UT> provided you also activate the logical switches in ocean.in:

Aout(idTsur) == T F ! shflux, ssflux surface net heat and salt flux
Aout(idUvel) == T ! u 3D U-velocity
Aout(idVvel) == T ! v 3D V-velocity
Aout(idTvar) == T F ! temp, salt temperature and salinity
Aout(idUTav) == T F ! utemp, usalt quadratic <u*t> T/S terms
Aout(idVTav) == T F ! vtemp, vsalt quadratic <v*t> T/S terms

You will need to construct <u't'> yourself from the saved <u*t> - <u>*<t>

Obviously, you must choose the same averaging time interval for averages and diagnostics, i.e. NAVG = NDIA.

Something to be aware of is that because of the time-changing vertical s coordinate your vertical integrals may be only approximate. ROMS offers the options to save the products of <dz*u*t> etc with these logical switches:

Aout(iHUTav) == T F ! Huontemp, ... T/S volume flux, <Huon*t>
Aout(iHVTav) == T F ! Hvomtemp, ... T/S volume flux, <Hvom*t>

These values then already have the vertical grid metrics included, so a vertical integral is achieved simply by a vertical summation. The horizontal grid metric dx = 1/m and dy = 1/n are also in these products ("Huon" means vertical spacing "H" times u "over n") so some care is needed here to use these values correctly.

This might not be an issue for you, but it was for me in my calculations of heat budgets on the New England Shelf: Wilkin, J. (2006):, The summer time heat budget and circulation of southeast New England shelf waters, Journal of Physical Oceanography, 36(11), 1997-2011 http://marine.rutgers.edu/%7Ewilkin/Wil ... PO2006.pdf I could not get an exact balance unless I took the time correlations of vertical grid spacing (which depends on sea level, zeta) and tidal velocity into account.

--
John.

[aneeshcs]

Thanks very much, John, for the detailed reply! Will do the online diagnostics for my model run.

[mashinde]

Hi,

I am also interested in heat budget. Which quantities do I need to use from dia* file and which formula i have to use ?

Thanks

[WencaiZhou]

Dear Prof. Wilikin and everyone,

I'm trying to analyze information from both, the AVERAGE file and the DIAGNOSTICS file for performing a salt budget, but there are some features that are a bit confusing. The first one I want to ask here is about the time averaging of the diagnostics variables.

All the internal documentation in the *.F files refers to the diagnostic variables as "time-averaged diagnostics". For example, the output.F says:
! If appropriate, process time-averaged diagnostics NetCDF file

and the wrt_diags.F file says:
!=======================================================================
! !
! This subroutine writes model time-averaged diagnostic fields into !
! diagnostics NetCDF file. !
! !
!=======================================================================

But I have not been able to find any section in any *.F file of the code that performs any time averaging on the diagonostics variables.
So, my question is:

1. Where is the time-averaging of the diagnostics variables performed? Or, why are these variables and files named as "time-averaged diagnostic" everywhere? Why are they called like that, what is the meaning?

Many thanks,
Wencai Zhou.

[wilkin]

The diagnostics calculations are throughout the code. Search for the CPP directive DIAGNOSTICS_TS (or DIAGNOSTICS_UV) and you will find snippets of code everywhere. For example ...

Code: Select all

queequeg:ROMS wilkin$ grep DIAGNOSTICS_TS Nonlinear/* | more
grep: Nonlinear/Biology: Is a directory
grep: Nonlinear/Sediment: Is a directory
Nonlinear/pre_step3d.F:# ifdef DIAGNOSTICS_TS
Nonlinear/pre_step3d.F:# ifdef DIAGNOSTICS_TS
Nonlinear/pre_step3d.F:#  ifdef DIAGNOSTICS_TS
Nonlinear/pre_step3d.F:#  ifdef DIAGNOSTICS_TS
Nonlinear/pre_step3d.F:# if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_UV
Nonlinear/pre_step3d.F:#  ifdef DIAGNOSTICS_TS
Nonlinear/step3d_t.F:# ifdef DIAGNOSTICS_TS
Nonlinear/step3d_t.F:# ifdef DIAGNOSTICS_TS
Nonlinear/step3d_t.F:# ifdef DIAGNOSTICS_TS
...

The time summation is performed by code such as (this taken from line 944 in step3d_t.F):

Code: Select all

#  ifdef DIAGNOSTICS_TS
              DiaTwrk(i,j,k,itrc,iTxadv)=DiaTwrk(i,j,k,itrc,iTxadv)-    &
     &                                   cff1
              DiaTwrk(i,j,k,itrc,iTyadv)=DiaTwrk(i,j,k,itrc,iTyadv)-    &
     &                                   cff2
              DiaTwrk(i,j,k,itrc,iThadv)=DiaTwrk(i,j,k,itrc,iThadv)-    &
     &                                   cff3
#  endif

The accumulated sums are converted to time averages in Utility/set_diags.F. Look for the code following the comment:

Code: Select all

!
!-----------------------------------------------------------------------
!  Convert accumulated sums into time-averages, if appropriate
!-----------------------------------------------------------------------
!

The length of the averaging interval is set in ocean.in by the keywords:

Code: Select all

! NDIA         Number of time-steps between writing time-averaged diagnostics
!                data into DIAGNOSTIC file.  Averaged date is written for all
!                fields. Set NDIA=0 to suppress writing of DIAGNOSTIC file.
!
! NDEFDIA      Number of time-steps between the creation of new time-averaged
!                diagnostics file.  If NDEFDIA=0, the model will only process
!                one DIAGNOSTICS file.  This feature is useful for long
!                simulations when DIAGNOSTIC file get too large; it creates
!                a new file every NDEFDIA time-steps.

[WencaiZhou]

Thank you very much Prof. Wilkin, now I understand it. Thanks for opening my eyes, now I see that all I asked for is in the set_diags.F and that the diagnostics variables are averaged in the same way as the average variables (in the set_avg.F).

Best regards,
Wencai

[WencaiZhou]

Dear Professor Wilkin,

I am wondering if you can help me further with another mysterious issue with the time-averaged variables. Basically, I'm not observing correspondence between variables in the history file and variables in the avg file. For example, salinity values of the same cell look unexpectedly different. I have posed this issue in a different thread:

viewtopic.php?f=14&t=4610

Could you please comment your opinion in that thread? Thanks a lot.

Wencai

[kunal.c]

Hi,

How do we get the values of K_h and K_v mentioned in the diffusivity term of the heat budget equation (the same equation which started the discussion and also attached for easy referral).
I suppose density can be derived by adding rho_anomaly with mean density. Please comment if I am wrong.
We take the value of C_p (=4000) from mod_scalars.F.
What is z_0 (the lower limit of the integration). Does it refer to the top surface or ocean bottom ?

[fabeobd]

Hi all,

I'm trying to have a closed temperature/salinity budgets using the diagnostics calculated via DIANOSTICS_TS. Assuming the temperature and salinity equations as defined by:

temp_rate = temp_hadv + temp_vadv + temp_vdiff + temp_hdiff + shflux
salt_rate = salt_hadv + salt_vadv + salt_vdiff + salt_hdiff + ssflux

Below the surface layer where the surface fluxes are not taken in account, I get a good match between the tendency (temp/salt_rate, LHS) and the ocean convergences due to advection/diffusion (RHS), see figure attached for horizontally-integrated budget over the model domain (bottom panels show the residue between LHS and RHS).

However, there is a significant mismatch in the upper layer. I divided the net surface heat flux (shflux) by Cp x rho to convert from W m-2 to degC m s-1, while for net surface salt flux are already in m s-1. Either accounting the surface fluxes (dashed blue line in the top panels) or not (magenta dashed line) don't help to close the budget (still order of magnitudes larger than the ocean convergences). Should I account for the thickness of the upper layer to convert degC m s-1 to degC s-1? The ocean_dia.nc file indicates the tendency diagnostics are in degC s-1, is that correct?

I also tried to closed the budget by integrating vertically the temp_vdiff term and comparing with the surface fluxes, again unsuccessfully (see second figure attached, annual mean). I used monthly avreaged diagnostics (for both avg/dia) in these calculations, and this is a circumpolar model of the Antarctic with 10km horizontal resolution (31 sigma levels). I understand that the time average interval might affect the zeta and the offline calculation of the layer thickness (used here to integrated the vertical diffusion component), so I'm running the model again to save these diagnostics with 5-day interval. Any other suggestion on what might be potentially wrong?

Thanks,
Fabio

[wilkin]

temp_rate = temp_hadv + temp_vadv + temp_vdiff + temp_hdiff + shflux

You shouldn't have shflux here because it was applied as the surface boundary condition to temp_vdiff and so is already included in the vertical divergence d/dz*(K_v*dT/dz) (i.e. shflux/(Cp*rho0) = K_v(dT/dz at z = zeta; notice the division is by parameter rho0 not the model in situ density).

If you vertically integrate temp_vdiff = d/dz*(K_v*dT/dz) between the limits z = -h and z = zeta then you should simply get K_v*dT/dz |z = zeta minus K_v*dT/dz |z = -h which is shflux/(Cp*rho0) minus 0.

You are correct that the diagnostic terms are in degC s-1. The units of shflux/(Cp*rho0) are degC m s-1 because of the vertical integral.

I don't know why your vertically integrated temp_vdiff in degC m s-1 does not balance shflux/(Cp*rho0) but you are correct that it should. Make sure you used the Cp and rho0 that ROMS assumes.

[fabeobd]

Hi John,

thanks for your response.

I have used rho0 (not in-situ rho) and heat capacity obtained from mod_scalars.F. I figured out an error in my dz calculation, now obtained by:

Code: Select all

Zo_w = (ds.hc * ds.s_w + ds.Cs_w * ds.h) / (ds.hc + ds.h)
z_w = ds.zeta + (ds.zeta + ds.h) * Zo_w

for tt in np.arange(0,12): # loop in time-dimension
    dz[tt,:,:,:] = numpy.diff(z_w,axis=2) # diff at vertical-dimension

By fixing the layer thickness calculation, the difference between the vertical integral of vdiff and shflux got closer (especially out of the continental shelf), but still there is a mismatch near the Antarctic coast (see fig. attached, top row).

This model includes the ice shelf module, so I'm wondering if this has anything to do with ice shelf - ocean fluxes not properly diagnosed. I will look further into this in the ice shelf code.

However, even if the vertical integral of vdiff and shflux get a perfect match, I'm puzzled by the plot of the horizontally-integrated tendencies. If the surface fluxes are included in vdiff, why the temp_rate and the sum of all other terms (temp_hadv + temp_vadv + temp_diff + temp_vdiff) do not match? If the surface fluxes at the upper layer are much larger than other terms (hdiff and advection), that shouldn't reflect in the temp_rate? I'm including a new version of the previous figure here without include the surface fluxes in the budgeting (so considering temp_rate = temp_hadv + temp_vadv + temp_vdiff + temp_hdiff).

Thanks,
Fabio

[wilkin]

Your plots seem to indicate the mismatch in your budgets is only at the Antarctic coast. So that does rather implicate the ice shelf as a possible issue - perhaps some term not added to the diagnostics, or that needs to be separately explicitly taken account of in your analysis.

I don't follow what you say about horizontally-integrated term balances. Only if your northern boundary is closed would the horizontal integral be zero. Is that what you mean?

[fabeobd]

Hi John,

thanks for your insights. I hope to find out what's missing in terms of the surface fluxes soon - I will keep it posted here.

Regarding the horizontally-integrations terms, my northern boundary is not closed so I wouldn't expect this integral to be zero. What I was expecting is that temp_rate (LHS) and advection+diffusion (RHS) balance each other at every layer. That occurs for every layer except the upper one.

So I did plot the vertically-integrated terms (attached), in an attempt to reproduce the analyses you did in your 2006 JPO paper (fig. 10), but annually averaged. As we know, the vertical diffusion term when depth-integrated should be equal to the net surface flux (panel C below). That flux entering or leaving the model domain should be transported by other terms via vertical/lateral fluxes (advection and diffusion), and result in a temperature/salinity tendency (temp_rate/salt_rate, panel A). The panel B shows the vertically-integrated LHS minus the sum of the RHS terms (C, D, E, F), which has a residue a few order of magnitude smaller than the individual terms (O(10^-5) degC/day). I think that makes sense and indicates the model has a closed budget, would you agree? So the upper layer is not direct balanced in the horizontal integral because the surface fluxes enter the ocean in that layer but are distributed in the vertical (and lateral) by individual components of the temperature/salinity fluxes - but they do balance when all terms are depth-integrated.

Cheers,
Fabio