#!/bin/ksh

# ICON
#
# ------------------------------------------
# Copyright (C) 2004-2026, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
# Contact information: icon-model.org
# See AUTHORS.TXT for a list of authors
# See LICENSES/ for license information
# SPDX-License-Identifier: BSD-3-Clause
# ------------------------------------------

# -------------------------------------------------------
# MARS request to read SST or SKT climatology from ERA40
#(climatology of monthly means of daily means),
# H_SNOW and RHO_SNOW. Calculate W_SNOW from h_snow*rho_snow*(lsm>=0.5)
# Interpolate to ICON grid.
# Instead of SST it is possible to read skin temperature, SKT.
#
# Usage:  mars4icon_e4_sst "months" sst|skt
#
# attention:
#   * at DWD requires ECMWF login by eccert
#   * requires atleast mars version 1.9.9
#
# history:
#  Helmut Frank
# -------------------------------------------------------

# Executable to do horizontal interpolation
#iconremap=~fprill/remap_cons/icontools/iconremap
iconremap=~for0exp/ICON/iconremap

# netCDF file with ICON grid
icongrid_dir=~for0exp/ICON/Testdata/ICONGPI_example
#icon_grid_file=${icongrid_dir}/extpar_${icongrid}.nc
icon_grid_file=${icongrid_dir}/input_grid.nc

# define data directory (if not set otherwise)
datadir=${SSTICEDIR:-${SCRATCH}/icon}
if [[ ! -d ${datadir} ]]; then mkdir -p $datadir; fi

typeset -Z2 month
month_list="${1:-01 02 03 04 05 06 07 08 09 10 11 12}"
temp=${2:-skt}   # sst = sea surface temperatur, or skt = skin temperature

typeset -u Temp=${temp}
typeset t_code=34
case ${Temp} in
    SKT) t_code=235 ;;
    SST) t_code=34  ;;
esac

set -A month_name null jan feb mar apr may jun jul aug sep oct noc dec

typeset mdate ifs_data icon_data mars_req

for month in ${month_list}
do
    mdate=${month_name[$month]}
    ifs_data=${datadir}/ifs_e4_${Temp}_CL${month}
    case ${temp} in
        skt) icon_data=${datadir}/icon_cl${month}_TS_Wsnow ;;
          *) icon_data=${datadir}/icon_cl${month}_Ts_Wsnow ;;
    esac

    if [[ -s ${ifs_data}.grb ]]; then
        print -- "ERA40 data already present:"
#       ls -l ${ifs_data}.grb
#       grib_list -v -stat ${ifs_data}.grb
    else

        mars_req=${TMPDIR}/MARS_IN.$$
        cat > ${mars_req} << MARS_IFS2ICON_EOF.$$
# MARS request for ifs2icon
retrieve,
        class    = e4,                                  # ERA40
        expver   = 1,
        domain   = g,
        stream   = moda,
        type     = cl,                                  # monthly mean over several years
        date     = ${mdate},
        target   = "${ifs_data}.tmp",
        param    = ${temp},                             # sst or skt
        repres   = gg,                                  # gaussian grid
        gaussian = regular,
        grid     = 80,                                  # for ERA40
        levtype  = surface
retrieve,
        param = sd, fieldset = hsnow                    # snow depth
retrieve,
        param = rsn, fieldset = rhosnow                 # snow density
retrieve,
        param = lsm, fieldset = lsm                     # land sea mask
compute,fieldset=wsnow, formula="hsnow*rhosnow*(lsm>=0.5)"
write,  fieldset=wsnow,target="${ifs_data}.tmp"
MARS_IFS2ICON_EOF.$$

        mars ${mars_req}
        rc_mars=$?

        if [ $rc_mars -ne 0 ]; then
            { print -- "Error $rc_mars executing mars"
              print -- "Input file to MARS:"
              cat ${mars_req}
            } >&2
             rm ${mars_req}
            exit $rc_mars
        fi

        print "\nRetrieved ERA40 analysis for ${mdate} from MARS"
        ls -al ${ifs_data}.tmp
        print -- "count,dataDate,dataTime:s,step,stepUnits:s,shortName,level,typeOfLevel,average,min,max"
        grib_get -P count,dataDate,dataTime:s,step,stepUnits:s,shortName,level,typeOfLevel,average,min,max ${ifs_data}.tmp

