Are significant digits of P24 in GRIB the same as double precision in ECHAM6/NETCDF?

[l975421700]

Hi All

I am not sure if this is the right place to ask this question, but I meet this problem while coding in ECHAM6. My understanding & questions are described below. If you have any ideas, please let me know.

1, In ECHAM6, the calculation is done by default in double precision (dp/wp), while the output for many variables/streams is in p16 for GRIB. I understand dp has 15 and sp has 7-8 significant digits. How about p16 and p24 for GRIB files? I googled a lot but found no answers.

2, If I want to output variables in double precision, shall I just specify "bits=24" when adding the variables to the stream? Is there anything else I need to do?

3, In the GRIB file, I can see the DTYPE of my desired variable is p24 with "cdo -sinfo". But when I import this data to python with xarray and cfgrib engine, it is float32. Even when I transfer it to float 64 with cdo, it has only around 9 significant digits. Does it mean p24 has only 9 significant digits?

4, Can I output float64 in netcdf format directly?

Regards, Qinggang

[pgierz]

Moin Qinggang,

First: this is the right place :-) while the esm tools team doesn't provide direct model or Fortran support, we do host the forum. It is the best "collective place" we could think of for a question/answer style forum.

precision with grib might be a bit of a tricky topic, especially given that echam I/O is not so straightforward. It may help if you could hypothetically write down what you would expect to happen on the analysis side. What would your python code produce if Fortran did what you wanted?

Some questions from my side to be able to help better:

• why is the bit precision so important?
• how are the variables defined Fortran internally (please post how it is allocated, please with file and line numbers)
• how is it saved? (Post the snippet that writes to the netcdf file please)

I can't promise anything: there is a good chance that I will redirect you to the mpi-met forum to ask there instead, but we can try :-)

[l975421700]

Hi Paul

Thanks. Below is my answer. It is indeed a general question about the precision of GRIB, I am not sure whether it is worthwhile to go into this specific case together?

1, Bit precision is important because I am testing quantity conservation. If grib p24 has only 9 significant digits, there would always be deviations, even 9 order of magnitude smaller compared to largest values. It prevents me from conclude whether it results from my code or from this transition.

2, In tagging branch, the quantity is defined in mo_memory_wiso.f90, line 96: REAL(dp), POINTER, PUBLIC :: tagmap(:,:,:). And it is read from readtagmap.f90. It is then changed in the file mo_physc.f90 line 2349.

3, It is saved in file mo_memory_wiso.f90:
CALL add (wiso,'tagmap', tagmap, code=90, longname='map of tagging area', &
bits=24, units='none', leveltype=W_ISO)

Regards, Qinggang

[denizural]

Hi Qinggang,

Can I output float64 in netcdf format directly?

As my previous supervisor said, "GRIB is an illness and NetCDF is the medicine". I don't know why someone would use GRIB in climate modelling instead of NetCDF.

CDO, Xarray, ...

I am not sure if they tell you the exact precision of the data or the allocation of the data in the current session. You can get the correct answer by going to the source code of the subroutine that writes files in Echam. If the calculations are performed in DP but written in SP but you need high accuracy data, then I would branch-off & fork the model and would implement my DP output routine there.
Also as an exercise, it would be an interesting thing to look at the corresponding grid-base differences and overall RMSE (SP vs DP output).

[l975421700]

Thanks.

That's good to know. When I print the variables in the subroutine, it is in double precision.

I thought I can output netcdf file directly by specifying in the yaml file: "echam:add_namelist_changes:namelist.echam:runctl:out_filetype:2". But in this case, I do not get my desired streams wiso and echam at all. See experiments /work/ollie/qigao001/output/echam-6.3.05p2-wiso/pi/pi_echam6_1d_021_3.31.

Regards, Qinggang

[denizural]

Hi Qinggang,
I think I am getting close to your problem. I will write a related answer but before that I noticed one important thing:
Your log file
/work/ollie/qigao001/output/echam-6.3.05p2-wiso/pi/pi_echam6_1d_021_3.31/log/pi_echam6_1d_021_3.31_echam_compute_20000101-20000131_15764526.log
is 6GB in size. Maybe you would like to check its contents. Excessive file I/O can reduce model performance.

