Add support for using relative humidity rather than specific humidity.

[russfiedler]

As mentioned in this thread COSIMA/libaccessom2#30
Initial development is here https://github.com/russfiedler/cice5/tree/rel_humid

[rmholmes]

Thanks @russfiedler. I've successfully compiled your code and am setting up two simulations to test it. However, I ran into a problem - to do a proper test I would need to reproduce the calculation of relative humidity done by sat_vapor_press_mod in order to create the forcing file to start with. This seems to be rather complicated (involving a look-up table) so I'm not confident in reproducing it in python or the like (or indeed, trying to call those fortran routines from elsewhere). Any suggestions?

Alternatively, I could do the calculation approximately with a simpler formula (like the one you had in there on your first attempt) - which will probably give me similar results.

[russfiedler]

@rmholmes Yeah, it probably is a bit of overkill for our purposes but I wanted to be mostly consistent. The actual functions that are interpolated are messy but not impossible. Just all those log10 terms rather natural logs. I think the main reason to use a lookup table was for speed in the dim dark distant days when these calulations were very slow. The full version of those lookup tables are in the MOM code in https://github.com/mom-ocean/MOM5/blob/master/src/shared/sat_vapor_pres/sat_vapor_pres_k.F90 Just ignore the lookup table and use the formulae there.

Actually, it looks like I could have used the routine compute_qs directly to calculate the specific humidity if the full module was there. So I'm not being totally consistent with what's in MOM.

Feel free to use the simpler formula that I pencilled in the. It should be close enough for your purposes.

[rmholmes]

@russfiedler the test is failing somewhere soon after the read-in of the relative humidity data. I haven't been able to glean much information from the logs. Only that atmosphere/log/matmxx.pe00000.log shows that the error is somewhere in the rhuss read-in stage (i.e. there is a field_update-data-file but no checksum-matmxx for this field). I've had a look through the code and nothing stands out - but I don't have any experience with these bits.

My test run is at /home/561/rmh561/access-om2/1deg_jra55_ryf_RHtest, with corresponding control at /home/561/rmh561/access-om2/1deg_jra55_ryf (same executable, just without the forcing.json and input_ice.nml changes). My code and inputs are in /g/data/e14/rmh561/access-om2/. Any hints/suggestions are appreciated. If left on my own my next step would be to compile with debug options.

[russfiedler]

@rmholmes From the looks of things you haven't changed the namcouple file to allow passing relative humidity to cice

Field 08 : 2m air humidity
##########
qair_ai qair_i -1 -1 1 NA EXPORTED <-------------------

This should be

relh_i relh_i -1 -1 1 NA EXPORTED

though I think it would be more consistent to change the name in forcing.json to relh_ai and the line becomes

relh_ai relh_i -1 -1 1 NA EXPORTED

There should have been an error message message somewhere but (like you) I can't find it. There was no backtrace for yatm either which was a pain.

[rmholmes]

@russfiedler I'm still having some problems with this. I've got it running but there is still a ~10Wm-2 difference in the latent heat flux in the two simulations. This has a significant impact on the global OHC, so I don't think it's acceptable (this plot is from the first month):

