Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

ACCESS-OM experiment (MOM5 with -DACCESS) does not conserve heat well #175

Closed
nichannah opened this issue Jun 20, 2017 · 16 comments
Closed

ACCESS-OM experiment (MOM5 with -DACCESS) does not conserve heat well #175

nichannah opened this issue Jun 20, 2017 · 16 comments

Comments

@nichannah
Copy link
Contributor

nichannah commented Jun 20, 2017 •
edited
Loading

From Fabio's email:

In the access.out file we can search for that segment:

Error in heat content change = 1.09666002438988800E+15 J.
Error in rate of heat content change = 3.66547652005596255E-05 W/m^2.

which is “Measures of global integrated tracer conservation over multiple time steps”. It should be printed in the .o files too, but I don’t have permittions to access these files in your experiment.

Following what Stephen Griffies said to me, the error in heat content change should be O(+11) Joules or less.

To ilustrate what that means, I plot the figure below with my ACCESS-OM (JRA-55 forcing) where the net heat flux crossing the surface is not equal of the time derivative of the ocean heat content:

access_om_heat_budget

and when I calculate the same for the Richard Matear’s model (MOM4-SIS, also forced with JRA55), both match perfectly:

mom4sis_heat_budget
@nichannah
Copy link
Contributor Author

See attached output containing budget diagnostics for a short 1 deg run.

access-om2.out.txt

nichannah pushed a commit to nichannah/MOM5 that referenced this issue Jun 20, 2017
Set melt consistently in SBC. mom-ocean#175.
1d9b24d
@StephenGriffies
Copy link
Contributor

StephenGriffies commented Jun 20, 2017 via email

@nichannah
Copy link
Contributor Author

nichannah commented Jun 21, 2017 •
edited
Loading

Pull request #176 makes improvements to SBCs are handled for ACCESS. The fixes are:

  1. salt flux and reference salinity are not used to calculate melt, instead melt comes from the coupler. This has several benefits:

    • the reference salinities between ice and ocean don't need to be the same (there have been problem with this in the past)
    • melt is defined in only one way, instead of two. Previously there could be a discrepancy between melt being passed into MOM via the melt variable and that included in PME.

  2. An indexing error was fixed. Previously melt was added to PME at one grid point and subtracted from PME_river at another.

This does not fix the large heat content errors.

nichannah added a commit that referenced this issue Jun 21, 2017
@nichannah
Merge pull request #176 from nicjhan/175-heat-conservation
b84248b 
Set melt consistently in SBC. #175.
@russfiedler
Copy link
Collaborator

Here's what I reported to Nic, Fabio and Simon yesterday:

I've had a look at auscom_ice_formation_new() and have come to the conclusion that the method is slightly faulty.

What happens is that starting at the deepest layer, kmxice, a test is made to see if ice can be formed. If so, the temperature is raised due to latent heat release. If layers above are warmer than the freezing point they can melt this ice and become cooler. They can also contribute more frazil if cool enough. Also, any accumulated ice during the coupling period can be melted again. This is all fine.

The problem is that the heat due to frazil is calculated via

PTR_FRAZIL = rho_cp * frazil_factor * ( PTR_TEMP - TEMP_BEFORE) * PTR_THICK
= rho_cp * frazil_factor * ( PTR_TEMP - TEMP_BEFORE) * dzt

This is wrong and dates way back to the original implementation of the algorithm. I looked at the POP code and the error has been copied from there. The correct amount of heat is

cp_ocean * frazil_factor * ( PTR_TEMP - TEMP_BEFORE) * Thickness%rho_dzt(iisc:iiec,jjsc:jjec,k)

That is, the real density should be used rather than a reference density.
See examples in mom/src/mom5/ocean_tracers/ocean_frazil.F90.

I also think that the changes of temperature should be made with respect to mass rather than just thickness/volume. This should make the diagnostic consistent.

This means changing the thickness pointer to rho_dzt. This means that POTICE, QICE and AQICE now have units degC*kg/m^2

 PTR_THICK => Thickness%rho_dzt(iisc:iiec,jjsc:jjec,k)

and

 PTR_FRAZIL = cp_ocean * frazil_factor * ( PTR_TEMP - TEMP_BEFORE) * PTR_THICK

and

 Ocean_sfc%frazil(i,j) = QICE(i,j) * frazil_factor * cp_ocean / float(ATIME)

Does this look right?

Also, for some reason, AQICE is being calculated in src/mom5/ocean_tracers/ocean_frazil.F90 but it doesn't seem to be used anywhere but in auscom_ice_formation_new. Nic, can you see it anywhere? It also has inconsistent units with what is in auscom_ice.F90. Maybe legacy code which should be removed?

@russfiedler
Copy link
Collaborator

russfiedler commented Jun 22, 2017 via email

@StephenGriffies
Copy link
Contributor

Russ, for the zstar model, rho_dzt = rho0*dzt with rho0=1035. So I think the two formula you wrote will give the same results. Thoughts?

Another possibility is that the frazil contribution to the heat budget may not be computed properly for cases where frazil kicks in for k>1. Ocean_sfc%frazil(i,j) is a 2d array. When frazil is happening for k>1 then there needs to be a vertical sum to accumulate the interior frazil heating impacts. Is that done properly?

It may be useful to run a test where frazil is allowed only in the k=1 cell.

Steve

@russfiedler
Copy link
Collaborator

russfiedler commented Jun 22, 2017 via email

@fabiobdias
Copy link

Hi Russ and Steve,

I've done the test with kmxice=1 but the error magnitude diagnosed online still the same:

Error in heat content change = 2.08752682505739200E+15 J.
Error in rate of heat content change = 6.82889551358594902E-05 W/m^2

@StephenGriffies
Copy link
Contributor

OK. Thanks for the test. As Russ hypothesizes, the problem must be in the melting part of the scheme...

@russfiedler
Copy link
Collaborator

russfiedler commented Jun 29, 2017 via email

@russfiedler
Copy link
Collaborator

With this change cp_over_lhfusion needs to be calculated in the code rather than being initialised in a declaration statement as cp_ocean is no longer a parameter.

@StephenGriffies
Copy link
Contributor

Wow, that is a gotcha! Good catch. I hope this is the smoking gun.

@russfiedler
Copy link
Collaborator

I just checked by scaling the frazil in one of Fabio's runs by cp_mom/cp_fms and it it matches the error to about 6 decimal places. Fabio is doing a quick run but I'm sure we've got it.

@fabiobdias
Copy link

Just to confirm that Russ really got it (yew!!!):

Error in heat content change = 1.63592161600000000E+09 J.
Error in rate of heat content change = 5.35156609725172690E-11 W/m^2.

I'm sending the new auscom_ice.F90 to Russ so he should update here soon.

@StephenGriffies
Copy link
Contributor

Great news! And wonderful that the code fix is rather simple to implement. This is really wonderful. Great job Russ!!

@nichannah
Copy link
Contributor Author

Well done @russfiedler! @fabiobdias I have not merged this fix. Can you send me the changes or do a PR? Thanks!

@spirillen spirillen mentioned this issue Feb 8, 2025
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@nichannah @StephenGriffies @fabiobdias @russfiedler