[l975421700]

Hi Deniz

Thanks. I am developing ECHAM6 so I intentionally wrote out many diagnostic information in that experiment. Such print statements were removed after I got the info.

Regards, Qinggang

[denizural]

Hi Qinggang,

I thought I can output netcdf file directly by specifying in the yaml file: "echam:add_namelist_changes:namelist.echam:runctl:out_filetype:2". But in this case, I do not get my desired streams wiso and echam at all.

ESM-Tools is actually setting it correctly. If you have a look at the namelist.echam, the runctl section has out_filetype = 2. The ECHAM log also verifies that:

Open output streams:

base filename /work/ollie/qigao001/output/echam-6.3.05p2-wiso/pi//pi_echam6_1d_021_3.31/run_20000101-20000131/work

file                                             stream                fileID  lpost lrest

pi_echam6_1d_021_3.31_200001.01_accw             accw                  63      T     T
pi_echam6_1d_021_3.31_200001.01_js_wiso          js_wiso               65      T     T
pi_echam6_1d_021_3.31_200001.01_la_wiso          la_wiso               67      T     F
pi_echam6_1d_021_3.31_200001.01_accw_wiso        accw_wiso             69      T     T
pi_echam6_1d_021_3.31_200001.01_echam.nc         sp                    71      T     F
pi_echam6_1d_021_3.31_200001.01_echam.nc         gl                    71      T     T
pi_echam6_1d_021_3.31_200001.01_echam.nc         g3b                   71      T     T
pi_echam6_1d_021_3.31_200001.01_echam.nc         scm                   71      T     F
pi_echam6_1d_021_3.31_200001.01_ma.nc            ma                    73      T     F
pi_echam6_1d_021_3.31_200001.01_surf.nc          surf                  75      T     T
pi_echam6_1d_021_3.31_200001.01_wiso.nc          wiso                  77      T     T
pi_echam6_1d_021_3.31_200001.01_sf_wiso.nc       sf_wiso               79      T     T
pi_echam6_1d_021_3.31_200001.01_co2.nc           co2                   81      T     T
pi_echam6_1d_021_3.31_200001.01_spim.nc          spim                  83      T     T
pi_echam6_1d_021_3.31_200001.01_glim.nc          glim                  85      T     T
pi_echam6_1d_021_3.31_200001.01_aclcim.nc        aclcim                87      T     T
pi_echam6_1d_021_3.31_200001.01_g3bim.nc         g3bim                 89      T     T
pi_echam6_1d_021_3.31_200001.01_glday.nc         glday                 91      T     T
pi_echam6_1d_021_3.31_200001.01_g3bday.nc        g3bday                93      T     T
pi_echam6_1d_021_3.31_200001.01_g3bid.nc         g3bid                 95      T     T
pi_echam6_1d_021_3.31_200001.01_jsbid.nc         jsbid                 97      T     T
pi_echam6_1d_021_3.31_200001.01_sp6h.nc          sp6h                  99      T     T
pi_echam6_1d_021_3.31_200001.01_g3b1hi.nc        g3b1hi               101      T     T

As you can see some Wiso outputs are NetCDF but some are not. So why is that? Here is the reason:

All Wiso variables from ECHAM are converted to NetCDF because you are setting it to 2 (NetCDF).

./echam/mo_memory_streams.f90:125:      CALL new_stream (surf_wiso ,'sf_wiso',ftyp,lrerun=.FALSE.,               post_suf='_sf_wiso')
./echam/mo_memory_streams.f90:142:        CALL set_stream (surf_wiso ,lrerun=.TRUE.  ,rest_suf= '_sf_wiso')

The other variables are actually coming from JSBACH: /work/ollie/qigao001/model_codes/echam-6.3.05p2-wiso/src/src_jsbach3/mo_jsbach_interface.f90

