README.md

Ingredients to make a reproducibility test of EC-Earth 3.1 and later.
Largely based on Paolo Davini's initial scripts to make a climate assessment of EC-Earth

Current limitations and constraints:
* The output must come from EC-Earth 3.1 or later
* The output must be available as MMO, MMA tar files and ICMGG, ICMSH grib files
* The output must consist in 12 month data per file
* You have to use the classic outclass (https://earth.bsc.es/gitlab/es/auto-ecearth3/-/issues/1547)
* The simulation must start with January and finish with December
* Run at least 5 members for each experiment
* The scripts are adapted to compare two and only two experiments at time (in `compare.sh`)

In proj.conf, set:
* NFIXYR=`"YEAR OF YOUR SIMULATION (YYYY)"`
* OCEAN_ini=official_repro
* OCEAN_ini_member=all_members
* ICE_ini=official_repro
* ICE_ini_member=all_members
* ATM_ini=official_repro
* ATM_ini_member=all_members
* ATM_ini_member_perturb=FALSE
* SAVEMMA=TRUE
* SAVEMMO=TRUE

# Checkout the git project from the BSC archive

```sh
 git clone https://earth.bsc.es/gitlab/es/reproducibility-test.git
 cd reproducibility-test
 git checkout master
```

# Edit the configuration file

```sh
 vi ./scripts/config.sh
```
setting variables:
- STARTDATE=`"YOUR START_DATE (YYYYMMDD)"`
- ECVERSION=`"3.3.1 OF LATTER"`
- maskfile=${ece3_constant_files}/ifs/`"RESOLUTION"`/ICMGGECE3INIT


# Run the reproducibility test
## A) Use the full_loop script
If you feel lucky, set and run the `full_repro_loop.sh` script
```sh
 ./full_repro_loop a3ko a3kp 1950 1950 5
```

## B) Go step by step
# Pre-process the EC-Earth output

```sh
cd ./scripts
./preprocess.sh ${exp} ${member} ${yearb} ${yeare}
```

[ Execute ./preprocess.sh to see the help ]

Process *all* your members of all experiments. This can be done in paralle
by using nohup commands for instance:

```sh
exp=m06e;
for i in 0 1 2 3 4 ;  do 
  sleep 2; 
  ./preprocess.sh ${exp} fc${i} 1850 1869 >& log/log_prpr_${exp}_fc${i} 
done
```

# Re-grid to the same grid as the observational reference

```
for member in 0 1 2 3 4 ;  do
  sleep 2;
  nohup ./postprocess.sh ${exp} fc${member} ${yearb} ${yeare}
done

```

# Compute the Reichler and Kim indices

```sh
./scores.sh ${exp} ${member} ${yearb} ${yeare}
```

NetCDF files and a text file will be available in the directory EC-Earth_diags.
This text file (PI3.txt) consists in scores. By looking at how scores vary from
member to member one can guess if an experiment is reproducible wrt another one.

Once steps 1/ -> 5/ have been done for all members, we can proceed
to ensemble analyses:

# Ensemble analyses

```sh
./ensemble.sh ${exp} ${yearb} ${yeare}
```

Process all your experiments with this script. And then produce the plots comparing the Reichler and Kim indices and the variables differences between two experiments:

# Run reproducibility diagnostics between two experiments

```sh
./compare.sh ${exp1} ${exp2} ${year1} ${year2} ${number_of_members}
```

This script launch three R scripts:

1. Create the plots for basic time series (2-m temperature, sea ice extent, precip)
2. The second on make the comparison based on the Reichler-Kim indices with the KS test (script R_scripts/KS_index.R). It produces the plot reichler_kim_scores_stat_${exp1}_${exp2}.eps
3. The third script produces the maps of the differences, with a KS test applied in each grid cell (script map_diff_experiments.R), for the variables 't2m','msl','qnet','tp','ewss','nsss' and 'SICE'. RData files are also saved to redo quickly the maps, changing for example the levels of the colour bars.

The results will be in scripts/EC-Earth_diags. 
If you prefere the plots to be in **pdf** format, just run the convert_and_merge_to_pdf script:

```sh
./convert_and_merge_to_pdf ${exp1} ${exp2}
```

Note that there is no proper way to handle fc?? naming for members (it should work fine with fc?). If you are using it,  you should do some changes in the R_scripts:

- R_scripts/basic_plots.R :
```
 -@@ -48,7 +48,7 @@ tp  <- array(NA, dim = c(n.exp,         n.memb, n.year, n.month))
   for (j.memb in seq(1, n.memb)) {
-    fcmemb <- paste("fc", j.memb - 1, sep = "")
+    fcmemb <- paste("fc0", j.memb - 1, sep = "")
     print(paste("  ", fcmemb))
     # SEA ICE
```

- R_scripts/map_diff_experiments.R :
```
@@ -41,7 +41,7 @@ for (jexp in 1:2) {
   for (jmemb in 1:nmemb) {

     memb=jmemb-1
-    fnc <- open.ncdf(paste(path,'/',exp_name[jexp],'/fc',memb,'/',var,'_mon_2x2.nc',sep=""))
+    fnc <- open.ncdf(paste(path,'/',exp_name[jexp],'/fc0',memb,'/',var,'_mon_2x2.nc',sep=""))
     lon <- get.var.ncdf(fnc,'longitude')
```