Hooking extreme event detection algorithms into esm_runscripts

[JanStreffing]

At the palo climate retreat we had a discussion on automatically triggering some of our existing extreme event detection algorithms for AWI-ESM2, AWI-CM3, and eventually AWI-ESM3 output. People that would be interested to hook their existing scripts into esm_tools were: @norelrimbu and @justuscontzen after his thesis.

Questions and suggestions from my side:

• Do we ensure that the detection is used only on runs that are long enough for them to make sense, or do we just trust the user to make an informed choice here?
• How do we to ensure the variables need are available? / How do we make sure we get warnings not errors if they are not?
• How are your scripts configured at the moment? If it is via setting variables in the script itself, we would first need to make modifications to either A) read command lines arguments or B) read directly from a yaml file.

I would suggest we try this process with a one or two frequently used algorithm and see how it goes.

@nwieters @mandresm @pgierz From your side it would be nice to have some input on which place to hook in. e.g. do we use the general postprocessing script interface or some place else. Also, do we better control via arguments or is it easier to read straight from the finished yaml.

[mandresm]

Hi Jan,

Depending on how expensive these algorithms are in terms of resources, one would need to use 1) the workflow manager to launch an independent job in slurm, or 2) make a plugin for it to connect to the same job as the simulation, the advantage being that the plugin allows coding in Python and accessing the information in config.

I don't think we should add by default extra stuff to AWICM3 and AWIESM2 (I think the user should make sure they trigger it themselves), however, I don't really understand the science behind it or the needs for it. As always, it's more of a matter of what the users want/need.

[norelrimbu]

Hello all,

As the background of people interested in this project is very different, would be useful to clarify first what we want to deliver at the end. We can create climate extreme index data sets with monthly, seasonal or annual resolution from different model simulations. They can be stored and used further for different analysis. We can create also tools (operators similar to cdo) to be used for extreme index calculations by the users when high resolution data, i.e. daily or higher, are available.
For understanding extreme climate variability, atmospheric blocking is essential. There are some tools on the market for blocking calculation , e.g. https://docs.esmvaltool.org/en/latest/recipes/recipe_miles.html, that can be used for comparison.
I propose to start with the calculation of one-dimensional Tibaldi-Molten,(TM) blocking index. It is largely used for blocking-frequency calculation. More important, the definition is very simple and only one field, i.e. daily or higher resolution 500 hPa geopotential height, is need as the input. I can provide an R function to be checked and adapted to be used in a complex computing environment.

Best,

Norel

[JanStreffing]

Hello Norel,
My understanding is that the goal of this effort here would be to adapt the methods that are already being used by you and @justuscontzen, so that they can be controlled and called by esm_tools tools as a postprocessing. I started this discussion, hoping to clarify:

1. which of your methods is most promising for hooking into esm_tools, both from a utility and ease of implementation perspective
2. get some input from the esm_tools team about the technical realization / where & how we should try this.

[norelrimbu]

Hello Jan,

OK, I understood perfectly well the goal. I wrote a lot of functions in R for analysis of extremes. Some of them are stored on /isibhv/projects/paleo_work/EXTREMES-NOREL/. You can have a short look , if you have time. I test these functions on reanalysis data, NCEP1, 20CRv2 and 20CRV3. It seems that they work properly. However, these functions should be reevaluated and optimized from a computing efficiency perspective.
So, from my point of view, we should start to 1-D blocking index calculation. If agree, I will explain, in detail the algorithm, as well as the associated R function. As a computing professionals, you can decide if my R function is appropriate or find a more efficient way to compute this index.
I will insert more comments into the R code and upload the function soon.

Regards,
Norel

[norelrimbu]

Sorry, the locations of the functions is /isibhv/projects/paleo_pool/EXTREMES-NOREL/

[norelrimbu]

Hello Jan,
As I mentioned before, I suggest starting with calculation of 1D blocking TM index. It is simple and largely used as a blocking frequency indicator. I added some more comments on my R function for TM index to be easily checked.
I attached the function as well as some comments that helps to understand the computing strategy. I am here for any further questions.
Best,
Norel

TM_index.pdf
bl01_1d.txt

[JanStreffing]

Good point. Expensive things should not stop us from resubmitting the next job. Anything longer than a minute or two should run in a separate job. @norelrimbu suggested having a look at https://github.com/oloapinivad/MiLES 1D and 2D blocking. I like that idea, as I used that package myself for a paper a while ago. It's written in R, but even as an R novice it was usable.