IF (lwiso) THEN   ! deniz: starting at line 1531

  CALL new_stream(theLand%IO_wiso_stream, 'js_wiso', filetype=theOptions%FileType, ztype=theOptions%FileZtype, &
                  lpost=theOptions%OutputModelState, lrerun=.TRUE.)
  CALL default_stream_setting(theLand%IO_wiso_stream, lpost=theOptions%OutputModelState, repr=LAND, lrerun=.TRUE., &
                  table=land_table)

  CALL new_stream(theLand%IO_diag_wiso_stream, 'la_wiso', filetype=theOptions%FileType, ztype=theOptions%FileZtype, &
                  lpost=.TRUE., lrerun=.FALSE.)
  CALL default_stream_setting(theLand%IO_diag_wiso_stream, lpost=.TRUE., repr=LAND, lrerun=.FALSE., table=ldiag_table)

  ! Stream that contains accumulated values weighted by land fraction, necessary to compute water and energy
  ! balance in ECHAM postprocessing
  CALL new_stream(theLand%IO_wiso_accw, 'accw_wiso', filetype=theOptions%FileType, ztype=theOptions%FileZtype, &
                  lpost=.NOT. theOptions%Standalone, lrerun=.NOT. theOptions%Standalone)
  CALL default_stream_setting(theLand%IO_wiso_accw, lpost=.NOT. theOptions%Standalone, lrerun=.NOT. theOptions%Standalone, &
                  repr=LAND, table=ldiag2_table)

So what is the value of filetype=theOptions%FileType, that seems to be responsible for the file output type.

Source: /work/ollie/qigao001/model_codes/echam-6.3.05p2-wiso/src/src_jsbach3/mo_jsbach.f90
227 file_type = GRIB

This is the default initialization and it is updated by the namelist read, that comes just after that.

I am sure that if you also set JSBACH output type to 2 via ESM-Tools, as you did for ECHAM, you will get NetCDF for other WISO outputs.

The log file tells you the output selected for JSBACH files. Eg. at line 877

0:  jsbach_config: JSBACH output file type:            1
0:  jsbach_config: JSBACH output file compression:            1

Could you give these a try.

[l975421700]

Hi Deniz

Thanks very much for the explanation. I tried to specify in the yaml file : jsbach>add_namelist_changes>namelist.jsbach>runctl>out_filetype: 2. But I still cannot get *01_wiso.nc output.

I still get the following in the log file:

0: jsbach_config: JSBACH output file type: 1
0: jsbach_config: JSBACH output file compression: 1

The experiment is installed here: /work/ollie/qigao001/output/echam-6.3.05p2-wiso/pi/pi_echam6_1d_029_3.38.

[denizural]

Hi Qinggang,

Here is an excerpt from namelist.jsbach:

&jsbach_ctl
    standalone = .false.
    ntiles = 11
    use_bethy = .true.
    use_phenology = .true.
    use_albedo = .true.
    with_yasso = .true.
    with_hd = .false.
    use_roughness_lai = .true.
    use_roughness_oro = .false.
    veg_at_1200 = .false.
    file_type = 1
    file_ztype = 1
    lpost_echam = .false.
    debug = .false.
/

As you can see unlike ECHAM, JSBACH does not have a runctl section. Sorry for misleading you in the previous message. It should be something like

jsbach:
    add_namelist_changes:
        namelist.jsbach:
            file_type: 2

Could you please try that.

[l975421700]

Thanks

When I execute esm_runscripts /work/ollie/qigao001/a_basic_analysis/d_model_simulations/2_echam6/1.1_pi/pi_echam6_1d.yaml -e pi_echam6_1d_029_3.38, I got the following error:

Traceback (most recent call last):
File "/home/ollie/qigao001/.local/bin/esm_runscripts", line 30, in
sys.exit(load_entry_point('esm-tools', 'console_scripts', 'esm_runscripts')())
File "/home/ollie/qigao001/esm_tools/src/esm_runscripts/cli.py", line 268, in main
Setup()
File "/home/ollie/qigao001/esm_tools/src/esm_runscripts/sim_objects.py", line 69, in call
self.prepcompute()
File "/home/ollie/qigao001/esm_tools/src/esm_runscripts/sim_objects.py", line 161, in prepcompute
self.config = prepcompute.run_job(self.config)
File "/home/ollie/qigao001/esm_tools/src/esm_runscripts/prepcompute.py", line 37, in run_job
config = evaluate(config, "prepcompute", "prepcompute_recipe")
File "/home/ollie/qigao001/esm_tools/src/esm_runscripts/helpers.py", line 62, in evaluate
framework_recipe, framework_plugins, config
File "/home/ollie/qigao001/esm_tools/src/esm_plugin_manager/esm_plugin_manager.py", line 141, in work_through_recipe
config = getattr(submodule, workitem)(config)
File "/home/ollie/qigao001/esm_tools/src/esm_runscripts/prepcompute.py", line 186, in modify_namelists
config[model] = Namelist.nmls_remove(config[model])
File "/home/ollie/qigao001/esm_tools/src/esm_runscripts/namelists.py", line 132, in nmls_remove
for key in list(change_entries):
TypeError: 'int' object is not iterable