I played around with the limits on the temperature input into escomp (this was causing some crashes that I fixed) and checked the units there. I think the problem is either coming from a problem with another input (e.g. I haven't checked the SLP units), or is actually arrising because the escomp calculation of e_sat is different to the simple formula I've used to create the forcing file.

I am off on two weeks leave tomorrow and will revisit this when I get back. In the mean-time, if you're keen to keep working on this my test runs are at:
/home/561/rmh561/access-om2/1deg_jra55_ryf, https://github.com/rmholmes/1deg_jra55_ryf/tree/RHtest_cont
/home/561/rmh561/access-om2/1deg_jra55_ryf_RHtest, https://github.com/rmholmes/1deg_jra55_ryf/tree/RHtest

[russfiedler]

@rmholmes SLP units are Pa and that's what it says in the files. Your changes to the temp limits look good. That's a pretty large bias that's coming out of the simulation. I wonder if there's a problem with the nonlinearity of SH <-> RH conversion and averaging though I guess that's been looked at. The simplified equation for q should be ok. I'll check the JRA-55 documentation to see if I can find out what they used.

[rmholmes]

@russfiedler did you make any progress on this?

[russfiedler]

@rmholmes No, I haven't had a chance to look at it. Having said that I guess we should look at using the same formulae for RH<->SH as here https://climate.mri-jma.go.jp/~htsujino/docs/JRA55-do/v1_5-manual/User_manual_jra55_do_v1_5.pdf Section 4B page 7. Should be easy to code up.

Hmm... maybe we should have these options in the code anyway since I think RH is calculated somewhere along the way.. I'll have to think about that. There's also some comments there about the latent heat of fusion and other issues. Basically use Gill (1982) rather than Large and Yeager (2004).

[russfiedler]

"3. Recommended methods for driving ocean–sea-ice models using the JRA55-do dataset
It is basically recommended that the users follow the experimental protocol of OMIP (Griffies et al., 2016) for driving their
ocean–sea-ice models using the JRA55-do dataset. Exception is the computation of properties of moist air, such as density, specific
heat, latent heat of vaporization and sublimation of water vapor, and saturation specific humidity. We advocate using a set of formulas
given by Gill (1982) for computing them, which is thought to be more accurate than the simple and cost-effective set of formulas used
for the current CORE/OMIP framework as given by Large and Yeager (2004). Appendix B compares the two formulas.
After we accumulate experiences with the JRA55-do dataset for driving ocean–sea-ice models, we will provide recommendation
on how ocean–sea-ice models should be suitably forced using this datas"

[rmholmes]

@russfiedler sorry for the radio silence.

It seems to me that there are two separate issues here:

1. What formula's should be used in the model internally for converting the JRA-55do dataset's provided specific humidity + other parameters into the quantities needed to force the model. As you've mentioned, what is currently being used in ACCESS-OM2 may not be what is recommended. This issue affects all ACCESS-OM2 simulations run up to now - and so I feel like it needs it's own github issue.

2. How we tweak the code so that we can provide RH as an input but end up with the same results as if we'd instead provided QH. For this problem I think the formula's are not that important - we just have to make sure that the same formula is used both in the model to convert the RH input to QH as I use in the pre-processing to create the RH input files from the QH files.

So, one approach to solve the present problem is to just use the simpler formula for RH->QH that I'm already using in the pre-processing (basically resurrecting your earlier solution). If that sounds ok to you, I'll have a go at it when I get some time.

[rmholmes]

Success! If I use the simple formula then I get a good match. There are still some small differences (note the color-scale limits) which are probably not unexpected given the model runs aren't exactly reproducible?:

I've put a pull request here: #59 so that we can discuss the changes. Please have a look and let me know what you think.

[russfiedler]

Good stuff. You wouldn't expect them to be the same because the RH<->SH conversion that you use is different to that in Gill 82 which is what is used for JRA-55 (bit hard to tell with different base logs and T(K) vs T(C). However, I think it's pretty clear that it's a much closer match than the GFDL form.

[rmholmes]

I still think we should expect them to be the same - up to the model runs not being reproducible because of machine precision/noise etc. All I am doing is taking the JRA55-do provided QH files, converting them to RH using this simple formula in pre-processing (producing an RH file), and then the online model code just converts this input RH straight back to QH using the same simple formula. It should be immaterial what happens after this (i.e. what formulas are subsequently used).

Unless there is some time interpolation operations that give you differences because the formula is non-linear? I don't think this is the case...

In any case, I think these differences are small enough not to worry about.

[russfiedler]

Aha, I thought you were taking JRA-55 supplied RH files rather than cooking them yourself. If you're supplying the new fields in single precision that's way more than enough to ensure chaos rules as expected! Even in double precision I wouldn't expect the conversions to be exactly reversible over the full range of values.