Support newer JRA55-do

[aekiss]

ACCESS-OM3 should support the latest JRA55-do v1.5.0 but our test configs currently use /g/data/ik11/inputs/cime/inputdata/2023-03-10/ocn/jra55/v1.3_noleap, which is apparently JRA-55, not JRA55-do, and an old version of JRA-55 to boot.

$ ncdump -h /g/data/ik11/inputs/cime/inputdata/2023-03-10/ocn/jra55/v1.3_noleap/JRA.v1.3.v_10.TL319.2007.171019.nc
...
// global attributes:
		:title = "v_10 -- JRA-55 based atmosphere data on TL319 grid" ;
		:Source = "Tsujino et al., 2017. JRA-55 Based Surface Dataset for Driving Ocean-Sea-Ice Models. Ocean Modell.\n",
			"Kobayashi et al., 2015. The JRA-55 Reanalysis. J. Met. Soc. Jap., 93, 5-48 " ;
		:Conventions = "CF1.0" ;
		:history = "Original file v_10.2007.18Oct2017.nc reformatted by whokim Thu Oct 19 23:59:42 MDT 2017" ;
		:notes = "JRA-55 TL319 data (v1.3, 2017 Oct 18) reformatted for use in CESM, removed all Feb 29\\\'s, calendar now noleap" ;

DATM says it supports "JRA-55 v1.3 forcing data"
https://escomp.github.io/CDEPS/versions/master/html/datm.html#supported-data-modes

CORE_IAF_JRA (datm_datamode_jra_mod.F90)
In conjunction with JRA-55 Project, provides the atmosphere forcing when coupling an active ocean model with observed atmospheric forcing. This mode and associated data sets implement the JRA-55 v1.3 forcing data.

It's unclear whether DATM would also support a newer -do version. More investigation needed.

[aekiss]

@ezhilsabareesh8 these input files in the CIME MOM6-CICE6 config apparently contain different materials deposited from the atmosphere to the surface. We don't need them for ACCESS-OM3 (and they aren't available from anyway) so you can leave them out.

    - /g/data/ik11/inputs/cime/inputdata/2023-03-10/atm/cam/chem/trop_mozart_aero/aero/aerosoldep_WACCM.ensmean_monthly_hist_1849-2015_0.9x1.25_CMIP6_c180926.nc # datm aerosol stream
    - /g/data/ik11/inputs/cime/inputdata/2023-03-10/lnd/clm2/ndepdata/fndep_clm_hist_b.e21.BWHIST.f09_g17.CMIP6-historical-WACCM.ensmean_1849-2015_monthly_0.9x1.25_c180926.nc # datm lnd stream