[denizural]

Hi Qinggang,
Sorry for the typo, jsbach_ctl is missing there. I just checked it in my workdir and it worked. So

add_namelist_changes:
    namelist.jsbach:
        jsbach_ctl:    # this was missing
            file_type: 2

[l975421700]

Hi Deniz

Thanks. Sorry to say it still does not work. The experiment is installed in the directory: /work/ollie/qigao001/output/echam-6.3.05p2-wiso/pi. I run three experiments to test:

1. pi_echam6_1d_032_3.39: I specified namelist.jsbach: jsbach_ctl: file_type: 2 as you suggested. There are some NETCDF outputs, but no meaningful variables are output (only dimensions).
2. pi_echam6_1d_033_3.39: I removed the above settings and confirmed that my yaml file and executable could output GRIB files.
3. pi_echam6_1d_034_3.39: I added file_ztype: 2 to the first experiment. It does not change anything.

Do you have any idea where is the problem?

Besides, assuming it works, will it give me float 64 variables in .nc files if I set 'bits=24' for the GRIB output.

Regards, Qinggang

[denizural]

Hi Qinggang,

I checked your experiment directory (/work/ollie/qigao001/output/echam-6.3.05p2-wiso/pi/pi_echam6_1d_034_3.39) and here are the wiso files:

$ find . -name "*wiso*nc" -exec du -h {} \;
3.5K    ./run_20000101-20000131/work/pi_echam6_1d_034_3.39_200001.01_js_wiso.nc
512     ./run_20000101-20000131/work/pi_echam6_1d_034_3.39_200001.01_sf_wiso.nc
3.5K    ./run_20000101-20000131/work/pi_echam6_1d_034_3.39_200001.01_accw_wiso.nc
3.5K    ./run_20000101-20000131/work/pi_echam6_1d_034_3.39_200001.01_la_wiso.nc
3.5K    ./run_20000101-20000131/work/pi_echam6_1d_034_3.39_200001.01_wiso.nc
3.5K    ./unknown/pi_echam6_1d_034_3.39_200001.01_js_wiso.nc
512     ./unknown/pi_echam6_1d_034_3.39_200001.01_sf_wiso.nc
3.5K    ./unknown/pi_echam6_1d_034_3.39_200001.01_accw_wiso.nc
3.5K    ./unknown/pi_echam6_1d_034_3.39_200001.01_la_wiso.nc
3.5K    ./unknown/pi_echam6_1d_034_3.39_200001.01_wiso.nc

As you said, they only contain dimensions but not actual output data.
I see that the JSBACH module responsible for writing the output data is ./mo_soil.f90.
I am responsible for the wiso development but maybe you can have a look at that module to see why the data is not written.
@pgierz, do you have any suggestion?

[l975421700]

Hi Deniz

Thanks a lot. I will look into it tonight and tomorrow. Possibly with Martin in Cambridge. If you get any hints in the meantime, please let me know.

[l975421700]

Hi Deniz

I just found that if I run the simulation for one month rather than one day, I can get netcdf output in the unknown folder for echam and wiso streams. The simulation is here: /work/ollie/qigao001/output/echam-6.3.05p2-wiso/pi/pi_echam6_1m_114_3.46. Do you have any idea about why? For development, I indeed prefer to run 1 day simulation with NC output.

In addition, DTYPE of variables in the wiso stream is always F32, even I set "bits=24" in the "mo_memory_wiso.f90" for some variables. How can I get FLOAT64 for NETCDF with esm_tools?

[denizural]

Hi Qinggang,

