ExoRT Port Issues

[jtruesdal]

I'm moving the email discussion for ExoRT debug to this thread.

Current Issues when running ExoRT in CESM3-planets

1. In the initial runs of the model we are seeing large heating rate values in the lowest levels of the model. When the temperature drops to unrealistic values the model dies after a few timesteps.
2. The radiation diagnostics are being written to the history tape but only appear to have values for the initial time step.
3. The current vertical structure for the atmosphere runs from 6mb at the surface up to very small pressures near the top 1.e-10 mb. This is far outside the limits of ExoRT and I am setting the variable camtop to level 31 ( .1 mb). It looks as though the radiation code loops over camtop to pverp for the flux calculations. I am thinking we should create a new initial condition with a vertical coordinate that has the same model levels but with a top of .1 mb which is within the range of ExoRT. We know that ExoRT has been tested with camtop = 1 and we can move back to that with a new IC file.
4. The radiation.F90 interface has place holders for writing out lw and sw diagnostics as well as spectral data. These routines need to be enabled and the diagnostics need to be debugged to see why we are not getting history values past the initial timestep.

[dan800]

@jtruesdal I think a more cautious approach might be beneficial here - simply call exoRT to calculate diagnostic SW and LW heating rates and continue to call the H-S Mars temperature relaxation routine. We will need to do that anyway to avoid the model wandering off in the region of the atmosphere where heating rates are not produced by exoRT (at least until we add more physics). Good to see in the code that we can now save out QRS and QRL. Inspecting those after one timestep (and comparing to the LMD model) will be useful.

[jtruesdal]

Great idea. I will look into adding that to our simple physics version. The flexibility will allow us to work in parallel on debugging ExoRT and adding gw/vdif/land_surface.

[FelixSainsbury]

So I have run the model, disabling the updates the to temperature profile by exoRT and have found that there is a significant difference between QRS and QRL in terms of their scale: For example, near the surface, QRS is on the order of 1e-9 (i know if should be low) and QRL is on the order of unity. Even moving to lower pressures, we are looking at 1e-8 vs 1e-3 - so it looks like something is going wrong in the scaling that is being passed to CESM......

[dan800]

That's not what I would expect considering the rapid heating of the atmosphere that happened when T tendencies were applied (i.e., you would expect QRS >> QRL). Is one term missing a rho*Cp factor?

[FelixSainsbury]

I'm having a look now, but something is going on with QRL as it is both heating and cooling the atmosphere -
- So we have a negative QRL at the equator, where the temperature goes highly negative, and a positive QRL at the poles where the temperature goes highly positive

[FelixSainsbury]

Currently I am just double checking that setting CAMTOP to 31 is not breaking something, but I will look at the components of QRL next to try and see where this is coming from

[storyofthewolf]

So I believe the positive QRL (long wave heating) is occurring where near-surface inversions occur. Essentially there is net positive longwave radiation accumulating in a grid box, because the inversion layer is radiating more radiation downward, than the lower/cooler layer is radiating upward

[FelixSainsbury]

It looks like running with camtop=31 and camtop =1 leads to the same results near the surface, and that the difference between QRS and QRL is apparent at low pressures as well - if we look at the top of the atmosphere, the peak QRS is about 48, whereas the peak QRL is 60000/-480000!

[dan800]

For comparison, here are the SW heating rates (K/s) from PCM-Mars calculated by Boding.

[storyofthewolf]

do you have thermal cooling rates?

[dan800]

All, I had a dig into where LTE starts to break down. According to a tech doc. by Miguel Lopez-Valverde in 2001 the merge region is between 82 and 97 km, which is about 0.08 to 0.01 Pa (note not hPa).

[jtruesdal]

Interesting the z height calculated in exort near the top (level 1) is 71km. I think a lower top for our might help stability for ExoRT.

[storyofthewolf]

Hi All, this is a random question about the .h1i. files. Is there is odd trick to these netcdf files that I am unaware of? They won't open in IDL nor a GUI viewer I use (panoply), but they do seem to open in python, and work with netcdf operators (NCO). It's a not huge deal, but I might have to rewrite my IDL script to python.

[jtruesdal]

No oddities that I'm aware of. They should be standard netcdf files. If you find something though I'm happy to pass it along to get it fixed in CAM.

[storyofthewolf]

IDL spits back NCDF_OPEN: Unable to open the file (NC_ERROR=-51), same with both IDL version 8.5, and IDL version 8.8
Panoply: ( a NASA GISS GUI Viewer), "The file does doe seem to be a valid CDM"
python seems to work fine.

