coupling_ice2echam.functions

#!/usr/bin/ksh

## @fn ice2echam() 
ice2echam() {
        echo "                =================================================================="
        echo "                 *** S T A R T I N G    ice2echam *** "; echo

        . ${FUNCTION_PATH}/../general/general_lists.functions
        . ${FUNCTION_PATH}/../general/coupling_general.functions

        ISM_TO_ECHAM_update_orography=${ISM_TO_ECHAM_update_orography:-1}
        ISM_TO_ECHAM_update_glacial_mask=${ISM_TO_ECHAM_update_glacial_mask:-1}
        ISM_TO_ECHAM_update_land_runoff=${ISM_TO_ECHAM_update_land_runoff:-1}
        HOSING_CORRECTION=${HOSING_CORRECTION:-0}

        :> nco_stderr_ice2echam; :> cdo_stderr_ice2echam
        
        counter=0
        COUNT_MAX=60
        while [ ${counter} -lt ${COUNT_MAX} ]
        do
            echo "          * inside while loop "
            if [ -f ${COUPLE_DIR}/ice_file_for_atmosphere.nc ]; then
                break
            fi

            echo; echo "        * File ${COUPLE_DIR}/ice_file_for_atmosphere.nc not found. Waiting for 10 seconds ..."
            sleep 10
            counter=$((counter+1))
        done
        
        read_names ice atmosphere
        define_echam_names
        if [ "$ISM_TO_ECHAM_update_glacial_mask" -eq 1 ]; then
                iterative_coupling_ice_echam6_update_glacial_mask
        fi
        if [ "$ISM_TO_ECHAM_update_orography" -eq 1 ]; then
                iterative_coupling_ice_echam6_update_orography
        fi
        if [ "$ISM_TO_ECHAM_update_land_runoff" -eq 1 ]; then
                iterative_coupling_ice_echam6_update_land_runoff
        fi

        echo "                 *** F I N I S H E D    ice2echam *** "
        echo "                =================================================================="
}

## @fn define_echam_names() 
define_echam_names() { 
        echam_glacial_mask_name=glac
}

## @fn iterative_coupling_ice_echam6_update_orography() 
iterative_coupling_ice_echam6_update_orography() { 
        echo "                Updating orography and gravity wave drag parameterization in echam6"
        update_orography_select_field
        if [ "$MACHINE" == "levante" ]; then
            hi_res_init_oro="/pool/data/ECHAM5/T511/T511_jan_surf.nc"
        fi
        if [ -f ${hi_res_init_oro} ]
        then
            regrid_to_echam ${COUPLE_DIR}/ice_orog_difference.nc remapbil T511
            regrid_to_echam ${COUPLE_DIR}/ice_orog_difference.nc remapbil 
            update_orography_apply_orography_anomaly_high_res
            update_orography_prepare_calnoro_high_res
        else 
            regrid_to_echam ${COUPLE_DIR}/ice_orog_difference.nc remapbil 
            update_orography_apply_orography_anomaly
            update_orography_prepare_calnoro
        fi
        update_orography_run_calnoro
        update_orography_generate_target_orography_file
        #set_in_jsbach_restart_file elevation ${COUPLE_DIR}/target_orography_echam6_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc OROMEA unpack
        update_orography_set_namelist_modifications_next_echam_run
        echo; echo "                 ...done."; echo
}

## @fn iterative_coupling_ice_echam6_update_glacial_mask() 
iterative_coupling_ice_echam6_update_glacial_mask() { 
        echo "                Updating Glacial Mask in echam6"
        set_glacial_mask_select_field
        set_glacial_mask_echam_restart_file $ifile
        if [ -d ${COUPLE_DIR}/jsbach_init_file_modifications ]; then
                rm -rf ${COUPLE_DIR}/jsbach_init_file_modifications
        fi
        mkdir ${COUPLE_DIR}/jsbach_init_file_modifications
        here=$(pwd)
        cd ${COUPLE_DIR}/jsbach_init_file_modifications
        set_glacial_mask_echam_make_dummy_jan_surf
        set_glacial_mask_jsbach_update_init_file
        cd $here
        rm -rf ${COUPLE_DIR}/jsbach_init_file_modifications
        set_glacial_mask_jsbach_update_restart_jsbach
        set_glacial_mask_jsbach_update_restart_veg
        ensure_echam_restart_aliases_for_m_streams
        ensure_nonempty_grib_rerun_files_for_resume_month
        
        # === FINAL FIX: Force lgfw=.false. in namelist.echam ===
        # Even with namelist modification files, lgfw sometimes persists as .true.
        # causing IO_read_stream failures for 'gfw'. We patch it directly here.
        if [ -n "${WORK_DIR}" ] && [ -f "${WORK_DIR}/namelist.echam" ]; then
            echo "                -   Patching ${WORK_DIR}/namelist.echam to force lgfw = .false."
            sed -i 's/lgfw[[:space:]]*=[[:space:]]*\.true\./lgfw = .false./Ig' "${WORK_DIR}/namelist.echam"
            grep -i "lgfw" "${WORK_DIR}/namelist.echam"
        fi
        
        echo; echo "                 ...done."; echo
}

## @fn iterative_coupling_ice_echam6_update_land_runoff() 
iterative_coupling_ice_echam6_update_land_runoff() { 
        echo "                Setting Glacial Runoff for incorporation in HD-Model via oasis3-mct"
        update_land_runoff_select_glacial_discharge
        regrid_to_echam ${COUPLE_DIR}/ice_discharge.nc
        update_land_runoff_fill_missing
        if [ "x${HOSING_CORRECTION}" == "x1" ]; then
            echo; echo "            * Correcting freshwater input by hosing "
            distribute_total_discharge_over_ocean
        else
            echo; echo "            * NOT correcting freshwater input by hosing "
        fi
        update_land_runoff_prepare_for_oasis
        update_land_runoff_set_namcouple_override
        update_land_runoff_set_namelist_modifications_next_echam_run
        update_land_runoff_set_namelist_modifications_next_jsbach_run
        echo; echo "                 ...done."; echo
}

## @fn update_orography_select_field() 
update_orography_select_field() {
        echo; echo "                *   selecting orography difference over ice run"
        cdo -s selvar,${ice_orography_difference_name} \
                ${COUPLE_DIR}/ice_file_for_atmosphere.nc \
                ${COUPLE_DIR}/ice_orog_difference.nc
        rmlist="$rmlist ${COUPLE_DIR}/ice_orog_difference.nc"
}