I also realized that they are moved to the unknown directory. This happens when ESM-Tools doesn't know what to do with that file type. This behavior can be modified in the echam.yaml file. After each simulation, ESM-Tools moves the files in the work/ directory to their corresponding target directories (eg. restarts, outdata, ...)

In addition, DTYPE of variables in the wiso stream is always F32, even I set "bits=24" in the "mo_memory_wiso.f90" for some variables. How can I get FLOAT64 for NETCDF with esm_tools?

Unfortunately, that has nothing to do with the memory allocated by the computer. The bits variable there is probably related to the GRIB format. As you can see in mo_memory_wiso.f90

REAL(dp), ALLOCATABLE, PUBLIC :: wisoqte   (:,:,:,:)
REAL(dp), ALLOCATABLE, PUBLIC :: wisoxlte  (:,:,:,:)
REAL(dp), ALLOCATABLE, PUBLIC :: wisoxite  (:,:,:,:)
...
ALLOCATE (wisoqte  (nproma,nlev,nwiso,ngpblks)) ;wisoqte  = 0.0_dp
ALLOCATE (wisoxlte (nproma,nlev,nwiso,ngpblks)) ;wisoxlte = 0.0_dp
ALLOCATE (wisoxite (nproma,nlev,nwiso,ngpblks)) ;wisoxite = 0.0_dp

These are always double precision floating point numbers. The strange thing is that netcdf show them as float instead of double.
I need to leave office now but I will edit this message later on

[denizural]

Hi Qinggang,
could you please show me which files / variables have single precision floats? I check the NetCDF files in your experiment and they are all double. Two things are important:

• How the variables are declared and allocated in Fortran. Eg. real(kind=dp) :: wiso_variable
• How NetCDF declares and writes them to the files via: nf_var_def (or nf90_var_def). It is important that these functions use NF_DOUBLE (or NF90_DOUBLE) for the double precision.

You can always verify this with the ncdump command. Eg. ncdump -h <myfile.nc>

[l975421700]

Hi Deniz

Thanks. I used the command cdo sinfo /work/ollie/qigao001/output/echam-6.3.05p2-wiso/pi/pi_echam6_1m_114_3.46/unknown/pi_echam6_1m_114_3.46_200001.01_wiso.nc. All vairbales have Dtype F32. This file is under the unknown directory.

If I use the command cdo sinfo /work/ollie/qigao001/output/echam-6.3.05p2-wiso/pi/pi_echam6_1m_114_3.46/outdata/echam/pi_echam6_1m_114_3.46_200001.01_accw_wiso. Variables can have Dtype P24s or P16s. What does this two means for NETCDF file? And why it is different for files in different directories.

The variables are always declared as dp/wp, i.e. double precision, in the model.

I have no idea yet how the model write variables into NETCDF file. Which subroutine shall I search?

Thanks for your time.

[denizural]

Hi Qinggang,
To query variables cdo sinfo works but sometimes it fails. ncdump -h always works. When the datatype is float it is single-precision (SP) file and double is double-precision (DP) floating point number.