I agree we should not use the esm_tools backend, just the options to hook in.

[pgierz]

I would go for the route of new slurm jobs via the workflow manager. In principle, post jobs could get more and more expensive the more stuff we add in over time. In my view it should be:

main job --> tidy --> main job: depends on tidy (on slurm)
                  --> post job: depends on tidy (on slurm)
                               --> viz job (on slurm): depends on post
                  --> archive job: depends on tidy 

Diagramming is better in Gitlab 😉. I think clever use of the workflow manager can handle that.

@JanStreffing, what exactly do you mean here:

I agree we should not use the esm_tools backend, just the options to hook in.

[JanStreffing]

In the past some functionality was put directly into the esm_tools python code. e.g. postprocessing oasis EXPOUT files to generate oasis restart files. I don't think we should do this here. I though this is what Miguel mean here:

I don't think we should add by default extra stuff to AWICM3 and AWIESM2 (I think the user should make sure they trigger it themselves),

Though reading it again, he could have also meant that things should be behind choose blocks.

[mandresm]

yeap, I would suggest a choose block, and then plug in whatever and however ones wants to plug in something

[JanStreffing]

@mandresm, @norelrimbu. I think we have collected enough information to give it a first try. Maybe we can find a 2 hour time slot somewhere, meet at AWI and see how far we can get.

To give a bit of structure to our kickoff meeting, here are a few points we may want to consider. Feel free to add more:

• How do we want to call the script. Likely as separate postprocessing job. Already discussed a bit above.
• When we do this via a postprocessing job, do we make one job per script or one job for all scripts?
• Which variables are configurable in the script?
• Do we want to pass these variables as arguments or have the script read the required info from the finished yaml file?
• For testing purposes we need to be able to run this particular step of the esm_tools chain. I think this is easy for Migule, but it has been years for me that I used this feature.

I have no major meetings etc. the next weeks. I leave the suggestion of a date and time to you.

[pgierz]

I'm in for this.

How about the coming Monday (10th July) at 13:30, after lunch break?

[pgierz]

One point to add: If we use an external language (e.g. R) are there any external libraries we need to consider? How do we ensure these are available? (In practice this is a trivial thing, but still something we should at least keep in mind)

[mandresm]

Next week only really works for me on Thursday

[mandresm]

My point to add to Jan's bullet points:

• How many people is going to use this and how often

[norelrimbu]

For my R function only ncdf4 R-package is necessary. But for more complex calculations, like PV blocking based indices, we need more R packages. Keep in mind that new version of the R packages are published often. It is not difficult, probably, to change the R codes in Fortran routines.

[pgierz]

Yes, running generic shell commands (which we could also use to manipulate R) directly from FORTRAN is simple, e.g.:

program HelloWorld
    implicit none
    character(len=100) :: command
    integer :: exit_status

    ! Print "Hello world" from Fortran
    print *, "Hello world"

    ! Run external command using the shell
    command = "echo 'Hello world'"
    exit_status = execute_command_line(command)

    if (exit_status /= 0) then
        print *, "Error executing the command:", trim(command)
    end if

contains

    function execute_command_line(command) result(exit_status)
        character(len=*), intent(in) :: command
        integer :: exit_status

        call system(command, exit_status)
    end function execute_command_line

end program HelloWorld

[mandresm]

Some documentation about it:
https://github.com/esm-tools/workshops/blob/main/202207_AWI_ESM-Tools/exercises/7.5_workflow_manager.md

[justuscontzen]

The R scripts I would like to add also use external libraries. Until now, I ran them by using a conda environment where I had installed the packages before:

#!/bin/bash
#SBATCH [...]

module load conda
source activate myEnv
srun Rscript c3_locest_merge.R $1 $2 $3 $SLURM_ARRAY_TASK_ID

Would something similar be possible also when running the scripts as part of esm_tools?

[pgierz]

In principle yes, the hooks allow you to run any arbitrary shell script. We might need to fiddle a bit to get the requirements automatically installed, but that should not be a show stopper in my view.

It would be useful if you have a environment definition yaml for your conda

[justuscontzen]

environment.zip
I don't need all the R libraries installed there, I will check which ones exactly are required.