update_orography_from_martins_script() {
cdo -L -remapbil,${RES_echam}grid /
            ${COUPLE_DIR}/ice_orog_difference_T511grid.nc / 
            ${COUPLE_DIR}/ice_orog_difference_${RES_echam}grid.nc
#cdo -L -setmisstoc,0. -ifthen slf.nc ${COUPLE_DIR}/ice_orog_difference_${RES_echam}grid.nc dummy.nc
#mv dummy.nc ${COUPLE_DIR}/ice_orog_difference_${RES_echam}grid.nc
ECHAM_g=9.80665
if [ -L ${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc ]; then
        rfile=$(readlink ${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc)
else
        rfile=${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc
fi
cdo -L -sp2gp -gp2sp -add -selvar,GEOSP ${rfile} -mulc,${ECHAM_g} /
        ${COUPLE_DIR}/ice_orog_difference_${RES_echam}grid.nc /
        ${COUPLE_DIR}/echam6_new_orography_before_calnoro_low_res.nc
mv geosp.nc ${output_dir}/input/echam6/GEOSP_${paleo_time}_${RES_echam}.nc
}

update_oromea_from_martins_script() {
ofile=${COUPLE_DIR}/target_orography_echam6_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc

cdo -L -remapnn,r720x360 -setgrid,T511grid -selvar,OROMEA ${hi_res_init_oro} oromea_pi_r720x360.nc
cdo -L -remapnn,r720x360 ${COUPLE_DIR}/ice_orog_difference_T511grid.nc anom_topo_r720x360.nc
#cdo -L -setmisstoc,0. -ifthen -remapnn,r720x360grid slf.nc anom_topo_r720x360.nc dummy.nc
#mv dummy.nc anom_topo_r720x360.nc
cdo -L -add oromea_pi_r720x360.nc anom_topo_r720x360.nc oromea_r720x360.nc
cdo -L -mul oromea_r720x360.nc -gtc,0. oromea_r720x360.nc dummy.nc; mv dummy.nc oromea_r720x360.nc
CALNORO_RES=${echam_res:1:3}                       # use defined ECHAM6 grid resolution for calnoro
cp oromea_r720x360.nc topodata.nc
echo ${CALNORO_RES} | ${calnoro_dir}/calnoro
cdo -L -s -f nc -setgrid,${echam_res}grid -invertlat sso_par_fil.srv calnoro_output.nc
cdo -L -chvar,var51,OROMEA,var52,OROSTD,var53,OROSIG,var54,OROGAM,var55,OROTHE,var56,OROPIC,var57,OROVAL calnoro_output.nc ${ofile}
cdo -L -mul -selvar,OROMEA ${ofile} -gtc,0. -selvar,OROMEA ${ofile} oromea.tmp
# save geosp and oro-xxx in different files. why oromea from original file @PG?
cdo -L -replace ${ofile} oromea.tmp dummy.nc; mv dummy.nc ${ofile}
rm oromea.tmp
rm oromea_pi_r720x360.nc anom_topo_r720x360.nc oromea_r720x360.nc topodata.nc sso_par_fil.srv calnoro_output.nc
}


## @fn update_orography_apply_orography_anomaly_high_res() 
update_orography_apply_orography_anomaly_high_res() {
        echo; echo "                *   applying surface height changes to high-resolution echam topography"
        ifile=${COUPLE_DIR}/ice_orog_difference_T511grid.nc
        ofile=${COUPLE_DIR}/echam6_new_orography_before_calnoro_high_res.nc
        ncrename -v ${ice_orography_difference_name},oromea $ifile
        cdo -s -setmisstoc,0 $ifile tmp; mv tmp $ifile
        if [[ -f ${ofile} ]]
        then
            rfile=${ofile}
            ncks -O -v oromea $rfile tmp_oromea
        else
            rfile=${hi_res_init_oro}
            ncks -O -v OROMEA $rfile tmp_oromea
            ncrename -v OROMEA,oromea tmp_oromea
        fi
        change_lonlat $ifile oromea
        # LA: need to backup echam6_new_orography_before_calnoro!!!
        ncbo -O --operation=add $ifile tmp_oromea $ofile 2>> nco_stderr_ice2echam

        update_orography_apply_orography_anomaly
}

## @fn update_orography_apply_orography_anomaly() 
update_orography_apply_orography_anomaly() {
        echo; echo "                *   applying surface height changes to low-resolution echam topography"
        ifile=${COUPLE_DIR}/ice_orog_difference_${RES_echam}grid.nc
        ofile=${COUPLE_DIR}/echam6_new_orography_before_calnoro_low_res.nc
        ncrename -v ${ice_orography_difference_name},oromea $ifile
        cdo -s -setmisstoc,0 $ifile tmp; mv tmp $ifile
        if [ -L ${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc ]; then
                rfile=$(readlink ${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc)
        else
                rfile=${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc
        fi
        ncks -O -v oromea $rfile tmp_oromea
        change_lonlat $ifile oromea
        ncbo -O --operation=add $ifile tmp_oromea $ofile 2>> nco_stderr_ice2echam
}

## @fn update_orography_prepare_calnoro_high_res() 
update_orography_prepare_calnoro_high_res() {
        echo; echo "                *   preparing for external program CALNORO"
        ifile=${COUPLE_DIR}/echam6_new_orography_before_calnoro_high_res.nc
        ofile=topodata.nc
        cdo -s -P 28 \
                 -remapnn,r720x360 \
                 -setgrid,T511grid \
                 -chname,oromea,OROMEA \
                 $ifile \
                 $ofile 2>>cdo_stderr_ice2echam
}

## @fn update_orography_prepare_calnoro() 
update_orography_prepare_calnoro() {
        echo; echo "                *   preparing for external program CALNORO"
        ifile=${COUPLE_DIR}/echam6_new_orography_before_calnoro_low_res.nc
        ofile=topodata.nc
        cdo -s -P 28 \
                 -remapnn,r720x360 \
                 -setgrid,${RES_echam}grid \
                 -chname,oromea,OROMEA \
                 $ifile \
                 $ofile 2>>cdo_stderr_ice2echam
}

## @fn update_orography_run_calnoro() 
update_orography_run_calnoro() {
        echo; echo "                *   running CALNORO"
        ofile=${COUPLE_DIR}/calnoro_output.nc
        if [ ! -f calnoro ]; then
                if [ ! -f ${FUNCTION_PATH}/../utils/calnoro ]; then
                        if [ -f ${FUNCTION_PATH}/../utils/install_calnoro.sh ]; then
                                echo "                   install_calnoro.sh found. Running installation ..."
                                currentdir=$(pwd)
                                cd ${FUNCTION_PATH}/../utils/
                                ./install_calnoro.sh ${MACHINE}
                                cd ${currentdir}
                                cp ${FUNCTION_PATH}/../utils/calnoro .
                        else
                                echo "                   !!!   FATAL ERROR: Compiled external utility CALNORO not found"
                                echo "                   !!!   Looking in: ${FUNCTION_PATH}/../utils/"
                                echo "                   !!!   FATAL ERROR: Please compile this program using the install_calnoro.sh script!"
                                exit 42
                        fi
                else
                        cp ${FUNCTION_PATH}/../utils/calnoro .
                fi
        fi
        # Make sure dimension order is correct:
        change_lonlat topodata.nc OROMEA
        
        module list
        if [ "$MACHINE" == "levante" ]; then
            export LD_LIBRARY_PATH=/sw/spack-levante/netcdf-fortran-4.5.3-jlxcfz/lib/:$LD_LIBRARY_PATH
        elif [ "$MACHINE" == "albedo" ]; then
            export LD_LIBRARY_PATH=/albedo/soft/sw/spack-sw/netcdf-fortran/4.5.4-yb7woqz/lib/:$LD_LIBRARY_PATH
	        tmp_mods=$( module list -t | tail -n +2 )
            module purge
        fi

        echo ${RES_echam/T/} | ./calnoro > ${COUPLE_DIR}/calnoro_stdout_stderr 2>&1 
        
        if [ "$MACHINE" == "albedo" ]; then
            module load $tmp_mods
        fi
        
        cdo -s -f nc -setgrid,$RES_echam -invertlat sso_par_fil.srv $ofile 2>>cdo_stderr_ice2echam
        #rm sso_par_fil.srv topodata.nc
}

## @fn update_orography_generate_target_orography_file() 
update_orography_generate_target_orography_file() { 
        echo; echo "                *   generating target orography file for linear stepwise change of echam6 orography"
        ifile=${COUPLE_DIR}/calnoro_output.nc
        ofile=${COUPLE_DIR}/target_orography_echam6_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc

        typeset -A orography_names
        orography_names[var51]=OROMEA
        orography_names[var52]=OROSTD
        orography_names[var53]=OROSIG
        orography_names[var54]=OROGAM
        orography_names[var55]=OROTHE
        orography_names[var56]=OROPIC
        orography_names[var57]=OROVAL

        for varname in "${!orography_names[@]}"; do
                ncrename -v ${varname},${orography_names[$varname]} ${ifile} tmp; mv tmp ${ifile}
        done

        ## Ensure that we use the OROMEA from the "raw" PISM-corrected orography,
        ## not what is given from the calnoro program


        ############################################################
        # Define the method for orography update (oro_update_mod):
        #   1: use low-resolution reference; take all sub-grid scale parameters from CALNORO; take GEOSP from low-res input
        #   2: use low-resolution reference; take all sub-grid scale parameters from CALNORO; take OROMEA and GEOSP from low-res input
        #                                                                                     (like it has been done before)
        #   3: use high-resolution reference; take all sub-grid scale parameterds from CALNORO; take OROMEA and GEOSP from high-res input
        #   4: update orography but keep sub-grid scale parameters from unit.24
        #   5: update orography but change sub-grid scale parameters only in ice domain where ice mask or change of ice mask (ice loss/ice gain) 

        if [ "x${oro_update_mod}" == "x1" ]; then # according to Martins script: OROMEA from calnoro; GEOSP from input
                
                ECHAM_g=9.80665
                cdo -s -sp2gp -gp2sp -mulc,${ECHAM_g} -chname,oromea,GEOSP \
                        -remapnn,${RES_echam}grid \
                        ${COUPLE_DIR}/echam6_new_orography_before_calnoro_low_res.nc \
                        ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc 2>>cdo_stderr_ice2echam

                # Make sure dimension order is correct:
                change_lonlat ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc GEOSP

                cdo -s -O merge \
                        ${ifile} \
                        ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc \
                        ${ofile} 2>>cdo_stderr_ice2echam
        elif [ "x${oro_update_mod}" == "x2" ]; then # old coupling scheme: OROMEA and GEOSP from low-res input
                cdo -s -replace \
                        ${ifile} \
                        -chname,oromea,OROMEA -remapnn,${RES_echam}grid ${COUPLE_DIR}/echam6_new_orography_before_calnoro_low_res.nc \
                        tmp 2>>cdo_stderr_ice2echam
                mv tmp ${ifile}

                # Take the GEOSP field from the "raw" PISM-corrected orography, and not
                # what is given from the calnoro program 
                
                ECHAM_g=9.80665
                cdo -s -sp2gp -gp2sp -mulc,${ECHAM_g} -chname,oromea,GEOSP \
                        -remapnn,${RES_echam}grid \
                        ${COUPLE_DIR}/echam6_new_orography_before_calnoro_low_res.nc \
                        ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc 2>>cdo_stderr_ice2echam

                # Make sure dimension order is correct:
                change_lonlat ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc GEOSP

                cdo -s -O merge \
                        ${ifile} \
                        ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc \
                        ${ofile} 2>>cdo_stderr_ice2echam
        elif [ "x${oro_update_mod}" == "x3" ]; then # new coupling scheme: OROMEA and GEOSP from high-res input
                cdo -s -replace \
                        ${ifile} \
                        -chname,oromea,OROMEA -remapnn,${RES_echam}grid ${COUPLE_DIR}/echam6_new_orography_before_calnoro_high_res.nc \
                        tmp 2>>cdo_stderr_ice2echam
                mv tmp ${ifile}

                # Take the GEOSP field from the "raw" PISM-corrected orography, and not
                # what is given from the calnoro program 
                
                ECHAM_g=9.80665
                cdo -s -sp2gp -gp2sp -mulc,${ECHAM_g} -chname,oromea,GEOSP \
                        -remapnn,${RES_echam}grid \
                        ${COUPLE_DIR}/echam6_new_orography_before_calnoro_high_res.nc \
                        ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc 2>>cdo_stderr_ice2echam

                # Make sure dimension order is correct:
                change_lonlat ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc GEOSP

                cdo -s -O merge \
                        ${ifile} \
                        ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc \
                        ${ofile} 2>>cdo_stderr_ice2echam
                 #
        elif [ "x${oro_update_mod}" == "x4" ]; then # new coupling scheme: OROMEA and GEOSP from high-res input but overwrite with unit.24 sub-grid parameters
                cdo -s -replace \
                        ${ifile} \
                        -chname,oromea,OROMEA -remapnn,${RES_echam}grid ${COUPLE_DIR}/echam6_new_orography_before_calnoro_high_res.nc \
                        tmp 2>>cdo_stderr_ice2echam
                mv tmp ${ifile}

                # Take the GEOSP field from the "raw" PISM-corrected orography, and not
                # what is given from the calnoro program 
                
                ECHAM_g=9.80665
                cdo -s -sp2gp -gp2sp -mulc,${ECHAM_g} -chname,oromea,GEOSP \
                        -remapnn,${RES_echam}grid \
                        ${COUPLE_DIR}/echam6_new_orography_before_calnoro_high_res.nc \
                        ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc 2>>cdo_stderr_ice2echam

                # Make sure dimension order is correct:
                change_lonlat ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc GEOSP

                cdo -s -O merge \
                        ${ifile} \
                        ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc \
                        ${ofile} 2>>cdo_stderr_ice2echam
                
                cdo -s replace ${ofile} ${INIT_DIR_echam}/unit.24 tmp_replace; mv tmp_replace ${ofile}
                 #
        
        elif [ "x${oro_update_mod}" == "x5" ]; then # masked: new coupling scheme: OROMEA and GEOSP from high-res input
                cdo -s -replace \
                        ${ifile} \
                        -chname,oromea,OROMEA -remapnn,${RES_echam}grid ${COUPLE_DIR}/echam6_new_orography_before_calnoro_high_res.nc \
                        tmp 2>>cdo_stderr_ice2echam
                mv tmp ${ifile}

                # Take the GEOSP field from the "raw" PISM-corrected orography, and not
                # what is given from the calnoro program 
                
                ECHAM_g=9.80665
                cdo -s -sp2gp -gp2sp -mulc,${ECHAM_g} -chname,oromea,GEOSP \
                        -remapnn,${RES_echam}grid \
                        ${COUPLE_DIR}/echam6_new_orography_before_calnoro_high_res.nc \
                        ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc 2>>cdo_stderr_ice2echam

                # Make sure dimension order is correct:
                change_lonlat ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc GEOSP

                cdo -s -O merge \
                        ${ifile} \
                        ${COUPLE_DIR}/echam6_new_GEOSP_for_target_orog.nc \
                        ${ofile} 2>>cdo_stderr_ice2echam
                
                if [ -L ${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc ]; then
                        rfile=$(readlink ${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc)
                else
                        rfile=${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc
                fi
                
                cdo -s -ifthen ${COUPLE_DIR}/ice_domain_current_T63grid.nc -gtc,0 \
                        -add ${COUPLE_DIR}/ice_mask_current_T63grid.nc -selname,GLAC ${INIT_DIR_echam}/unit.24 \
                        ${COUPLE_DIR}/mask_oro_update.nc
#                cdo -s selname,glac ${rfile} ${COUPLE_DIR}/ice_mask_before_T63grid.nc
                for VAR in "OROMEA" "OROSTD" "OROSIG" "OROGAM" "OROTHE" "OROPIC" "OROVAL"
                do
                    cdo -s ifthenelse -setmisstoc,0 ${COUPLE_DIR}/mask_oro_update.nc \
                            -selname,${VAR} ${ofile} \
                            -selname,${VAR} ${INIT_DIR_echam}/unit.24 ${VAR}.orotmp

                done
                cdo -s merge *.orotmp ORO_merged #; mv ORO_merged ${ofile}; rm *.orotmp
                cdo -s replace ${ofile} ORO_merged tmp; mv tmp ${ofile}
                 #
        fi
        ############################################################
        ln -sf ${COUPLE_DIR}/target_orography_echam6_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc ${COUPLE_DIR}/latest_target_orography.nc
        # Destroy# the time dimension; just in case
	#if ncdump -h# $ofile | grep -q "time"; then
	#	ncwa -a time $ofile tmp; mv tmp $ofile
	#fi
        echo; echo -e "\t\t-   finished with target orography generation"
}

## @fn update_orography_set_namelist_modifications_next_echam_run() 
update_orography_set_namelist_modifications_next_echam_run() {
        echam_variable_modify_file=${COUPLE_DIR}/echam_namelist_switches_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.dat  
        echo; echo "                *   setting echam variables for update ororgapy to be read from file in next run"
              echo "                -   echam_variable_modify_file=$echam_variable_modify_file"
        if [ ! -s $echam_variable_modify_file ]; then
                :> $echam_variable_modify_file
        fi
        echo "submodelctl___lupdate_orog___nml_entry=.TRUE." >> $echam_variable_modify_file
        echo "submodelctl___lupdate_orog___nml_file=namelist.echam" >> $echam_variable_modify_file
        add_to ${echam_variable_modify_file} echam_namelist_switches.dat config echam
}

## @fn set_glacial_mask_select_field() 
set_glacial_mask_select_field() {
        echo; echo "                *   selecting glacial mask"
        cdo \
            -selvar,${ice_glacial_mask_name} \
            ${COUPLE_DIR}/ice_file_for_atmosphere.nc \
            ${COUPLE_DIR}/ice_mask_current.nc 2>> cdo_stderr_ice2echam
	mv ${COUPLE_DIR}/ice_mask_current.nc tmp
        ncdump tmp | sed 's/byte/float/g' | ncgen -o tmp.nc; mv tmp.nc tmp
        mv tmp ${COUPLE_DIR}/ice_mask_current.nc
        rmlist="$rmlist ${COUPLE_DIR}/ice_mask_current.nc"

}

## @fn set_glacial_mask_echam_restart_file() 
set_glacial_mask_echam_restart_file() { 
        if [ -L ${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc ]; then 
                ECHAM_RESTART_FILE=$(readlink ${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc)
        else
                ECHAM_RESTART_FILE=${RESTART_DIR_echam}/restart_${EXP_ID}_echam.nc
        fi
        ifile=${COUPLE_DIR}/ice_mask_current.nc
        dfile=${COUPLE_DIR}/ice_domain_current.nc
        cdo -s setrtoc,0,1,1 $ifile $dfile 

        echo; echo "                *   Setting ice mask in: $ECHAM_RESTART_FILE"
              echo "                                   from: $ifile"
              echo "                           using region: $dfile"

        # Define the entire area where pism is defined:
        ncrename -v ${ice_glacial_mask_name},ice_mask   $ifile
        ncrename -v ${ice_glacial_mask_name},ice_domain $dfile

        regrid_to_echam $ifile
        regrid_to_echam $dfile

        ifile=${ifile%.*}_${RES_echam}grid.nc
        dfile=${dfile%.*}_${RES_echam}grid.nc

        cdo -s -setmisstoc,0 $ifile tmp 2>> cdo_stderr_ice2echam; mv tmp $ifile 
        cdo -s -gec,${ECHAM_GLACIAL_THRESHOLD} $ifile tmp 2>> cdo_stderr_ice2echam; mv tmp $ifile   # PG: Ensure the mask is binary after regridding.

        change_lonlat $ifile ice_mask
        change_lonlat $dfile ice_domain

        # Restrict the ice mask on the echam grid to only be valid where echam land
        cdo -s -setmisstoc,0 -ifthen -selvar,slm $ECHAM_RESTART_FILE $ifile tmp 2>> cdo_stderr_ice2echam; mv tmp $ifile

        cp $ECHAM_RESTART_FILE ${ECHAM_RESTART_FILE%.*}_backup_before_glacial_mask_replacement.nc
        # FIX: Skip ECHAM restart modification to avoid malloc error
        # The original code using ncdump/ncgen/ncks corrupts the NetCDF structure
        # JSBACH restart updates are sufficient for vegetation recovery
        echo "                -   Skipping ECHAM restart glac modification (using backup as-is)"
        echo "                    JSBACH restart files will be updated separately"
        
        ncrename -v ice_mask,glac $ifile
}

## @fn set_glacial_mask_echam_make_dummy_jan_surf() 
set_glacial_mask_echam_make_dummy_jan_surf() {
        echo; echo "                *   Generating a dummy unit.24 file with new ice orography and mask"
        # Generates a dummy unit.24 file by modifying a copy of the
        # The echam restart file has already been modified with the "correct"
        # glacial mask", construct a dummy one from the restart
	#CLEANUP cp ${BASE_DIR}/${EXP_ID}/awicm/input/echam/unit.24 dummy_jan_surf.nc
        cp ${INIT_DIR_echam}/unit.24 dummy_jan_surf.nc
        # vvvvvvvvvvvvv  LARS的代码加在这里 vvvvvvvvvvvvv
        if [[ ${DOMAIN_pism} == "greenland" ]]; then
                cdo setclonlatbox,0,0,360,0,90 -selvar,GLAC dummy_jan_surf.nc tmp
                cdo -replace dummy_jan_surf.nc tmp dummy_jan_surf.nc2
                mv dummy_jan_surf.nc2 dummy_jan_surf.nc
                rm tmp
        fi
        # ^^^^^^^^^^^^  LARS的代码加在这里 ^^^^^^^^^^^^
       	# Fix Lu: ---------------------------
        if [[ ${DOMAIN_pism} == "nhem" ]]; then
                cdo setclonlatbox,0,0,360,0,90 -selvar,GLAC dummy_jan_surf.nc tmp
                cdo -replace dummy_jan_surf.nc tmp dummy_jan_surf.nc2
                mv dummy_jan_surf.nc2 dummy_jan_surf.nc
                rm tmp
        fi
        # Fix Lu: ---------------------------
        # FIX 2025-01-24: Check if target_orography file exists before using it
        # The target_orography file is generated in update_orography step which runs AFTER update_glacial_mask
        # On first coupling iteration, this file may not exist yet
        if [ -f ${COUPLE_DIR}/target_orography_echam6_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc ]; then
                cdo -s replace dummy_jan_surf.nc \
                        ${COUPLE_DIR}/target_orography_echam6_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc \
                        tmp && mv tmp dummy_jan_surf.nc
        elif [ -f ${COUPLE_DIR}/latest_target_orography.nc ]; then
                # Try using the latest target orography from previous iteration
                echo "                -   Using latest_target_orography.nc from previous iteration"
                cdo -s replace dummy_jan_surf.nc \
                        ${COUPLE_DIR}/latest_target_orography.nc \
                        tmp && mv tmp dummy_jan_surf.nc
        else
                echo "                -   WARNING: No target_orography file found, skipping orography update in dummy_jan_surf.nc"
                echo "                -   This is expected on the first coupling iteration"
        fi
        # Fix Albedo where there are glaciers:
        #echam_albedo_on_glaciers=0.7
        echam_albedo_on_glaciers=${ECHAM_ALBEDO_ON_GLACIERS}
        cdo -s \
                -setmisstoc,${echam_albedo_on_glaciers} -setrtomiss,-1,2 \
                ${COUPLE_DIR}/ice_mask_current_${RES_echam}grid.nc \
                input_if_mask_true.nc
        cdo -s \
                -selvar,ALB dummy_jan_surf.nc input_if_mask_false.nc
        cdo -s -ifthenelse \
                ${COUPLE_DIR}/ice_mask_current_${RES_echam}grid.nc \
                input_if_mask_true.nc \
                input_if_mask_false.nc \
                dummy_jan_surf_ALB.nc
        cdo -s replace \
               dummy_jan_surf.nc \
               -chname,${echam_glacial_mask_name},ALB dummy_jan_surf_ALB.nc \
               tmp; mv tmp dummy_jan_surf.nc

        # Fix lakes where there are glaciers
        cdo -s \
                -setrtoc,0,1,0 -selvar,ALAKE dummy_jan_surf.nc input_if_mask_true.nc
        cdo -s \
                -selvar,ALAKE dummy_jan_surf.nc input_if_mask_false.nc
        cdo -s -ifthenelse \
                ${COUPLE_DIR}/ice_mask_current_${RES_echam}grid.nc \
                input_if_mask_true.nc \
                input_if_mask_false.nc \
                dummy_jan_surf_ALAKE.nc
        cdo -s replace \
                dummy_jan_surf.nc \
                dummy_jan_surf_ALAKE.nc tmp #-chname,${echam_glacial_mask_name},ALAKE dummy_jan_surf_ALAKE.nc \
        mv tmp dummy_jan_surf.nc

        # === BUG FIX: 退冰区域GLAC应该从1变为0 ===
        # 原代码bug: 退冰时保持原来的GLAC=1，导致desert_fpc=1，植被无法生长
        # 修复: 直接使用PISM的ice_mask作为GLAC
        #       ice_mask=1 -> GLAC=1 (有冰川)
        #       ice_mask=0 -> GLAC=0 (无冰川，允许植被生长)
        cdo -s \
                -setmisstoc,1 -setrtomiss,-1,2 \
                ${COUPLE_DIR}/ice_mask_current_${RES_echam}grid.nc \
                input_if_mask_true.nc
        # 退冰区域应该设为0，而不是保持原来的GLAC值
        cdo -s \
                -setmisstoc,0 -setrtomiss,-1,2 \
                ${COUPLE_DIR}/ice_mask_current_${RES_echam}grid.nc \
                input_if_mask_false.nc

        cdo -s -ifthenelse \
                ${COUPLE_DIR}/ice_mask_current_${RES_echam}grid.nc \
                input_if_mask_true.nc \
                input_if_mask_false.nc \
                dummy_jan_surf_GLAC.nc

        cdo -s replace \
                dummy_jan_surf.nc \
                -chname,${echam_glacial_mask_name},GLAC dummy_jan_surf_GLAC.nc \
                tmp; mv tmp dummy_jan_surf.nc

        # === AUTO-UPDATE unit.24 GLAC mask (input + run/input + work) ===
        # Why this is needed:
        # - ECHAM reads ${WORK_DIR}/unit.24 at startup
        # - esm_runscripts may restore ${WORK_DIR}/unit.24 from run/input/echam/unit.24
        # - Therefore BOTH must have the same (updated) GLAC mask as used in jsbach.nc (and PISM ice mask)
        #
        # IMPORTANT:
        # We update ONLY the GLAC variable in unit.24 via "cdo replace" to avoid relying on timestamps/metadata
        # and to keep other unit.24 fields untouched.
        echo "                -   Updating GLAC in unit.24 files (input + run/input + work)"
        cdo -s -chname,${echam_glacial_mask_name},GLAC dummy_jan_surf_GLAC.nc unit24_glac_update.nc

        # 1) Experiment-level input (used as template in some setups)
        if [ -f "${INIT_DIR_echam}/unit.24" ]; then
                cdo -s replace ${INIT_DIR_echam}/unit.24 unit24_glac_update.nc tmp_unit24 && mv tmp_unit24 ${INIT_DIR_echam}/unit.24
                echo "                    --> Updated: ${INIT_DIR_echam}/unit.24"
        else
                echo "                    WARNING: ${INIT_DIR_echam}/unit.24 not found"
        fi

        # 2) Current run work directory (ECHAM reads this)
        if [ -n "${WORK_DIR}" ] && [ -f "${WORK_DIR}/unit.24" ]; then
                cdo -s replace ${WORK_DIR}/unit.24 unit24_glac_update.nc tmp_unit24 && mv tmp_unit24 ${WORK_DIR}/unit.24
                echo "                    --> Updated: ${WORK_DIR}/unit.24"
        else
                echo "                    WARNING: WORK_DIR not set or ${WORK_DIR}/unit.24 not found"
        fi

        # 3) Current run input copy (may be used to restore protected files into work/)
        if [ -n "${WORK_DIR}" ] && [ -f "${WORK_DIR}/../input/echam/unit.24" ]; then
                cdo -s replace ${WORK_DIR}/../input/echam/unit.24 unit24_glac_update.nc tmp_unit24 && mv tmp_unit24 ${WORK_DIR}/../input/echam/unit.24
                echo "                    --> Updated: ${WORK_DIR}/../input/echam/unit.24"
        else
                echo "                    WARNING: run/input/echam/unit.24 not found (WORK_DIR/../input/echam/unit.24)"
        fi
        # === END AUTO-UPDATE ===
}

## @fn set_glacial_mask_jsbach_update_init_file() 
set_glacial_mask_jsbach_update_init_file() {
        echo; echo "                *   Generating a new JSBACH init file with updated land cover types"
        # Get the Program and required .mod files
        if [ -f ${FUNCTION_PATH}/../utils/jsbach_init_file ]; then
                cp ${FUNCTION_PATH}/../utils/jsbach_init_file .
        elif [ -f ${FUNCTION_PATH}/../utils/install_jsbach_init_file.sh ]; then
                echo "                   install_jsbach_init_file.sh found. Running installation ..."
                currentdir=$(pwd)
                cd ${FUNCTION_PATH}/../utils/
                ./install_jsbach_init_file.sh ${MACHINE}
                cd ${currentdir}
                cp -v ${FUNCTION_PATH}/../utils/jsbach_init_file .
        else
                echo -e "\t\t E R O R R --- binary missing, don't know how to proceed"
                exit
        fi
        cp ${FUNCTION_PATH}/../utils/_build/mo_kinds.mod .
        cp ${FUNCTION_PATH}/../utils/_build/mo_vegparams.mod .
        cleanup_list="$cleanup_list jsbach_init_file mo_kinds.mod mo_vegparams.mod"
        # Here we place a slimmed-down version of the Veronika Gayler script:
        res_atm=${RES_echam}
        res_oce=GENERIC
        ntiles=11                   # number of jsbach tiles

        dynveg=true                 # setup for dynamic vegetation 
        c3c4crop=true               # differentiate between C3 and C4 crops

        lpasture=true               # distinguish pastures from grasses
        read_pasture=LUH2v2h        # LUH: read pastures and crops from LUH states as in CMIP5
                                    # LUH2v2h: read pastures and crops from LUH2 states as in CMIP6
                                    # false: no separate input file for crops and pastures
        pasture_rule=true           # allocate pastures primarily on grass lands

        year_ct=1850                # year the cover_types are derived from
        year_cf=1850                # year cover fractions are derived from

        landcover_series=false      # generate a series of files with cover_types of 
                                    # year_ct and fractions from year_cf to year_cf2
        year_cf2=1859               # only used with landcover_series

        echam_fractional=false      # initial file for echam runs with fractional
                                    # land sea mask
        masks_file=default          # file with land sea mask (default: use echam land sea mask)

        mv dummy_jan_surf.nc ${RES_echam}${res_oce}_jan_surf.nc
        echo "                -   Using unit.24 (jan surf) file: ${RES_echam}GENERIC_jan_surf.nc"
        echo "                    --> (This is the dummy file after ice update)"
        ln -sf ${INIT_DIR_echam}/unit.91 ${res_atm}${res_oce}_VGRATCLIM.nc
        echo "                -   Using unit.91 (VGRATCLIM) file: ${INIT_DIR_echam}/unit.91"
        ln -sf ${INIT_DIR_echam}/unit.90 ${res_atm}${res_oce}_VLTCLIM.nc
        echo "                -   Using unit.90 (VLTCLIM) file: ${INIT_DIR_echam}/unit.90"

        cleanup_list="$cleanup_list ${RES_echam}_jan_surf.nc"

        pool=${POOL_DIR_echam}/ECHAM6/${RES_echam}
        pool_land=${POOL_DIR_jsbach}/JSBACH/prepare/${RES_echam}
        #------------------------------------------------------------------------------
        # prepare the namelist
        #------------------------------------------------------------------------------

        [[ ${res_oce} = "" ]] && lcouple=.false. || lcouple=.true.
        [[ ${dynveg} = true ]] && ldynveg=.true. || ldynveg=.false.
        [[ ${c3c4crop} = true ]] && lc3c4crop=.true. || lc3c4crop=.false.
        [[ ${read_pasture} != false ]] && lread_pasture=.true. || lread_pasture=.false.
        [[ ${landcover_series} = false ]] && year_cf2=${year_cf}

        desert_only=.false.         # setup for a desert-only experiment
        grass_only=.false.          # setup for a grass-only experiment
        woods_only=.false.          # setup for a woods-only experiment

        cat > namelist <<EOF
&INITCTL
  res_atm="${res_atm}"
  res_oce="${res_oce}"
  ntiles=${ntiles}
  nlct=21
  year_ct=${year_ct}
  year_cf=${year_cf}
  lcouple=${lcouple}
  ldynveg=${ldynveg}
  lc3c4crop=${lc3c4crop}
  lpasture=${lpasture}
  lread_pasture=${lread_pasture}
  pasture_tag="${read_pasture}"
  lpasture_rule=.${pasture_rule}.
  echam_fractional=.${echam_fractional}.
  masks_file="${masks_file##*/}"
  desert_only=${desert_only}
  grass_only=${grass_only}
  woods_only=${woods_only}
  cover_fract_only=.${landcover_series}.
  info=.true.
/
EOF
        echo -e "\t\t - This namelist was passed to the jsbach_init program:"
        cat namelist
        ln -sf ${pool}/${res_atm}_TSLCLIM2.nc   .

        if [ $(echo ${res_atm} | cut -c1) != T ]; then
                res_atmg=${res_atm}
        else
                res_atmg=${res_atm}gauss
        fi
        if [[ ${read_pasture} = false ]]; then
          ln -sf ${pool_land}/vegtype_${year_cf}_${res_atmg}_pa14.nc\
                              vegtype_${year_cf}_${res_atm}gauss_pa14.nc
          if [[ ${year_cf} != ${year_ct} ]]; then
            ln -sf ${pool_land}/vegtype_${year_ct}_${res_atmg}_pa14.nc \
                              vegtype_${year_ct}_${res_atm}gauss_pa14.nc
          fi
        fi
        ln -sf ${pool_land}/vegmax_6_${res_atm}_0-360.nc  \
                            vegmax_6_${res_atm}.nc
        ln -sf ${pool_land}/${res_atm}_topo_75.nc     .
        ln -sf ${pool_land}/albedo_${res_atm}.nc      .
        ln -sf ${pool_land}/C3C4_mask_${res_atmg}.nc \
                            C3C4_mask_${res_atm}gauss.nc
        ln -sf ${pool_land}/potveg_${res_atm}.nc      .
        if [[ ${c3c4crop} = true ]]; then
          if [[ ${read_pasture} = LUH2v2h ]]; then
            ln -sf ${pool_land}/C3C4_crop_LUH2v2h_${res_atm}.nc \
                                C3C4_crop_${res_atm}.nc
          else
            ln -sf ${pool_land}/C3C4_crop_${res_atm}.nc \
                                C3C4_crop_${res_atm}.nc
          fi
        fi
        if [[ ${read_pasture} = LUH2v2h ]]; then
          if [[ ${dynveg} = true ]]; then
            LUH_states=LUH2v2h_states_${res_atm}_dynveg.nc
          else
            LUH_states=LUH2v2h_states_${res_atm}_all-oceans_no-dynveg.nc
          fi
        elif [[ ${read_pasture} = LUH ]]; then
          LUH_states=LUH_states_${res_atm}.nc
        fi
        if [[ ${read_pasture} != false ]]; then
          if [[ -f ${pool_land}/${LUH_states} ]]; then
            ln -sf ${pool_land}/${LUH_states} .
          elif  [[ -f ${pool_land}/${LUH_states}.gz ]]; then
            cp ${pool_land}/${LUH_states}.gz .
            gunzip ${LUH_states}.gz
          fi
          #typeset -Z4 yr=${year_cf}
          yr=${year_cf}
          while [[ ${yr} -le ${year_cf2} ]]; do
            cdo -s selyear,${yr} ${LUH_states} LUH_states_${yr}_${res_atm}.nc
            (( yr = yr + 1 ))
          done
        fi
        [[ ${masks_file} != default && ! -f ${masks_file##*/} ]] && ln -sf ${masks_file} .

        ln -sf ${pool_land}/soil_parameters_${res_atm}.nc .
        echo "                -   Running jsbach_init_file (working dir: `pwd`)"
        chmod 755 jsbach_init_file
        ./jsbach_init_file > ${COUPLE_DIR}/jsbach_init_file_stdout_stderr 2>&1
        if [[ $? -eq 0 ]]; then
                echo -e "\t\t -   ...finished with external program jsbach_init_file!"
        else
                echo "error in jsbach_init_file"
                exit 1
        fi
        # BUG: How do we know what the output file is called?
        mv jsbach_T63GENERIC_11tiles_5layers_1850_dynveg.nc \
                ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc  
        
        # === Fix tiny values (1e-10) in cover_fract ===
        # The jsbach_init_file program generates 1e-10 values for tiles that should be 0.
        # These tiny values cause JSBACH's fpc_to_cover_fract_pot validation to fail.
        echo -e "\t\t -   Fixing tiny cover_fract values in jsbach.nc..."
        if [ -f ${FUNCTION_PATH}/../utils/fix_jsbach_tiny_values.py ]; then
            /sw/spack-levante/mambaforge-22.9.0-2-Linux-x86_64-kptncg/bin/python3 \
                ${FUNCTION_PATH}/../utils/fix_jsbach_tiny_values.py \
                ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc \
                >> ${COUPLE_DIR}/jsbach_init_file_stdout_stderr 2>&1
            if [ $? -eq 0 ]; then
                echo -e "\t\t -   ...tiny values fixed successfully!"
            else
                echo -e "\t\t -   WARNING: fix_jsbach_tiny_values.py failed, continuing anyway"
            fi
        else
            echo -e "\t\t -   WARNING: fix_jsbach_tiny_values.py not found, skipping"
        fi
        # === End fix tiny values ===
        
        echo :> ${COUPLE_DIR}/jsbach_init_override_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.dat
        
        ############################################################################
        # LA UPDATE
        ln -sf ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc ${COUPLE_DIR}/latest_jsbach_init_file.nc
        #cp ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc ${WORK_DIR}/jsbach.nc
        ############################################################################
        
        echo "LAND_BOUNDARY_CONDITIONS_jsbach=${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc" > ${COUPLE_DIR}/jsbach_init_override_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.dat
        add_to ${COUPLE_DIR}/jsbach_init_override_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.dat jsbach_init_override.dat config jsbach
}

set_glacial_mask_jsbach_update_restart_jsbach() {
        echo; echo "                *   updating land cover fractions in jsbach restart file >>jsbach<<"
        
        # === DISABLED: pack/unpack process corrupts the restart file ===
        # The unpack -> modify -> pack workflow introduces floating-point precision
        # errors that cause JSBACH to fail with "sum of cover_fract_pot /= 1" error.
        # The jsbach.nc file (cover_fract) is already updated by jsbach_init_file.
        # Dynamic vegetation will handle the actual vegetation dynamics.
        echo -e "\t\t- SKIPPED: jsbach restart modification disabled (pack/unpack corrupts file)"
        newest_jsbach_restart_file=$(readlink ${RESTART_DIR_jsbach}/restart_${EXP_ID}_jsbach.nc)
        ln -sf ${newest_jsbach_restart_file} ${COUPLE_DIR}/latest_jsbach_restart_file.nc
        
        # === Fix tiny values in jsbach restart file ===
        # The jsbach restart may contain 1e-10 values from previous model runs
        echo -e "\t\t- Fixing tiny values in jsbach restart file..."
        if [ -f ${FUNCTION_PATH}/../utils/fix_jsbach_tiny_values.py ]; then
            /sw/spack-levante/mambaforge-22.9.0-2-Linux-x86_64-kptncg/bin/python3 \
                ${FUNCTION_PATH}/../utils/fix_jsbach_tiny_values.py \
                ${newest_jsbach_restart_file} --restart \
                >> ${COUPLE_DIR}/jsbach_restart_fix_stdout_stderr 2>&1
            if [ $? -eq 0 ]; then
                echo -e "\t\t- ...jsbach restart tiny values fixed!"
            else
                echo -e "\t\t- WARNING: jsbach restart fix failed, continuing anyway"
            fi
        fi
        # === End fix ===
        
        # === FIX 2025-01-24: Update cover_fract in jsbach restart to match jsbach.nc ===
        # Without this fix, cover_fract in restart shows glacier (Tile 1=1) while
        # jsbach.nc shows non-glacier (e.g., Tile 7=1), causing JSBACH to crash with
        # "sum of cover_fract_pot /= 1" error in fpc_to_cover_fract_pot function.
        
        # FIX 2025-01-25: Create reference file HERE if it doesn't exist yet
        # This file is needed to identify truly deglaciated areas
        JSBACH_REFERENCE_FILE=${COUPLE_DIR}/jsbach_original_reference.nc
        if [ ! -f ${JSBACH_REFERENCE_FILE} ]; then
            echo -e "\t\t- Creating reference file: ${JSBACH_REFERENCE_FILE}"
            cp ${FORCING_DIR_jsbach}/jsbach.nc ${JSBACH_REFERENCE_FILE}
        fi
        
        echo -e "\t\t- Updating cover_fract in jsbach restart to match jsbach.nc..."
        if [ -f ${FUNCTION_PATH}/../utils/fix_jsbach_restart_cover_fract.py ]; then
            # IMPORTANT: Use jsbach_original_reference.nc to identify truly deglaciated areas
            # This avoids incorrectly modifying non-glacier Tile1 points (like tropical evergreen)
            /sw/spack-levante/mambaforge-22.9.0-2-Linux-x86_64-kptncg/bin/python3 \
                ${FUNCTION_PATH}/../utils/fix_jsbach_restart_cover_fract.py \
                ${newest_jsbach_restart_file} \
                ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc \
                ${JSBACH_REFERENCE_FILE} \
                >> ${COUPLE_DIR}/jsbach_restart_cover_fract_fix_stdout_stderr 2>&1
            if [ $? -eq 0 ]; then
                echo -e "\t\t- ...cover_fract in jsbach restart updated successfully!"
            else
                echo -e "\t\t- WARNING: cover_fract fix failed, continuing anyway"
            fi
        else
            echo -e "\t\t- WARNING: fix_jsbach_restart_cover_fract.py not found"
        fi
        # === End cover_fract fix ===
        
        # === FIX 2025-01-24: Sync jsbach restart to WORK_DIR ===
        # esm_tools copies restart files to WORK_DIR before coupling scripts run,
        # so our modifications in restart directory are not reflected in WORK_DIR.
        if [ -n "${WORK_DIR}" ] && [ -d "${WORK_DIR}" ]; then
            echo -e "\t\t- Syncing fixed jsbach restart to WORK_DIR..."
            if [ -f "${newest_jsbach_restart_file}" ]; then
                cp -v "${newest_jsbach_restart_file}" "${WORK_DIR}/restart_${EXP_ID}_jsbach.nc"
                echo -e "\t\t- jsbach restart synced to WORK_DIR successfully"
            fi
        fi
        # === End WORK_DIR sync ===
        
        return
        # === END DISABLED ===
        
        echo -e "\t\t-   Making a temporary directory to work in"
        startdir=$(pwd)
        if [ -d ${COUPLE_DIR}/update_land_cover_fractions ]; then
                # Throw away the temporary working directory from last time to
                # avoid merge problems in the next coupling iteration:
                rm -rf ${COUPLE_DIR}/update_land_cover_fractions
        fi
        mkdir -p ${COUPLE_DIR}/update_land_cover_fractions; cd ${COUPLE_DIR}/update_land_cover_fractions
        echo -e "\t\t-   determining where the glacial mask has advanced and retreated"

        # === Fix for iterative coupling: use fixed reference file instead of updating jsbach.nc ===
        JSBACH_REFERENCE_FILE=${COUPLE_DIR}/jsbach_original_reference.nc
        if [ ! -f ${JSBACH_REFERENCE_FILE} ]; then
            echo -e "\t\t-   Creating reference file for future comparisons: ${JSBACH_REFERENCE_FILE}"
            cp ${FORCING_DIR_jsbach}/jsbach.nc ${JSBACH_REFERENCE_FILE}
        fi
        # === End of fix ===

        echo -e "\t\t-   New JSBACH Init File: ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc"
        echo -e "\t\t-   Reference JSBACH File: ${JSBACH_REFERENCE_FILE}"
        cdo -s -sub \
            -selvar,glac ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc \
            -selvar,glac ${JSBACH_REFERENCE_FILE} \
            ${COUPLE_DIR}/update_land_cover_fractions/jsbach_glac_diff.nc 2>> ${COUPLE_DIR}/update_land_cover_fractions/cdo_stderr_ice2echam
        
        # Make two masks: one for glacial retreat, and one for glacial advance
        cdo -s gtc,0 ${COUPLE_DIR}/update_land_cover_fractions/jsbach_glac_diff.nc ${COUPLE_DIR}/update_land_cover_fractions/new_glaciers.nc
        cdo -s gtc,0 -mulc,-1 ${COUPLE_DIR}/update_land_cover_fractions/jsbach_glac_diff.nc ${COUPLE_DIR}/update_land_cover_fractions/new_non_glaciers.nc
        sum_new_glaciers=$(cdo -s -output -fldsum -gtc,0 ${COUPLE_DIR}/update_land_cover_fractions/jsbach_glac_diff.nc | tr -d "[:blank:]")
        sum_new_land=$(cdo -s -output -fldsum -gtc,0 -mulc,-1 ${COUPLE_DIR}/update_land_cover_fractions/jsbach_glac_diff.nc | tr -d "[:blank:]")
        echo -e "\t\t- Total number of new glacial cells: >>>$sum_new_glaciers<<<"
        echo -e "\t\t- Total number of new deglacial cells: >>>$sum_new_land<<<"

        if [[ $sum_new_glaciers -eq 0 ]] && [[ $sum_new_land -eq 0 ]]; then
                echo -e "\t\t- >>>INFO<<< No new land or new glaciers, this routine will be skipped..."
                newest_jsbach_restart_file=$(readlink ${RESTART_DIR_jsbach}/restart_${EXP_ID}_jsbach.nc)
                ln -sf ${newest_jsbach_restart_file} ${COUPLE_DIR}/latest_jsbach_restart_file.nc
                return
        fi

        newest_jsbach_restart_file=$(readlink ${RESTART_DIR_jsbach}/restart_${EXP_ID}_jsbach.nc)
        # Make a backup copy
        #cp $(readlink ${RESTART_DIR_jsbach}/restart_${EXP_ID}_jsbach.nc) ${RESTART_DIR_jsbach}/restart_${EXP_ID}_jsbach_backup_$(date +%Y%m%d%H%M%S).nc
        cp ${RESTART_DIR_jsbach}/restart_${EXP_ID}_jsbach.nc ${RESTART_DIR_jsbach}/restart_${EXP_ID}_jsbach_backup_$(date +%Y%m%d%H%M%S).nc
        echo -e "\t\t- Unpacking JSBACH restart file: $newest_jsbach_restart_file"
        # Ensure that the unpack binary is here:
        if [ ! -f unpack ]; then
                if [ -f ${FUNCTION_PATH}/../utils/unpack ]; then
                        #cp -v ${FUNCTION_PATH}/../utils/unpack ${COUPLE_DIR}/
                        cp -v ${FUNCTION_PATH}/../utils/unpack .
                elif [ -f ${FUNCTION_PATH}/../utils/install_jsbach_pack_unpack.sh ]; then
                        echo "                   install_jsbach_pack_unpack.sh found. Running installation ..."
                        currentdir=$(pwd)
                        cd ${FUNCTION_PATH}/../utils/
                        ./install_jsbach_pack_unpack.sh ${MACHINE}
                        cd ${currentdir}
                        #cp -v ${FUNCTION_PATH}/../utils/unpack ${COUPLE_DIR}/
                        cp -v ${FUNCTION_PATH}/../utils/unpack .
                else
                        echo -e "\t\t E R O R R --- unpack binary missing, don't know how to proceed"
                        exit
                fi
        fi
        
        echo " *** unpack binary found here: $(pwd) " 2>&1 > ${COUPLE_DIR}/unpack_script_stdout_stderr
        
        # Copy the files here:
        cp ${newest_jsbach_restart_file} ./
        cp ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc ./
        echo -e "\t\t- Running unpack binary with arguments:"
        echo -e "\t\t- Script: ${FUNCTION_PATH}/../utils/unpack_file.ksh"
        echo -e "\t\t- Restart File: $(basename ${newest_jsbach_restart_file})"
        echo -e "\t\t- Mask file: $(basename ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc)"
        echo "$(module list)"
        ${FUNCTION_PATH}/../utils/unpack_file.ksh \
                $(basename $newest_jsbach_restart_file) \
                $(basename ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc) \
                2>&1 > ${COUPLE_DIR}/unpack_script_stdout_stderr

        echo -e "\t\t- Listing files after program execution:"
        ls -ratl
        restartfile=2d_$(basename ${newest_jsbach_restart_file})
        # 'cdo showvar' does not work for variables with more than three dimensiones
        echo "ncdump -h ${restartfile} | grep 'double ' | cut -f 2 -d ' ' | cut -f1 -d '('"
        variables=$(ncdump -h ${restartfile} | grep 'double ' | cut -f 2 -d ' ' | cut -f1 -d '(')
        echo $variables

        # Rename the "tiles" dimension to "time", since cdo can only handle 4D
        # variables if one is called time. This will be corrected again below
        ncrename -O -d tiles,time 2d_$(basename ${newest_jsbach_restart_file}) 
        for var in ${variables}; do
                ncks -O --no_abc \
                        -v ${var} \
                        2d_$(basename ${newest_jsbach_restart_file}) \
                        2d_restart_${EXP_ID}_jsbach_${var}.nc

                if [ ! -f 2d_restart_${EXP_ID}_jsbach_${var}.nc ]; then
                        echo -e "\t\t !!! WARNING !!! -- no file for ${var}"
                        echo -e "\t\t                    continue loop"
                        continue
                fi

                # Keep the directory clean: make one directory for each variable:
                here=$(pwd); mkdir -p ${var}; mv 2d_restart_${EXP_ID}_jsbach_${var}.nc ${var}; cd ${var}
                echo -e "\t\t- Updating JSBACH restart file -- ${var}\n"
                case ${var} in
                        cover_fract | cover_fract_pot | veg_ratio)
                                # === FIX: Skip modification of these variables ===
                                # The pack/unpack process introduces floating-point precision errors
                                # that cause JSBACH to fail with "sum of cover_fract_pot /= 1" error.
                                # Dynamic vegetation will handle vegetation recovery in deglaciated areas.
                                echo -e "\t\t- Skipping ${var} modification (pack precision issues, dynveg will handle)"
                                cd $here
                                continue
                                ;;
                        cover_fract_DISABLED | cover_fract_pot_DISABLED | veg_ratio_DISABLED)
                                # Original code kept here for reference:
                                echo -e "\t\t- Need to set correct cover fractions for retreating and advancing glaciers for $var"
                                # Split to all the land tiles
                                #
                                # NOTE: We use splitrec because we need to
                                # split fake "timesteps" and not "layers" since
                                # cdo can't handle 4D without time...remember
                                # that tiles was renamed to time above.
                                cdo -s splitrec 2d_restart_${EXP_ID}_jsbach_${var}.nc 2d_restart_${EXP_ID}_jsbach_${var}_tile
                                mergelist=""
                                for tile in $(ls 2d_restart_${EXP_ID}_jsbach_${var}_tile0000??.nc); do
                                        case ${tile} in
                                                2d_restart_${EXP_ID}_jsbach_${var}_tile000001.nc )
                                                        echo -e "\t\t Glacier landpoints for ${var} on ${tile}"
                                                        # New glaciers must have exactly 1 or keep the old value
                                                        cdo -s ifthenelse \
                                                                ${COUPLE_DIR}/update_land_cover_fractions/new_glaciers.nc \
                                                                ${COUPLE_DIR}/update_land_cover_fractions/new_glaciers.nc \
                                                                ${tile} \
                                                                ${tile}.tmp
                                                        # New land must get 1e-10 or keep the old value
                                                        cdo -s ifthenelse \
                                                                ${COUPLE_DIR}/update_land_cover_fractions/new_non_glaciers.nc \
                                                                -mulc,1e-10 ${COUPLE_DIR}/update_land_cover_fractions/new_non_glaciers.nc \
                                                                ${tile}.tmp \
                                                                ${tile}

                                                        ;;
                                                2d_restart_${EXP_ID}_jsbach_${var}_tile000007.nc )
                                                        echo -e "\t\t Tundra for ${var} on ${tile}"
                                                        # New Glaciers need to be exactly 0 everywhere except on tile 1
                                                        cdo -s ifthenelse \
                                                                ${COUPLE_DIR}/update_land_cover_fractions/new_glaciers.nc \
                                                                -mulc,0 ${COUPLE_DIR}/update_land_cover_fractions/new_glaciers.nc \
                                                                ${tile} \
                                                                ${tile}.tmp
                                                        # New land gets the maximum allowed value
                                                        # NOTE: This is only valid for 11 tiles!!
                                                        # BUG(?): What about the negative values that appear here?
                                                        cdo -s -subc,1e-9 \
                                                                ${COUPLE_DIR}/update_land_cover_fractions/new_non_glaciers.nc \
                                                                ${tile}_maxveg.nc
                                                        cdo -s ifthenelse \
                                                                ${COUPLE_DIR}/update_land_cover_fractions/new_non_glaciers.nc \
                                                                ${tile}_maxveg.nc \
                                                                ${tile}.tmp \
                                                                ${tile}
                                                        ;;
                                                *)
                                                        echo -e "\t\t Nothing special for ${var} on ${tile}"
                                                        # New glacial tiles get exactly 0
                                                        cdo -s ifthenelse \
                                                                ${COUPLE_DIR}/update_land_cover_fractions/new_glaciers.nc \
                                                                -mulc,0 ${COUPLE_DIR}/update_land_cover_fractions/new_glaciers.nc \
                                                                ${tile} ${tile}.tmp
                                                        # New land cells get 1e-10: smallest allowed minium
                                                        cdo -s -ifthenelse \
                                                                ${COUPLE_DIR}/update_land_cover_fractions/new_non_glaciers.nc \
                                                                -mulc,1e-10 ${COUPLE_DIR}/update_land_cover_fractions/new_non_glaciers.nc \
                                                                ${tile}.tmp \
                                                                ${tile}
                                                        ;;
                                        esac
                                        # Add the fake time axis after the case
                                        # statement, since the ifthenelse
                                        # clause might do strange things with
                                        # the time axis otherwise:
                                        tile_number=$(echo $tile | rev | cut -d_ -f1 | rev | sed s/tile// | sed s/\.nc//)
                                        cdo -settaxis,${tile_number}-01-01,00:00:00 $tile ${tile}.time; mv ${tile}.time ${tile}
                                        mergelist="$mergelist $tile"
                                done
                                cdo -s -O mergetime $mergelist ${var}.tmp
                                # Ensure the name is correct
                                cdo -s setvar,$var ${var}.tmp 2d_restart_${EXP_ID}_jsbach_${var}.nc
                                ;;
                        soil_moisture)
                                echo -e "\t\t- Need to change soil moisture for $var"
                                ;;
                        layer_moisture)
                                echo -e "\t\t- Need to change layer moisture for $var"
                                # NOTE: The layer moisture can't exceed the
                                # holding capacity, which changes under (new)
                                # glacier cells. Since the glacial mask was
                                # used to construct the jsbach.nc forcing file,
                                # we use those values for new glacial cells:
                                #
                                # NOTE 2: the layer moisture variable is called
                                # layer_moist in jsbach init file, not layer
                                # moisture. Therefore, the setvar command needs
                                # to rename it
                                cdo -s -setvar,${var} -ifthenelse \
                                        ${COUPLE_DIR}/update_land_cover_fractions/new_non_glaciers.nc \
                                        -selvar,layer_moist ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc \
                                        2d_restart_${EXP_ID}_jsbach_${var}.nc 2d_restart_${EXP_ID}_jsbach_${var}.tmp
                                cdo -s -setvar,${var} -ifthenelse \
                                        ${COUPLE_DIR}/update_land_cover_fractions/new_glaciers.nc \
                                        -selvar,layer_moist ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc \
                                       2d_restart_${EXP_ID}_jsbach_${var}.tmp 2d_restart_${EXP_ID}_jsbach_${var}.nc
                                ;;
                esac
                cd $here
        done
        mergelist=""
        for var in ${variables}; do
                ncrename -d .time,tiles ${var}/2d_restart_${EXP_ID}_jsbach_${var}.nc 2>&1 >> ${COUPLE_DIR}/nco_stderr_ice2echam
                mergelist="$mergelist ${var}/2d_restart_${EXP_ID}_jsbach_${var}.nc"
        done
        echo -e "\t\t- Creating merge list"
        # For debugging:
        #echo "Done constructing merge list and renaming time to tiles, these files will be merged"
        #for f in $mergelist; do
        #        ls -l $f
        #done
        # Merge with ncks, since cdo has issues:
        echo -e "\t\t- Merging with ncks, as cdo has problems"
        ofile=2d_restart_${EXP_ID}_jsbach_merged.nc
        for f in $mergelist; do
                ncks -A -h $f $ofile
        done
        echo -e "\t\t- Done merging, now we need to fix the tiles dimension to be limited again rather than unlimited"
        # Turn "tiles" back to a limited dimension:
        ncks --fix_rec_dmn tiles $ofile -o tmp; mv tmp $ofile
        echo -e "\t\t- Finished with modification of variables, now we need to repack"
        
        # Ensure that the pack binary is here:
        if [ ! -f pack ]; then
                #if [ -f ${FUNCTION_PATH}..//utils/_build/pack ]; then
                #        cp -v ${FUNCTION_PATH}..//utils/_build/pack .
                if [ -f ${FUNCTION_PATH}/../utils/pack ]; then
                        cp -v ${FUNCTION_PATH}/../utils/pack .
                else
                        echo -e "\t\t E R O R R --- pack binary missing, don't know how to proceed"
                        exit
                fi
        fi
        # Copy the files here:
        cp ${newest_jsbach_restart_file} .
        cp ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc .
        echo -e "\t\t- Running pack binary with arguments:"
        echo -e "\t\t- Script: ${FUNCTION_PATH}/../utils/pack_file.ksh"
        echo -e "\t\t- Restart File: $(readlink -f ${ofile})"
        echo -e "\t\t- Mask file: $(readlink -f ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc)"
        echo -e "\n\t\t - NOTE: Errors with >VARDIMIDS (value 2) is out of range (1:1)< can be ignored (maybe?)"
        ${FUNCTION_PATH}/../utils/pack_file.ksh \
                $(basename $ofile) \
                $(basename ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc)

        echo -e "\t\t- Listing files after program execution:"
        ls -ratl

        mv 1d_$(basename $ofile) ${newest_jsbach_restart_file}
        ln -sf ${newest_jsbach_restart_file} ${COUPLE_DIR}/latest_jsbach_restart_file.nc
        # Go back to the previous directory
        cd ${startdir}
}

set_glacial_mask_jsbach_update_restart_veg() {
        echo; echo "                *   updating active and bare plant fractions in jsbach restart file >>veg<<"
        echo -e "\t\t-   Making a temporary directory to work in"
        startdir=$(pwd)
        if [ -d ${COUPLE_DIR}/update_veg ]; then
                # Throw away the temporary working directory from last time to
                # avoid merge problems in the next coupling iteration:
                rm -rf ${COUPLE_DIR}/update_veg
        fi
        mkdir -p ${COUPLE_DIR}/update_veg; cd ${COUPLE_DIR}/update_veg

        # TODO: This next block is repeated code, it should be made it's own function
        echo -e "\t\t-   determining where the glacial mask has advanced and retreated"

        # === Fix for iterative coupling: use fixed reference file instead of updating jsbach.nc ===
        JSBACH_REFERENCE_FILE=${COUPLE_DIR}/jsbach_original_reference.nc
        if [ ! -f ${JSBACH_REFERENCE_FILE} ]; then
            echo -e "\t\t-   Creating reference file for future comparisons: ${JSBACH_REFERENCE_FILE}"
            cp ${FORCING_DIR_jsbach}/jsbach.nc ${JSBACH_REFERENCE_FILE}
        fi
        # === End of fix ===

        echo -e "\t\t-   New JSBACH Init File: ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc"
        echo -e "\t\t-   Reference JSBACH File: ${JSBACH_REFERENCE_FILE}"
        cdo -s -selvar,glac -sub \
            ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc \
            ${JSBACH_REFERENCE_FILE} \
            ${COUPLE_DIR}/update_veg/jsbach_glac_diff.nc 2>> ${COUPLE_DIR}/update_veg/cdo_stderr_ice2echam
        # Make two masks: one for glacial retreat, and one for glacial advance
        cdo -s gtc,0 ${COUPLE_DIR}/update_veg/jsbach_glac_diff.nc ${COUPLE_DIR}/update_veg/new_glaciers.nc
        cdo -s gtc,0 -mulc,-1 ${COUPLE_DIR}/update_veg/jsbach_glac_diff.nc ${COUPLE_DIR}/update_veg/new_non_glaciers.nc
        sum_new_glaciers=$(cdo -s -output -fldsum -gtc,0 ${COUPLE_DIR}/update_veg/jsbach_glac_diff.nc | tr -d "[:blank:]")
        sum_new_land=$(cdo -s -output -fldsum -gtc,0 -mulc,-1 ${COUPLE_DIR}/update_veg/jsbach_glac_diff.nc | tr -d "[:blank:]")
        echo -e "\t\t- Total number of new glacial cells: >>>$sum_new_glaciers<<<"
        echo -e "\t\t- Total number of new deglacial cells: >>>$sum_new_land<<<"

        if [[ $sum_new_glaciers -eq 0 ]] && [[ $sum_new_land -eq 0 ]]; then
                echo -e "\t\t- >>>INFO<<< No new land or new glaciers, this routine will be skipped..."
                newest_jsbach_restart_file=$(readlink ${RESTART_DIR_jsbach}/restart_${EXP_ID}_veg.nc)
                ln -sf ${newest_jsbach_restart_file} ${COUPLE_DIR}/latest_veg_restart_file.nc
                return
        fi

        newest_jsbach_restart_file=$(readlink ${RESTART_DIR_jsbach}/restart_${EXP_ID}_veg.nc)
        # Make a backup copy
        cp ${RESTART_DIR_jsbach}/restart_${EXP_ID}_veg.nc ${RESTART_DIR_jsbach}/restart_${EXP_ID}_veg_backup_$(date +%Y%m%d%H%M%S).nc
        echo -e "\t\t- Unpacking JSBACH restart file: $newest_jsbach_restart_file"
        # Ensure that the unpack binary is here:
        if [ ! -f unpack ]; then
                #if [ -f ${FUNCTION_PATH}..//utils/_build/unpack ]; then
                #        cp -v ${FUNCTION_PATH}..//utils/_build/unpack ./
                if [ -f ${FUNCTION_PATH}/../utils/unpack ]; then
                        cp -v ${FUNCTION_PATH}/../utils/unpack ./
                else
                        echo -e "\t\t E R O R R --- unpack binary missing, don't know how to proceed"
                        exit
                fi
        fi
        # Copy the files here:
        cp ${newest_jsbach_restart_file} ./
        cp ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc ./
        echo -e "\t\t- Running unpack binary with arguments:"
        echo -e "\t\t- Script: ${FUNCTION_PATH}/../utils/unpack_file.ksh"
        echo -e "\t\t- Restart File: $(basename ${newest_jsbach_restart_file})"
        echo -e "\t\t- Mask file: $(basename ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc)"
        ${FUNCTION_PATH}/../utils/unpack_file.ksh \
                $(basename $newest_jsbach_restart_file) \
                $(basename ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc) \
                2>&1 > ${COUPLE_DIR}/unpack_veg_script_stdout_stderr

        echo -e "\t\t- Listing files after program execution:"
        ls -ratl
        restartfile=2d_$(basename ${newest_jsbach_restart_file})
        # 'cdo showvar' does not work for variables with more than three dimensiones
        echo "ncdump -h ${restartfile} | grep 'double ' | cut -f 2 -d ' ' | cut -f1 -d '('"
        variables=$(ncdump -h ${restartfile} | grep 'double ' | cut -f 2 -d ' ' | cut -f1 -d '(')
        echo $variables

        for var in ${variables}; do
                ncks -O --no_abc \
                        -v ${var} \
                        2d_$(basename ${newest_jsbach_restart_file}) \
                        2d_restart_${EXP_ID}_veg_${var}.nc

                # The following can be turned into an exit if needed:
                if [ ! -f 2d_restart_${EXP_ID}_veg_${var}.nc ]; then
                        echo -e "\t\t !!! WARNING !!! -- no file for ${var}"
                        echo -e "\t\t                    continue loop"
                        continue
                fi

                # Keep the directory clean: make one directory for each variable:
                here=$(pwd); mkdir -p ${var}; mv 2d_restart_${EXP_ID}_veg_${var}.nc ${var}; cd ${var}
                echo -e "\t\t- Updating JSBACH restart file -- ${var}\n"
			case ${var} in
				desert_fpc )
					# Skip desert_fpc modification - let dynamic vegetation handle it
					# The original cdo ifthenelse and ncap2 methods both caused file corruption
					echo -e "\t\t- Skipping desert_fpc modification (dynamic vegetation will handle)"
					;;
				bare_fpc )
					# Skip bare_fpc modification - cdo ifthenelse causes dimension errors
					# that result in pack failure and missing variables in final restart file
					# Dynamic vegetation will handle bare soil fraction
					echo -e "\t\t- Skipping bare_fpc modification (cdo causes dimension errors, dynveg will handle)"
					;;
				act_fpc )
					# Skip act_fpc modification - cdo splitlevel/merge causes dimension errors
					# that result in pack failure and missing variables in final restart file
					# Dynamic vegetation will handle active plant cover fraction
					echo -e "\t\t- Skipping act_fpc modification (cdo causes dimension errors, dynveg will handle)"
					;;
			esac
                cd $here
        done
        mergelist=""
        for var in ${variables}; do
                mergelist="$mergelist ${var}/2d_restart_${EXP_ID}_veg_${var}.nc"
        done
        echo -e "\t\t- Creating merge list"
        # For debugging:
        #echo "Done constructing merge list and renaming time to tiles, these files will be merged"
        #for f in $mergelist; do
        #        ls -l $f
        #done
        # Merge with ncks, since cdo has issues:
        echo -e "\t\t- Merging with ncks, as cdo has problems"
        ofile=2d_restart_${EXP_ID}_veg_merged.nc
        for f in $mergelist; do
                ncks -A -h $f $ofile
        done
        echo -e "\t\t- Finished with modification of variables, now we need to repack"
        # Ensure that the pack binary is here:
        if [ ! -f pack ]; then
                if [ -f ${FUNCTION_PATH}/../utils/pack ]; then
                        cp -v ${FUNCTION_PATH}/../utils/pack .
                else
                        echo -e "\t\t E R O R R --- pack binary missing, don't know how to proceed"
                        exit
                fi
        fi
        # Copy the files here:
        cp ${newest_jsbach_restart_file} .
        cp ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc .
        echo -e "\t\t- Running pack binary with arguments:"
        echo -e "\t\t- Script: ${FUNCTION_PATH}/../utils/pack_file.ksh"
        echo -e "\t\t- Restart File: $(readlink -f ${ofile})"
        echo -e "\t\t- Mask file: $(readlink -f ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc)"
        echo -e "\n\t\t - NOTE: Errors with >VARDIMIDS (value 2) is out of range (1:1)< can be ignored (maybe?)"
        ${FUNCTION_PATH}/../utils/pack_file.ksh \
                $(basename $ofile) \
                $(basename ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc)

        echo -e "\t\t- Listing files after program execution:"
        ls -ratl

        mv 1d_$(basename $ofile) ${newest_jsbach_restart_file}
        ln -sf ${newest_jsbach_restart_file} ${COUPLE_DIR}/latest_veg_restart_file.nc
        
        # === Fix desert_fpc in deglaciated areas ===
        # The Python script directly modifies the 1D veg restart file to set desert_fpc=0
        # in deglaciated areas, allowing vegetation to grow. This avoids pack/unpack issues.
        echo -e "\t\t- Running fix_desert_fpc_safe.py to enable vegetation recovery"
        if [ -f ${FUNCTION_PATH}/../utils/fix_desert_fpc_safe.py ]; then
            # Load python environment with numpy/netCDF4 (mambaforge default has both)
            /sw/spack-levante/mambaforge-22.9.0-2-Linux-x86_64-kptncg/bin/python3 ${FUNCTION_PATH}/../utils/fix_desert_fpc_safe.py \
                ${newest_jsbach_restart_file} \
                ${JSBACH_REFERENCE_FILE} \
                ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc
            if [ $? -eq 0 ]; then
                echo -e "\t\t- desert_fpc fix applied successfully"
            else
                echo -e "\t\t- WARNING: desert_fpc fix failed, but continuing (dynveg may still work)"
            fi
        else
            echo -e "\t\t- WARNING: fix_desert_fpc_safe.py not found, skipping desert_fpc fix"
        fi
        # === End of desert_fpc fix ===
        
        # === FIX 2025-01-24: Sync pot_fpc with cover_fract ===
        # The pot_fpc in veg restart must match cover_fract in jsbach restart
        # Otherwise JSBACH will crash with memory errors (malloc invalid size)
        echo -e "\t\t- Syncing pot_fpc in veg restart with cover_fract in jsbach restart..."
        if [ -f ${FUNCTION_PATH}/../utils/fix_veg_restart_pot_fpc.py ]; then
            newest_jsbach_restart_file_main=$(readlink ${RESTART_DIR_jsbach}/restart_${EXP_ID}_jsbach.nc)
            /sw/spack-levante/mambaforge-22.9.0-2-Linux-x86_64-kptncg/bin/python3 \
                ${FUNCTION_PATH}/../utils/fix_veg_restart_pot_fpc.py \
                ${newest_jsbach_restart_file} \
                ${newest_jsbach_restart_file_main} \
                >> ${COUPLE_DIR}/update_veg/pot_fpc_sync_stdout_stderr 2>&1
            if [ $? -eq 0 ]; then
                echo -e "\t\t- pot_fpc synced successfully"
            else
                echo -e "\t\t- WARNING: pot_fpc sync failed, continuing anyway"
            fi
        else
            echo -e "\t\t- WARNING: fix_veg_restart_pot_fpc.py not found"
        fi
        # === End pot_fpc sync ===
        
        # === FIX 2025-01-26: Initialize act_fpc for deglaciated points ===
        # When ice melts, act_fpc remains 0 at deglaciated points because vegetation
        # hasn't grown yet. This causes fpc_to_cover_fract_pot to produce near-zero
        # cover_fract_pot, triggering JSBACH fatal error "sum of cover_fract_inst /= 1".
        # Solution: Initialize act_fpc using pot_fpc values at deglaciated points.
        echo -e "\t\t- Initializing act_fpc for deglaciated points..."
        if [ -f ${FUNCTION_PATH}/../utils/fix_veg_restart_act_fpc.py ]; then
            /sw/spack-levante/mambaforge-22.9.0-2-Linux-x86_64-kptncg/bin/python3 \
                ${FUNCTION_PATH}/../utils/fix_veg_restart_act_fpc.py \
                ${newest_jsbach_restart_file} \
                >> ${COUPLE_DIR}/update_veg/act_fpc_fix_stdout_stderr 2>&1
            if [ $? -eq 0 ]; then
                echo -e "\t\t- act_fpc initialized successfully"
            else
                echo -e "\t\t- WARNING: act_fpc initialization failed, continuing anyway"
            fi
        else
            echo -e "\t\t- WARNING: fix_veg_restart_act_fpc.py not found"
        fi
        # === End act_fpc initialization ===
        
        # === FIX 2025-01-24: Sync veg restart to WORK_DIR ===
        # esm_tools copies restart files to WORK_DIR before coupling scripts run,
        # so our modifications in restart directory are not reflected in WORK_DIR.
        # We must explicitly copy the fixed file to WORK_DIR.
        if [ -n "${WORK_DIR}" ] && [ -d "${WORK_DIR}" ]; then
            echo -e "\t\t- Syncing fixed veg restart to WORK_DIR..."
            if [ -f "${newest_jsbach_restart_file}" ]; then
                cp -v "${newest_jsbach_restart_file}" "${WORK_DIR}/restart_${EXP_ID}_veg.nc"
                echo -e "\t\t- veg restart synced to WORK_DIR successfully"
            fi
        fi
        # === End WORK_DIR sync ===
        
        # Go back to the previous directory
        cd ${startdir}
}

## @fn ensure_echam_restart_aliases_for_m_streams()
## ECHAM/JSBACH sometimes expects rerun/restart files with "*m" stream suffixes (jsbachm, vegm, surfm, co2m, yassom, ...)
## but our iterative coupling produces/updates the non-m restart files (restart_${EXP_ID}_jsbach.nc, restart_${EXP_ID}_veg.nc, ...).
## If the *m files are missing, ECHAM/JSBACH may initialize vegetation fractions with tiny defaults (1e-10) and crash.
##
## FIX 2025-01-25: Added logic to copy monthly restart files (jsbachm, landm, vegm, etc.) from previous
## run directory. These files contain variables like box_veg_ratio, veg_fract_correction, gfw that are
## not present in the standard jsbach/veg restart files. Without this, ECHAM crashes with:
##   IO_read_stream: failed to read from rerun file: box_veg_ratio
##   malloc(): invalid size (unsorted)
ensure_echam_restart_aliases_for_m_streams() {
        if [ -z "${WORK_DIR}" ] || [ ! -d "${WORK_DIR}" ]; then
                echo "                -   WARNING: WORK_DIR not set or not a directory; cannot create restart alias files"
                return 0
        fi
        echo "                -   Ensuring rerun/restart alias files for *m streams exist in WORK_DIR"
        cd "${WORK_DIR}"

        # === FIX 2025-01-25: First try to copy monthly restart files from previous run ===
        # The *m restart files (jsbachm, landm, vegm, etc.) contain accumulated/diagnostic variables
        # that are NOT present in the base restart files. We must copy them from the previous run.
        #
        # Derive experiment root directory from WORK_DIR:
        # WORK_DIR = /path/to/exp/run_awiesm_YYYYMMDD-YYYYMMDD/work
        # exp_root = /path/to/exp
        local current_run_dir=$(dirname "${WORK_DIR}")
        local exp_root=$(dirname "${current_run_dir}")
        
        local prev_run_root=""
        for cand in $(ls -dt ${exp_root}/run_awiesm_* 2>/dev/null); do
                # Skip current run directory
                if [ "$(realpath "${cand}")" = "$(realpath "${current_run_dir}")" ]; then
                        continue
                fi
                # Check if this candidate has the monthly restart files we need
                if [ -d "${cand}/work" ]; then
                        # Look for jsbachm file (with or without timestamp)
                        local jsbachm_file=$(ls ${cand}/work/restart_${EXP_ID}*_jsbachm.nc 2>/dev/null | tail -1)
                        if [ -n "${jsbachm_file}" ] && [ -s "${jsbachm_file}" ]; then
                                prev_run_root="${cand}"
                                break
                        fi
                fi
        done

        if [ -n "${prev_run_root}" ] && [ -d "${prev_run_root}/work" ]; then
                echo "                -   Found previous run with monthly restarts: ${prev_run_root}"
                # Copy monthly restart files from previous run
                # NOTE: Do NOT copy jsbachm and vegm! These files contain accumulated monthly mean values
                # (e.g., cover_fract sum = 450 instead of 1), not instantaneous restart values.
                # Using them as restart causes JSBACH to fail with "sum of cover_fract_inst /= 1" error.
                # The fallback logic below will correctly link jsbachm -> jsbach.
                for mstream in landm surfm co2m yassom spm glm g3bm; do
                        if [ ! -e "restart_${EXP_ID}_${mstream}.nc" ]; then
                                # Find the latest file (may have timestamp in name)
                                local src_file=$(ls ${prev_run_root}/work/restart_${EXP_ID}*_${mstream}.nc 2>/dev/null | tail -1)
                                if [ -n "${src_file}" ] && [ -s "${src_file}" ]; then
                                        cp "${src_file}" "restart_${EXP_ID}_${mstream}.nc"
                                        echo "                    --> copied ${mstream} from previous run"
                                fi
                        fi
                done
        else
                echo "                -   No previous run with monthly restarts found; will use fallback aliases"
        fi
        # === End FIX 2025-01-25 ===

        # Core component restarts: create *m aliases pointing to existing non-m files (ONLY if not already present)
        for base in co2 surf jsbach veg yasso; do
                if [ -f "restart_${EXP_ID}_${base}.nc" ] && [ ! -e "restart_${EXP_ID}_${base}m.nc" ]; then
                        ln -sf "restart_${EXP_ID}_${base}.nc" "restart_${EXP_ID}_${base}m.nc"
                        echo "                    --> linked restart_${EXP_ID}_${base}m.nc -> restart_${EXP_ID}_${base}.nc"
                fi
        done

        # Atmosphere diagnostic streams: link to main echam restart if specific files are missing
        if [ -f "restart_${EXP_ID}_echam.nc" ]; then
                for s in spm glm g3bm; do
                        if [ ! -e "restart_${EXP_ID}_${s}.nc" ]; then
                                ln -sf "restart_${EXP_ID}_echam.nc" "restart_${EXP_ID}_${s}.nc"
                                echo "                    --> linked restart_${EXP_ID}_${s}.nc -> restart_${EXP_ID}_echam.nc"
                        fi
                done
                for mm in monthly_means_2d monthly_means_3d_g3b monthly_means_3d_gl monthly_means_3d_sp; do
                        if [ ! -e "restart_${EXP_ID}_${mm}.nc" ]; then
                                ln -sf "restart_${EXP_ID}_echam.nc" "restart_${EXP_ID}_${mm}.nc"
                                echo "                    --> linked restart_${EXP_ID}_${mm}.nc -> restart_${EXP_ID}_echam.nc"
                        fi
                done
        fi

        # Some setups reference landm; if missing, link it to jsbach restart as a best-effort fallback.
        # NOTE: This is a FALLBACK only - the real landm file should have been copied above.
        if [ -f "restart_${EXP_ID}_jsbach.nc" ] && [ ! -e "restart_${EXP_ID}_landm.nc" ]; then
                ln -sf "restart_${EXP_ID}_jsbach.nc" "restart_${EXP_ID}_landm.nc"
                echo "                    --> linked restart_${EXP_ID}_landm.nc -> restart_${EXP_ID}_jsbach.nc (FALLBACK)"
        fi
}

## @fn ensure_nonempty_grib_rerun_files_for_resume_month()
## ECHAM/JSBACH in this setup uses GRIB1 rerun files like:
##   ${EXP_ID}_YYYYMM.01_echam, ${EXP_ID}_YYYYMM.01_jsbach, ...
## On a resume at YYYY-01-01 the YYYY01.01_* files may exist but be 0-byte (created early, never written),
## which leads to rerun read failures (st/svo/lsp/sd, etc.) and can crash with malloc() corruption.
## Workaround: before compute starts, seed the YYYY01.01_* files with the previous month (YYYY-1)12.01_* contents.
ensure_nonempty_grib_rerun_files_for_resume_month() {
        if [ -z "${WORK_DIR}" ] || [ ! -d "${WORK_DIR}" ]; then
                echo "                -   WARNING: WORK_DIR not set or not a directory; cannot seed GRIB rerun files"
                return 0
        fi

        # Derive current YYYYMM from CHUNK_START_DATE_echam (format: YYYY-MM-DDT...)
        year=${CHUNK_START_DATE_echam:0:4}
        month=${CHUNK_START_DATE_echam:5:2}
        yyyymm="${year}${month}"

        # Previous month YYYYMM (handles Jan -> Dec of previous year)
        if [ "${month}" = "01" ]; then
                prev_year=$((10#${year} - 1))
                prev_month="12"
        else
                prev_year=$((10#${year}))
                prev_month=$(printf "%02d" $((10#${month} - 1)))
        fi
        prev_yyyymm="${prev_year}${prev_month}"

        # Find a previous run directory that actually contains the previous-month GRIB files we need.
        # (Don't rely on string matching the current run name; instead, search by file presence.)
        prev_run_root=""
        for cand in $(ls -dt ${BASE_DIR}/${EXP_ID}/run_awiesm_* 2>/dev/null); do
                if [ -d "${cand}/work" ] && [ -s "${cand}/work/${EXP_ID}_${prev_yyyymm}.01_echam" ]; then
                        prev_run_root="${cand}"
                        break
                fi
        done

        echo "                -   Seeding GRIB rerun files for ${yyyymm} from previous month ${prev_yyyymm}"
        echo "                    current run: $(dirname "${WORK_DIR}")"
        echo "                    previous run: ${prev_run_root}"

        if [ -z "${prev_run_root}" ] || [ ! -d "${prev_run_root}/work" ]; then
                echo "                    WARNING: previous run directory not found; cannot seed GRIB rerun files"
                return 0
        fi

        # Streams to seed (core + monthly means); add more if needed
        for stream in echam jsbach veg land surf co2 yasso accw ism monthly_means_2d monthly_means_3d_g3b monthly_means_3d_gl monthly_means_3d_sp; do
                src="${prev_run_root}/work/${EXP_ID}_${prev_yyyymm}.01_${stream}"
                dst="${WORK_DIR}/${EXP_ID}_${yyyymm}.01_${stream}"
                prev_alias="${WORK_DIR}/${EXP_ID}_${prev_yyyymm}.01_${stream}"

                # Only seed if dst is missing or zero-byte, and src exists and is non-empty
                if [ ! -s "${dst}" ] && [ -s "${src}" ]; then
                        # IMPORTANT:
                        # Do NOT rewrite the internal GRIB timestamp here. For this setup, ECHAM resumes at YYYY-01-01
                        # but reads rerun fields at the last time step of the previous run (e.g. YYYY-1-12-31 23:52).
                        # Re-encoding timestamps can cause IO_read_stream failures and malloc() crashes.
                        
                        # Use cat instead of cp to avoid potential attribute/filesystem issues
                        cat "${src}" > "${dst}"
                        echo "                    --> seeded ${dst} from ${src} (size: $(stat -c %s ${dst} 2>/dev/null || ls -sh ${dst}))"
                        
                        # Verify copy succeeded
                        if [ ! -s "${dst}" ]; then
                            echo "                    !!! ERROR: Failed to seed ${dst} (0 bytes). Disk full?"
                            # Try one more time with cp
                            cp "${src}" "${dst}"
                            if [ ! -s "${dst}" ]; then
                                echo "                    !!! FATAL: Unable to create non-empty rerun file."
                                exit 1
                            fi
                        fi
                fi

                # Create an alias with previous-month filename in WORK_DIR pointing to the seeded current-month file.
                # Some ECHAM setups may still look for the previous-month filename during resume.
                if [ -s "${dst}" ] && [ ! -e "${prev_alias}" ]; then
                        ln -sf "${dst##*/}" "${prev_alias}"
                        echo "                    --> linked ${prev_alias} -> ${dst##*/}"
                fi

                # Also seed the .codes file if present (helps readability / variable table)
                src_codes="${src}.codes"
                dst_codes="${dst}.codes"
                prev_alias_codes="${prev_alias}.codes"
                if [ ! -s "${dst_codes}" ] && [ -s "${src_codes}" ]; then
                        cp -v "${src_codes}" "${dst_codes}"
                        echo "                    --> seeded ${dst_codes} from ${src_codes}"
                fi
                if [ -s "${dst_codes}" ] && [ ! -e "${prev_alias_codes}" ]; then
                        ln -sf "${dst_codes##*/}" "${prev_alias_codes}"
                        echo "                    --> linked ${prev_alias_codes} -> ${dst_codes##*/}"
                fi
        done
}
## @fn update_land_runoff_select_glacial_discharge() 
update_land_runoff_select_glacial_discharge() {
        echo; echo "                *   selecting glacial discharge from ice_file_for_atmosphere.nc"
        cdo -s -selvar,${ice_mass_loss_varname} ${COUPLE_DIR}/ice_file_for_atmosphere.nc ${COUPLE_DIR}/ice_discharge.nc
}


## @fn update_land_runoff_fill_missing() 
update_land_runoff_fill_missing() { 
        echo; echo "                *   filling non-pism areas with 0"
        ifile=${COUPLE_DIR}/ice_discharge_${RES_echam}grid.nc
        ofile=${COUPLE_DIR}/gfw_atmo.nc
        change_lonlat $ifile $ice_mass_loss_varname
        # multiply by -1 to get signs correct (in ice sheet, ice mass gain is positive, loss is negative
        # in the atmosphere, loss should be positive, gain should be negative)
        # PG: Fix for way too much water:
        # Divide by minus 1000; oasis gives m/s to FESOM, not kg/m2/s.
        cdo -s mulc,-1 -divc,1000 -setmisstoc,0 $ifile $ofile
}

distribute_total_discharge_over_ocean() {
#        echo; echo "                *   integrate discharge over echam grid"
        ifile1=${COUPLE_DIR}/ice_discharge_${RES_echam}grid.nc
        ifile2=${COUPLE_DIR}/ice_discharge_${RES_echam}grid.gridarea.nc
#        ofile=${COUPLE_DIR}/ice_discharge_${RES_echam}grid.fldsum.nc
        cdo -s gridarea ${ifile1} ${ifile2}
#        cdo -s fldsum -mul ${ifile1} ${ifile2} ${ofile}
        CELL_AREA_FESOM_FILE=${CELL_AREA_FESOM_FILE:-"fesom.mesh.diag.nc"}
        echo; echo "            * Start apply_hosing_correction.py with: $ifile1 $ifile2 ${MESH_PATH_FESOM}/${CELL_AREA_FESOM_FILE} ${HOSING_FILE_LANDICE_LOSS} "
        module load python3
        echo " * HOSING_FILE_LANDICE_LOSS = ${HOSING_FILE_LANDICE_LOSS}"
        if [ -f "${HOSING_FILE_LANDICE_LOSS}/landice_yearly_mass_loss.out" ]; then
            rm ${HOSING_FILE_LANDICE_LOSS}/landice_yearly_mass_loss.out
        fi

        python ${FUNCTION_PATH}/../utils/apply_hosing_correction.py "$ifile1" "$ifile2" "${MESH_PATH_FESOM}/${CELL_AREA_FESOM_FILE}" "${HOSING_FILE_LANDICE_LOSS}"
        module unload python3
        echo; echo "            * Finished apply_hosing_correction.py "
#
#        
#        ifile=${COUPLE_DIR}/ice_discharge_${RES_echam}grid.fldsum.nc
#        ofile=${COUPLE_DIR}/ice_discharge_${RES_echam}grid.fldsum.flux.nc
#        echo; echo "                *   evenly spread discharge over ocean with total node area = ${TOTAL_NODE_AREA}"
#        cdo -f divc,${TOTAL_NODE_AREA} ${ifile} ${ofile} 
#        TOTAL_DISCHARGE=$( ncdump ${ofile} | awk '/total_ice_mass_loss_flux =/{getline; print}' | tr -d ";" | xargs )
#        
#        echo; echo "                *   write to ${HOSING_FILE_LANDICE_LOSS}"i
#        echo 1 > ${HOSING_FILE_LANDICE_LOSS}
#        echo ${TOTAL_DISCHARGE} >> ${HOSING_FILE_LANDICE_LOSS}
}

## @fn update_land_runoff_prepare_for_oasis() 
update_land_runoff_prepare_for_oasis() {
        echo; echo "                *   preparing file for oasis"
        ifile=${COUPLE_DIR}/gfw_atmo.nc
        ofile=${COUPLE_DIR}/gfw_atmo_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc 
        ncrename -v total_ice_mass_loss_flux,gfw_atmo $ifile tmp; mv tmp $ifile
        ncrename -d lon,nx $ifile tmp; mv tmp $ifile
        ncrename -d lat,ny $ifile tmp; mv tmp $ifile
        ncks -A ${FUNCTION_PATH}/../utils/axis_366.nc $ifile
        ncpdq -a time,ny,nx $ifile tmp; mv tmp $ifile
        mv $ifile $ofile
        cp ${COUPLE_DIR}/gfw_atmo_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc ${WORK_DIR}/gfw_atmo.nc 
        ln -sf ${COUPLE_DIR}/gfw_atmo_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc ${COUPLE_DIR}/latest_gfw_atmo.nc 
} 

## @fn update_land_runoff_set_namcouple_override() 
update_land_runoff_set_namcouple_override() {
        echo; echo "                *   making a file with instructions to override namcouple"
        switch_file=${COUPLE_DIR}/namcouple_override_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.dat
        echo "" > $switch_file
        echo oasis3mct_add_general_info_to_namcouple 12 2 fesom echam6 \$RUNTIME_awicm >> $switch_file
}

## @fn update_land_runoff_set_namelist_modifications_next_echam_run() 
update_land_runoff_set_namelist_modifications_next_echam_run() {
        echam_variable_modify_file=${COUPLE_DIR}/echam_namelist_switches_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.dat  
        echo; echo "                *   setting echam variables for land runoff to be read from file in next run"
              echo "                -   echam_variable_modify_file=$echam_variable_modify_file"
        if [ ! -s $echam_variable_modify_file ]; then
                :> $echam_variable_modify_file
        fi
        # NOTE:
        # Enabling lgfw requires a consistent restart/rerun state for the gfw field.
        # In this experiment setup the gfw state is not present in the standard rerun/restart files,
        # and turning this on causes ECHAM to abort at resume with:
        #   IO_read_stream: failed to read from rerun file: gfw  -> malloc() crash
        # Therefore we keep lgfw disabled here. If gfw coupling is needed in the future,
        # implement a proper gfw rerun/restart provision first.
        echo "lgfw=0" > ${COUPLE_DIR}/oasis3mct_config_switches.dat
        # BUG: Moved to the function below due to bad logic in echam.functions
        # add_to ${echam_variable_modify_file} echam_namelist_switches.dat config echam
        add_to ${COUPLE_DIR}/oasis3mct_config_switches.dat oasis3mct_config_switches.dat config oasis3mct
}

## @fn update_land_runoff_set_namelist_modifications_next_jsbach_run() 
update_land_runoff_set_namelist_modifications_next_jsbach_run() { 
        jsbach_variable_modify_file=${COUPLE_DIR}/jsbach_namelist_switches_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.dat  
        echo; echo "                *   setting jsbach variables for land runoff to be read from file in next run"
              echo "                -   jsbach_variable_modify_file=$jsbach_variable_modify_file"
        if [ ! -s $jsbach_variable_modify_file ]; then
                :> $jsbach_variable_modify_file
        fi
        echo "submodelctl___lgfw___nml_entry=.FALSE." >> $jsbach_variable_modify_file
        # Yes, this is correct below: we want to modify the ECHAM namelist, but
        # need it to connect to JSBACH:
        echo "submodelctl___lgfw___nml_file=namelist.echam" >> $jsbach_variable_modify_file
        echo "hydrology_ctl___lgfw___nml_entry=.FALSE." >> $jsbach_variable_modify_file
        echo "hydrology_ctl___lgfw___nml_file=namelist.jsbach" >> $jsbach_variable_modify_file
        add_to ${jsbach_variable_modify_file} jsbach_namelist_switches.dat config jsbach
}

## @fn regrid_to_echam() 
regrid_to_echam() { 
        FILENAME=$1
        remapper=$2
        RES=${3:-${RES_echam}}
    if [[ "x$remapper" == "x" ]]; then
		remapper=remapcon2
	fi
	echo "                -   regridding $FILENAME to echam ${RES} grid with $remapper"
        cdo -s -P 28 \
                -$remapper,${RES}grid \
                -setgrid,${COUPLE_DIR}/ice.griddes \
                $FILENAME ${FILENAME%.*}_${RES}grid.nc
}

## @fn set_in_jsbach_restart_file() 
set_in_jsbach_restart_file() { 
        TARGET_VARNAME=$1
        SOURCE_FILENAME=$2
        SOURCE_VARNAME=$3
        unpack=$4
        RESTART_FILE=$(readlink ${RESTART_DIR_jsbach}/restart_${EXP_ID}_jsbach.nc)
        echo "                *   replacing: $TARGET_VARNAME "
        echo "                           in: $RESTART_FILE "
        echo "                         from: $SOURCE_VARNAME "
        echo "                           in: $SOURCE_FILENAME"
        if [ ! -z $unpack ]; then
                echo "                -   with unpacking from landpoint grid to latlon grid"
                python ${FUNCTION_PATH}/../utils/unpack_jsbach.py $TARGET_VARNAME $RESTART_FILE 
                ################################################################################ 
                # $ cdo -h replace
                # NAME
                #         replace - Replace variables
                #
                # SYNOPSIS
                #         replace  infile1 infile2 outfile
                #
                # DESCRIPTION
                #         The replace operator replaces variables in infile1 by variables from infile2 and write
                #         the result to outfile. Both input datasets need to have the same number of timesteps.
                ################################################################################ 
                cdo -s replace \
                        ${TARGET_VARNAME}_lonlat_grid.nc \
                        -chname,${SOURCE_VARNAME},${TARGET_VARNAME} \
                        -selvar,${SOURCE_VARNAME} ${SOURCE_FILENAME} \
                        ${TARGET_VARNAME}_lonlat_grid_replaced.nc
                echo "                -   and repacking from latlon grid to landpoint grid"
                # FIXME: The next line is probably wrong, or incomplete.
                python ${FUNCTION_PATH}/../utils/pack_jsbach.py ${TARGET_VARNAME} ${TARGET_VARNAME}_lonlat_grid_replaced.nc ${RESTART_FILE}
        else
               cdo -s replace \
                        ${RESTART_FILE} \
                        -chname,${SOURCE_VARNAME},${TARGET_VARNAME} \
                        -selvar,${SOURCE_VARNAME} ${SOURCE_FILENAME} \
                        tmp
               cp ${RESTART_FILE} ${RESTART_FILE%.*}_backup_before_replace_${TARGET_VARNAME}_$(date +%Y%m%d%H%M%S).nc
               mv tmp ${RESTART_FILE}
        fi
}

## @fn change_lonlat() 
change_lonlat() {
        file=$1
        varname=$2
        if ncdump -h $file | grep "${varname}" | grep -q "lon, lat"; then
                echo "                -   changing lat/lon dimension order for $varname in $file"
                ncpdq -a lat,lon $file tmp; mv tmp $file 2>> nco_stderr_ice2echam
        fi
}