(We do have a dust input /g/data/ik11/inputs/access-om2/input_bgc_20220224/mom_01deg for the BGC version of ACCESS-OM2, but that is read directly by MOM5 so we shouldn't use it in your modified DATM.)

[ezhilsabareesh8]

Thanks @aekiss. In ACCESS-OM3, the rainfall and snowfall fluxes are calculated from the precipitation fluxes in these lines in datm_datamode_jra_advance subroutine. The rainfall and snowfall fluxes are directly advertised here.

I declared the rainfall and snowfall fluxes separately as

real(r8), pointer :: strm_prrn(:) => null() ! Rainfall flux
real(r8), pointer :: strm_prsn(:) => null() ! Snowfall flux

and added the following lines to read the variables.

call shr_strdata_get_stream_pointer( sdat, 'Faxa_prrn' , strm_prrn , rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return

call shr_strdata_get_stream_pointer( sdat, 'Faxa_prsn' , strm_prsn , rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return

In access-om3, the rainfall and snowfall fluxes are calculated from the precipitation flux based on the air temperature (Sa_tbot). Refer here.

I changed these lines as

if (Sa_tbot(n) < tKFrz ) then ! assign precip to rain/snow components

      Faxa_rainl(n) = 0.0_R8
      
      Faxa_snowl(n) = strm_prsn(n).               ! Snowfall flux

else

      Faxa_rainl(n) = strm_prrn(n)                  ! Rainfall flux
      
      Faxa_snowl(n) = 0.0_R8

endif

I am keeping the condition if (Sa_tbot(n) < tKFrz ), that is, if the air temperature (Sa_tbot) is less than the freezing temperature of fresh water (tKFrz = 273.15 K) then the nodes are assigned to snowfall flux otherwise it will be assigned to Rainfall flux. Is this correct?

[aekiss]

Hi @ezhilsabareesh8 since we already have separate snow and rain data, I think we should advertise them without alteration, i.e. just set

Faxa_snowl(n) = strm_prsn(n)              ! Snowfall flux
Faxa_rainl(n) = strm_prrn(n)              ! Rainfall flux

with no if-else.

There is no need to adjust things based on air temperature - that was only needed in the original setup because they weren't already separate. Setting one or the other to zero based on air temperature will make it impossible to have a combination of rain and snow at the same time, even if that's what was in the JRA55-do data.

[aekiss]

See PR #45

[aekiss]

As a first step we should support JRA55-do v1.4 RYF from /g/data/ik11/inputs/JRA-55/RYF/v1-4. This requires code changes (PR #45) and also config changes in a JRA55do1.4.0-RYF fork of the MOM6-CICE6 configuration https://github.com/COSIMA/MOM6-CICE6, following
https://github.com/COSIMA/1deg_jra55_ryf/blob/master/atmosphere/forcing.json
(see #34)

When the RYF is working we can then we support IAF JRA55do v1.5.0 from /g/data/qv56/replicas/input4MIPs/CMIP6/OMIP/MRI/MRI-JRA55-do-1-5-0 in a JRA55do1.5.0-IAF fork of the MOM6-CICE6 configuration https://github.com/COSIMA/MOM6-CICE6 - see
https://github.com/COSIMA/1deg_jra55_iaf/blob/jra55v150/config.yaml
and
https://github.com/COSIMA/1deg_jra55_iaf/blob/jra55v150/atmosphere/forcing.json

[dougiesquire]

Is this going to be an issue:

From the CDEPS documentation

A model cannot loop over partial years, for example, from 1950-Feb-10 through 1980-Mar-15.

The ACCESS-OM2 RYF forcing is 1 May 1990 - 30 April 1991, right?

[aekiss]

Good question, but this won't be an issue with RYF. Although the data comes from 1 May 1990 - 30 April 1991, the time data in the netcdf files has been changed to nominally run from 1900-01-01 to 1900-12-31, ie it acts like there's a jump at the start of May in that year, e.g.:

$ ncdump -t -v time /g/data/ik11/inputs/JRA-55/RYF/v1-4/RYF.rlds.1990_1991.nc
netcdf RYF.rlds.1990_1991 {
dimensions:
	time = UNLIMITED ; // (2920 currently)
	lat = 320 ;
	bnds = 2 ;
	lon = 640 ;
variables:
	double time(time) ;
...
		time:units = "days since 1900-01-01" ;
		time:calendar = "noleap" ;
...
data:

 time = "1900-01-01", "1900-01-01 03", "1900-01-01 06", "1900-01-01 09",
    "1900-01-01 12", "1900-01-01 15", "1900-01-01 18", "1900-01-01 21",
...
    "1900-12-31", "1900-12-31 03", "1900-12-31 06", "1900-12-31 09",
    "1900-12-31 12", "1900-12-31 15", "1900-12-31 18", "1900-12-31 21" ;
}

This is on a noleap calendar, and both 1900 and 1991 are non-leap years.

[ezhilsabareesh8]

Issue #45 - Error in JRA55do-datamode Test Run Due to dtmax/dtmin Ratio Violation

   atm : model date        10101           0 
   atm : model date        10101        3600
  (shr_strdata_readstrm) reading file ub: /g/data/ik11/inputs/JRA-55/RYF/v1-4/RYF.rsds.1990_1991.nc       3
  (shr_strdata_advance) ERROR: for stream            4
  (shr_strdata_advance) ERROR: dday =            0
  (shr_strdata_advance) ERROR: dtime, dtmax, dtmin, dtlimit =
  0.125000000000000   0.250000000000000  0.125000000000000 1.50000000000000
  (shr_strdata_advance) ERROR: ymdLB, todLB, ymdUB, todUB =        10101 5400       10101       16200
  ERROR: (shr_strdata_advance) ERROR dt limit for stream, see atm.log output

The test run of JRA55do-datamode for a one-month period has encountered an error, which appears to be caused by the dtmax/dtmin ratio exceeding the specified dtlimit of 1.5 (dtlimit documentation). The error log indicates that there is an issue with the data stream, particularly when reading the file RYF.rsds.1990_1991.nc.

Upon examining the data, it was observed that the delta time (dt) is uniformly set to 0.125, corresponding to 3-hourly data, available from 0:00 hours to 21:00 hours each day. The problem might arise from the code reading data from 0:00 hours to the next day's 0:00 hours, omitting the data between 21:00 and 0:00. This discrepancy can lead to an incorrect calculation of dtmax as 0.25 (equivalent to 6 hours) and dtmin as 0.125, which results in a dtmax/dtmin ratio greater than the dtlimit of 1.5.

To resolve this issue, further investigation is required to understand how the code reads the data and to adjust it accordingly. The goal is to ensure that data from 21:00 to 0:00 is appropriately included in the reading process and that the dtmax/dtmin ratio remains within the specified limit.

[aekiss]

Thanks @ezhilsabareesh8. The increment in the time axis is uniformly 0.125hr, so this could not cause the ratio issue.

However this 12-month data stream is repeatedly cycled in a multi-year run and we need to check how the cycling is handled in the code.

The time axis in the .nc file starts at 00:00hr on 1 Jan and ends at 21:00hr on 31 Dec, but maybe the code expects it to end at midnight on 31 Dec (ie cover a complete year)? If so, maybe it compares the second-last (not the last) element of the time axis to the first element to work out the time increment when the data stream cycles around. With our data this would be 18:00hr on 31 Dec, resulting in an increment of 0.25hr to 00:00 on 1 Jan.

[ezhilsabareesh8]

Thanks @aekiss. The previous error encountered in JRA55do-datamode due to ratio issue has been resolved by making modifications to the datm.streams.xml file. The issue was related to the offset setting, which was adjusted to zero.

I found that the time settings in the JRA55 dataset differed from those in the JRA55-do dataset. Specifically, the JRA55 data started at 1:30 hrs and ended at 22:30 hrs, while the JRA55-do data started at 0:00 hrs and ended at 21:00 hrs.

The problem arose due to a predefined offset of 5400 seconds (1.5 hrs) in the SWDN (surface_downwelling_shortwave_flux_in_air) variable within the datm.streams.xml, leading to an increase in the time dtime to 6 hours. This increase caused the dtmax/dtmin ratio to exceed the specified dtlimit.

In the code, the cycling is managed based on the stream's time of the day (tod) lower and upper bounds (refer dshr_strdata_mod.f90), which are obtained from the input files.

[ezhilsabareesh8]

I was able to run the JRA55do-datamode for a one-month period. The grids of the input mesh file and the JRA55-do dataset are different. Specifically, the input mesh file in datm.streams.xml is in a displaced-pole grid, while the data grid of the JRA55-do dataset is in a lat-lon grid.

Surprisingly, despite the grid discrepancies, the simulation ran without any error, possibly due to both grids having the same number of total elements (204800). However, the output data shows discrepancies in the land-sea mask in the arctic, as illustrated in the attached figure, which maybe due to the grid incompatibility.

[dougiesquire]

Thanks @ezhilsabareesh8. As we discussed, I think it's just luck that this ran because the number of grid points in the JRA55 and JRA55-do grids happens to be the same. I'll generate a new mesh file for the JRA55-do grid today. Hopefully that fixes up the weirdness in the northern hemisphere.

[aekiss]

Great to see this runs, despite the incompatible grids.

Specifically, the input mesh file in datm.streams.xml is in a tripolar grid, while the data grid of the JRA55-do dataset is in a displaced pole grid.

Actually, is the mesh file derived from the CESM displaced-pole grid, and the JRA55-do is on a lat-lon grid, right?

[ezhilsabareesh8]

Thanks @aekiss. Yes, that is correct, I have updated my comment.

[dougiesquire]

@ezhilsabareesh8, I've created a mesh file for the JRA55-do grid. You can find it here: /g/data/tm70/ds0092/model/inputs/access-om3/JRA55do-ESMFmesh.nc. Let me know how you go with it.

[dougiesquire]

Sorry, this should've clicked sooner, but the issue with the plots above is because we haven't properly accounted for the MOM grid when plotting. We need to use the correct lons/lats from the static files, e.g.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

static = xr.open_dataset(
    "/scratch/tm70/ek4684/access-om3/archive/MOM6-CICE6_incompatible_mesh/output000/GMOM_JRA.mom6.static.nc"
)
ds = xr.open_dataset(
    "/scratch/tm70/ek4684/access-om3/archive/MOM6-CICE6_incompatible_mesh/output000/GMOM_JRA.mom6.sfc_0001_01.nc",
    use_cftime=True
)

SST = ds["tos"].isel(time=0).sortby(["xh","yh"]).assign_coords(
    {"geolon": static["geolon"], "geolat": static["geolat"]}
)

plt.figure(figsize=(8, 4))
ax = plt.axes(projection=ccrs.Robinson())
SST.plot.pcolormesh(
    ax=ax,
    x="geolon", y="geolat",
    transform=ccrs.PlateCarree(),
    vmin=-10, vmax=30, extend='both',
    cmap="magma", 
    cbar_kwargs = {
        'label': 'SST (°C)',
        'fraction': 0.03,
        'aspect': 15,
        'shrink': 0.7
    }
)
ax.coastlines();

produces

[aekiss]

Ah, good catch, I hadn't realised that hadn't been done

[ezhilsabareesh8]

Thanks @dougiesquire. I was able to run the JRA55do datamode with the /g/data/tm70/ds0092/model/inputs/access-om3/JRA55do-ESMFmesh.nc mesh file. Using the longitude and latitude for plotting fixed the issue with the plots.

[micaeljtoliveira]

My understanding is that this is now done in ACCESS-NRI/access-om3-configs#6. Please reopen this issue if that's not the case.