# Replace date by DWD date for climatological data
        grib_set -s dataDate=1111${month}11,dataTime=1100,centre=edzw ${ifs_data}.tmp ${ifs_data}.grb
        rm -f ${ifs_data}.tmp ${mars_req}

    fi

    ls -l ${ifs_data}.grb
    print -- "count,dataDate,dataTime:s,step,stepUnits:s,shortName,level,typeOfLevel,average,min,max"
    grib_get -P count,dataDate,dataTime:s,step,stepUnits:s,shortName,level,typeOfLevel,average,min,max ${ifs_data}.grb

#
# Interpolate to ICON grid
#
    nml_remap=${TMPDIR}/sst2icon_remap.nml
    nml_ifld=${TMPDIR}/sst2icon_infield.nml

cat > ${nml_remap} <<EOF_NML
! Step 1: interpolation regular grid-> ICON
&remap_nml
 in_grid_filename  = "${ifs_data}.grb"
 in_filename       = "${ifs_data}.grb"
 in_type           = 1
 out_grid_filename = "${icon_grid_file}"
 out_filename      = "${icon_data}.tmp"
 out_type          = 2
 out_filetype      = 2
 lcompute_vn = .FALSE.
 lcompute_vt = .FALSE.
/
EOF_NML

cat > ${nml_ifld} <<EOF_NML
! DEFINITIONS FOR IFS INPUT DATA
!
&input_field_nml  ! temperature
 inputname      = "${Temp}"
 outputname     = "T_S"
 code           = ${t_code}
 type_of_layer  = "surface"
 missval        = 270
/
&input_field_nml  ! snow
 inputname      = "sd"
 outputname     = "H_SNOW"
 code           = 141
 type_of_layer  = "surface"
/
EOF_NML

    print -- "\nMap ${ifs_data}.grb to ICON grid ${icon_grid_file}"
    ${iconremap} --remap_nml=${nml_remap} --input_field_nml=${nml_ifld} -v
#
#   Rename H_SNOW (from sd in MARS) to W_SNOW
#
    rules_file=${TMPDIR}/rules4ifs2icon.$$
    set -A datum $(date +'%Y %m %d %H %M %S')
    cat > ${rules_file} <<REOF
if ( shortName is "H_SNOW") {
    shortName = "W_SNOW";
}
 localDefinitionNumber=254;
 grib2LocalSectionPresent =1;
 setLocalDefinition =1;
 localCreationDateYear = ${datum[0]};
 localCreationDateMonth = ${datum[1]};
 localCreationDateDay = ${datum[2]};
 localCreationDateHour = ${datum[3]};
 localCreationDateMinute = ${datum[4]};
 localCreationDateSecond = ${datum[5]} ;

REOF

#   grib_set -w shortName=H_SNOW -s shortName=W_SNOW ${icon_data}.tmp ${icon_data}
    grib_convert ${rules_file} ${icon_data}.tmp ${icon_data}
    rm ${icon_data}.tmp ${rules_file}

    print -- "\nClimatological data interpolated to ICON grid ${icon_grid_file}"
    ls -l ${icon_data}
#   grib_list -v -stat ${icon_data}
    print -- "count,dataDate,dataTime:s,step,stepUnits:s,shortName,level,typeOfLevel,average,min,max"
    grib_get -P count,dataDate,dataTime:s,step,stepUnits:s,shortName,level,typeOfLevel,average,min,max ${icon_data}

done

exit 0