I also verified that the variables in the unknown directory are single-precision files. Normally, the NetCDF functions to look for are nf_var_def / nf90_vardef and nf_put_var / nf90_put_var. ECHAM has few layers of abstraction over these low-level IO functions. Eg. IO_DEF_VAR (defined in mo_io.f90 calls that. You can also check that this subroutine writes the coordinates in NF_DOUBLE:

./mo_io.f90:386:          CALL IO_DEF_VAR (fileID, io_name, NF_DOUBLE, ndim, &
./mo_io.f90:422:      CALL IO_DEF_VAR (fileID, 'lat', NF_DOUBLE, 1, IO_dims, io_vlat_id)
./mo_io.f90:424:      CALL IO_DEF_VAR (fileID, 'lon', NF_DOUBLE, 1, IO_dims, io_vlon_id)
./mo_io.f90:428:      CALL IO_DEF_VAR (fileID, 'landpoint', NF_DOUBLE, 1, IO_dims, io_vland_id)
./mo_io.f90:433:      CALL IO_DEF_VAR (fileID, 'vct_a', NF_DOUBLE, 1, IO_dims, io_vcta_id)
./mo_io.f90:434:      CALL IO_DEF_VAR (fileID, 'vct_b', NF_DOUBLE, 1, IO_dims, io_vctb_id)

The subroutine IO_PUT_VAR_DOUBLE is also responsible for the file IO. According to my observations the variables need to be written in DP but ncdump (or cdo sinfo) shows that variables (not the dimensions) are written in SP.

The strange thing is variables are written with DP to restart files but with SP to output files. It looks like restart files are written like this (call stack)

SUBROUTINE stepon
    CALL write_streams                     ! mo_io.f90
        CALL io_write_stream
            CALL IO_PUT_VAR_DOUBLE         ! mo_netcdf.f90
                CALL NF_PUT_VAR_DOUBLE     ! this is the actual NetCDF call    

I think the problem arises from CDI library which is used in file output:

SUBROUTINE stepon
    CALL out_streams                 ! mo_output.f90
        CALL out_stream
            CALL write_var
                CALL StreamWriteVar  ! comes from CDI library

CDI documentation says that it writes data in DP: https://code.mpimet.mpg.de/projects/cdi/embedded/cdi_fman.html#x1-510004.1.11

You can start your debugging by placing print statements before call to StreamWriteVar. Here are the lines from mo_output.f90:

1676         IF (nprocio == 0) THEN
1677           CALL StreamWriteVar(fileID, varID, xyz(:,:,:,1), nmiss)
1678         ELSE
1679 #ifdef HAVE_CDIPIO
1680           !PRINT *, 'pe ', p_pe, ', varid=', varid, 'size(xyz) = ', size(xyz), &
1681           !     'xt_idxlist_get_num_indices = ', xt_idxlist_get_num_indices(gp(SIZE(xyz,3))%idxlist), &
1682           !     'dim=', info%dim, SIZE(xyz,3)
1683           CALL streamWriteVarPart(fileID, varID, xyz(:,:,:,1), nmiss, &
1684           gp(SIZE(xyz,3))%idxlist)
1685 #endif
1686         ENDIF

[l975421700]

Hi Deniz

Thanks a lot for your efforts. I also checked nf_def_var , IO_DEF_VAR , IO_PUT_VAR_DOUBLE , NF_PUT_VAR_DOUBLE , but found no clues.

I can find the lines in mo_output.f90. But may I ask, what shall I print specifically? This part of code is totally new to me.

[denizural]

Hi Qinggang,
these lines write the data to the disk so you have the opportunity to observe them before they are written to files. I am also not an ECHAM expert 🙂 However, this gives you a good starting point for debugging the model. Ie. You can dump the data to a text file and to netcdf and compare them if there are significant differences.

I also have a short NetCDF IO program on Ollie: /home/ollie/dural/programming/netcdf/sfc_pres_temp_wr.f90. You can see how NetCDF is written there.

I hope you find these helpful. Sorry to tell you that although we work with ECHAM, we are not ECHAM developers so deep technical help on ECHAM is beyond your assistance.

[l975421700]

Hi Deniz

Thanks a lot, those are helpful suggestions. I will have a look into it.

[pgierz]

Moin all,

@pgierz, do you have any suggestion?

Yes! So a few things:

• This is very much an echam internal problem, so to the question:

How can I get FLOAT64 for NETCDF with esm_tools.

All esm-tools really does, if you boil it down, is copy around input and output files and fire off a batch system job for you. It's complicated because we try to support many different HPC systems and models, but output routines from ECHAM has nothing to do specifically here with respect to using esm-tools or some home-made scripts.

• Grb vs NetCDF: I am not quite sure if that will solve your problem. If the array is already not conserving mass or energy in the program memory, changing the output format will not really help you to reveal that problem. A better idea in my opinion would be to check that before you even flush the array to disk: do a global sum, subtract, and see if the answer is exactly zero, and have the model crash if that is not the case. If it is already wrong in memory, it'll be wrong on disk as well, and seeing which format it is saved to will not help you

• If CDI is causing you headaches, the best solution would be to ask at the MPI-Met forum. CDI is effectively just a Fortran and C interface to the same calculations that happen in CDO, and is proprietary to MPI-Met. We unfortunately cannot help there.