It's alright, I have been trying to move my work over to python anyway, since more people/students use it, and it altogether seems a more valuable skill moving forward than legacy IDL. So I will convert my profile read script, it will just take me a little longer, that's all.

[storyofthewolf]

On a separate note, now that I am getting back into ExoRT.... the public version of ExoRT now has O2 and O3 in it, which my colleague Russell Dietrick added and has a paper in submission where he tested it in 1D radiative calculations, and some simple 3D CESM1.2.1/ExoCAM simulations, with good results compared to RRTMG. However, O2/O3 are now pulled in from the physics buffer via the prescribed_ozone_file. In Russell's implementation with ExoCAM, he has a zero'd out version of the ozone file. If for some reason, O3 climatological values are lingering in the Mars model, this could cause issues. Just something to be aware of.

(I haven;t even had a chance to use Russell's version yet...)

[storyofthewolf]

Felix, could you make the output file again, but this time with 'P0'? The equation to convert the hybrid sigma coordinates to pressures requires both PS and P0. Also, could you change the output so that there is a single timestep per file (mfilt=1). This will just make it a little faster work the file. Thanks!

[storyofthewolf]

actually, just to make it more streamlined can you add "TS" and "Q" as well? I no that water is supposed to be zero. The more information I can get directly from the netcdf, and don't have to hard code or kludge, the better. I was hoping to have the calculation be quick and easy today, but IDL keeps throwing back an error at me, so python it is...

[FelixSainsbury]

I'm not sure why P0 is not being included, it always has in my outputs before.... maybe @dan800 knows why it might not be being output. As for the other variables, I have uploaded a new file to onedrive: https://leeds365-my.sharepoint.com/:u:/g/personal/phyfs_leeds_ac_uk/EfnkWRm821BPu4mIg06ninMBYOkSQxP5RXdgYhlSXCCBSg?e=Kfatmd

[storyofthewolf]

Hi Felix, I don't think this is the write file, this link appears to take me to source.sh.
that is puzzling, I am looking back through my old files and namelists, and it appears P0 is usually, automatically included, even when empty_htapes=.true. and P0 isn't included in fincl.

[FelixSainsbury]

Ah sorry, onedrive has been playing up. This should worK https://leeds365-my.sharepoint.com/:u:/g/personal/phyfs_leeds_ac_uk/EecI5UGHW0RGvnTedHx_W2IBS2c_hOGbc7Lfn1Ccpdtm7g?e=OhEC44

And yes I have the same issue, P0 seems to automatically have been included in my previous runs!

[dan800]

I ran into an issue with making a file in Python to be read in by CESM. I had to convert it to cdf5:

nccopy -k cdf5 mars_cam_vcoords_L49_c221221_netcdf3.nc mars_cam_vcoords_L49_c221221_netcdf3_cdf5.nc

Perhaps some conversion to / from 3 or 5 might make the h1i file work for IDL?

[storyofthewolf]

Oh interesting thanks! Missed this. I will try this tomorrow.

[FelixSainsbury]

Are you making netcdf3 files with python? I make netCDF 4 files and they are read fine by CESM.

[dan800]

I think the question might be whether Eric's version of IDL can read in newer versions of netcdf. @FelixSainsbury, if you are extracting fields from a CESM run, you may want to try making a second file with
dsout.to_netcdf(outfile, format="NETCDF3_64BIT")
Just in case there's an issue with reading more recent versions.

[FelixSainsbury]

Lets see if a simpler solution works, I have converted the output to netCDF3 - does this file work in IDL? https://leeds365-my.sharepoint.com/:u:/g/personal/phyfs_leeds_ac_uk/EUdSnx8Xa1FHicBcAkAWXnEBtLG6TIivfCeiLGW01g1UcQ?e=FrKr4m

[storyofthewolf]

This file works! Thanks for your patience guys! Seems like between my laptop, my cluster (NASA's discover computer), my ID on both laptop and cluster, my python, and ExoRT, I need to sort of the netcdf. I was using netcdf4... Anyways, I'll work on sorting that out, but the new file works with my IDL. I wrote a python script that can read in the files, but the netcdf it creates now doesn't work in ExoRT (do to a netcdf issue). I'll clean all that up, but regardless I move forward with some science. I will post results in 1D soon! For the moment I am setting P0 = PS. The flux and heating rate fields are not populated in this file though.

I did reconfirm that ExoRT offline gives exactly duplicate results for clearsky fluxes and heating rates compared 3D simulations run with ExoCAM.

[storyofthewolf]

Hi All, I think the issue is the surface temperature assumption of 190 K, is way inappropriate for the overlaying atmospheres. Here is the temperature profile (from a random column selected). This contributes to the long wave cooling rates getting crazy. Note, that positive long wave delta T/time is still found underneath the inversion layers, as I commented above on Felix's plot.

[jtruesdal]

Happy to help. Since its critical to get the new code changes I would recommend a fresh checkout of the entire code not just updating components. You should run the git-fleximod update command twice after download and you shouldn't see any warnings or errors on the second run of update. I could also send you a link to a tar ball separately if needed.

The initial file should have co2 and n2 on it and be about 2.1M in size:
-rw-r--r-- 1 jet cgdcsm 1223816 Oct 29 18:09 /fs/cgd/csm/inputdata/atm/cam/inic/se/CAM_Mars_ne3np4.HS.no-topo_241016.nc

Here is my user_nl_cam file. You'll have a different path for ncdata. Note camtop is now 1 since that has been the way exort has been mostly run in the past and its unclear whether setting it to something else works properly. The relaxation parameterization is replacing static energy tendencies from levels 1 to 32 so from cam's perspective it should be receiving reasonable data from the radiation_tend call at the upper levels. exo_rad_step is 1 so radiation called every time step. I'll also email you my atm.log file which you can compare to your log. Sometimes you can spot obvious differences in the setup that way using a visual difference application.

se_statefreq = 1
nhtfrq = 1,1
interpolate_output = .true.,.true.,.true.
interpolate_nlat = 18,18,18
interpolate_nlon = 36,36,36
ncdata = '/fs/cgd/csm/inputdata/atm/cam/inic/se/CAM_Mars_ne3np4.HS.no-topo_241016.nc'
se_nsplit = 8
se_rsplit = 4
se_hypervis_subcycle = 8
se_hypervis_subcycle_q = 1
mfilt = 144
write_nstep0 = .true.
avgflag_pertape(1) = 'I'
gravit = 3.72
sday = 88642.0
mwdry = 43.34
cpair = 735.0
rearth = 3.38992e6
co2vmr = 0.97
ch4vmr= 0.
n2ovmr= 0.
f11vmr= 0.
f12vmr= 0.
n2vmr= 0.03
rad_climate = 'A:CO2:CO2','A:Q:H2O','A:N2:N2'
water_species_in_air = 'Q'
do_exo_gw = .false.
do_exo_rt = .true.
do_exo_rt_clearsky = .true.
do_exo_rt_optimize_bands = .true.
do_exo_rt_spectral = .true.
do_exo_synchronous = .false.
exo_rad_camtop = 1
exo_rad_step = 1
mars_relaxation = .true.
mars_relax_bot_p = 20
mars_relax_fincl = ''
mars_relax_linear = .false.
mars_relax_tau_sec = 1800.
mars_relax_tau_bot_sec = 864000.
mars_relax_tau_top_sec = 172800.
mars_relax_top_p = 1e-20
mars_use_hs_uv = .true.

[FelixSainsbury]

So despite there being no differences in our setup (previously there were small differences), my model still diverges from yours - the only thing I can think of that might be happening is something is different with the compiler, leading to tiny differences that then compound.

[FelixSainsbury]

Not that I currently think spending more time on this is useful - we need to work out why it still explodes!

[jtruesdal]

As you indicated the models will diverge numerically because of roundoff differences due to the different machines, compilers, libraries, etc. So the actual global values reported for each time step will be different. I provided the log to help see if there were differences in the setup which it seems you were able to correct. Are you now able to run 8 or 9 time steps now before the model blows up?

[FelixSainsbury]

Mope, which is what I meant when I said it diverges. My closest run to your setup makes it around 2 timesteps before the temperature goes negative, with the error actually blaming the HS package:

ERROR: shr_assert_in_domain: state%t has invalid value -164.444361054217
at location: 4 47
Expected value to be greater than 0.000000000000000E+000
ERROR: Invalid value produced in physics_state by package held_suarez.

I suspect that this is because we are getting rather different forcings:

FORCING DIAGNOSTICS

                             MIN                    MAX

• FT -0.730512079816920E-01 0.185332650731878E+00

• FT -0.285235543128373E-02 0.624193444509490E-02
  FM -0.182327819064238E-03 0.166082123510446E-03
  FQ 0.000000000000000E+00 0.765641897165216E-29
  FCO2 0.000000000000000E+00 0.000000000000000E+00
  @@ -2621,75 +2605,75 @@

My model is the "-" in the diff. Have you made any other code changes that wernt reflected in the version on git - e.g. files in your SourceMods directory? Otherwise my setup is the same as yours, including a clean download of the source from git - the only modifications I made at the time of this run where to make CESM run on arc4, our local cluster which uses SGE. Anyway I have sent you my log in case you can spot anything I have missed!

However I agree with Dan, at the moment it is more important to understand what is going on with exoRT, and why the near surface flux is leading to model exploding, either in two timesteps or nine! This is especially true given that arc4 is due to be decommissioned in the next two months, and uses rather old compilers!

[dan800]

Here is a chunk of python that recreates the mean zonal mean surface temperature used in the H-S case. If you want to hard-wire a surface T, you can try this:

import numpy as np

# create latitude array and tref array
nlat = 36
lat = np.arange(nlat)*5-87.5
tsurface = np.zeros(nlat)

amp = np.zeros(5)
amp[0] = 1.837e+02 
amp[1] = -3.138e+01 
amp[2] = -9.052e+00 
amp[3] = -9.720e-01
amp[4] = 2.390e+00

y = 2.*(np.radians(lat)+0.5*np.pi)

tsurface= amp[0] \
    + amp[1]*np.cos(y) \
    + amp[2]*np.cos(2*y) \
    + amp[3]*np.cos(3*y) \
    + amp[4]*np.cos(4*y)

print(tsurface)

[FelixSainsbury]

I did have a look over radiation.f90 and also did not see an error, but saying that I was somewhat relying on the comments being right (vis-a-vis units such as J/kg/s etc). When i get a chance I will take another look over, but so far I have not seen any evidence that it is in K/day.....

[FelixSainsbury]

@jtruesdal How did you set the surface temperature, I tried state%t(i,pver) and the lower atmosphere still ended up too warm. Am I using the wrong temperature in state?

[dan800]

Is the surface temperature being set by the land model? I thought @jtruesdal had said the model was being configured with a data land model.

In general, in CESM, it's not advised to directly modify/overwrite state but to apply a tendency to arrive at the state you want. Just the way it's done.

[FelixSainsbury]

@dan800 I mean in exoRT - when value is @jtruesdal feeding aerad_driver?

[jtruesdal]

I update the cam_in value for TS using the lowest level atmospheric temperature. This cam_in%ts is fed into ExoRT. The update is done in the mars_surface_init routine inside mars_cam.F90 and is called before each radiation call. Because the atmospheric state temperature is updated after each radiation call, TS should track the atmospheric changes. In the same routine I hard code values of .5 for the direct/diffuse lw/sw albedos which should probably be more like .3 for desert conditions. These will also change the surface forcing and affect the air temps.

As for Dan's question, we are running a simple physics setup for Mars which doesn't include any boundary forcing, just the atmosphere. I've spoken with our land people about how to modify the CESM land model CTSM to provide a desert planet and am currently working on adding that while you all are looking at the ExoRT errors. I'm also asking checking into providing a data land component for CAM-Mars.

[FelixSainsbury]

Just to keep everyone in the loop, I have had another pass at radiation.F90 and the functions it calls over the last few (work) days - so far i cannot see anything obviously wrong but the next step would be to go into the maths - @storyofthewolf is there a single paper that gives all the equations being solved?

[dan800]

I think it's better to verify that what is being produced from Eric's single column model is the same as the inline version.

[FelixSainsbury]

That would also make sense - I was just thinking of what I can do on my end.

[jtruesdal]

@storyofthewolf is there something I can do to help you validate the results against the column model?

[storyofthewolf]

Hey guys, my apologies. I had a conference last the previous week, and got 50+ messages behind here, and haven’t been able to catch up.

…

On Mon, Nov 25, 2024 at 8:42 AM jtruesdal ***@***.***> wrote: @storyofthewolf <https://github.com/storyofthewolf> is there something I can do to help you validate the results against the column model? — Reply to this email directly, view it on GitHub <#2 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEYFAHYLUA2O5TWJBVZRTF32CNAONAVCNFSM6AAAAABQ4PJGGCVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCMZXGQZDCMY> . You are receiving this because you were mentioned.Message ID: ***@***.***>