diff --git a/couplings/echam/coupling_echam2ice.functions b/couplings/echam/coupling_echam2ice.functions index d65d54ebf..3b2fc3948 100644 --- a/couplings/echam/coupling_echam2ice.functions +++ b/couplings/echam/coupling_echam2ice.functions @@ -87,7 +87,7 @@ iterative_coupling_echam6_ice_make_forcing() { fi #add_to $(pwd)/atmosphere_file_for_ice.nc atmosphere_file_for_ice.nc couple cp $(pwd)/atmosphere_file_for_ice.nc ${COUPLE_DIR}/atmosphere_file_for_ice.nc - + cp $(pwd)/atmosphere_file_for_ice.nc ${COUPLE_DIR}/atmosphere_file_for_ice_${END_YEAR_echam}.nc echo; echo " ...done."; echo } diff --git a/couplings/echam/coupling_ice2echam.functions b/couplings/echam/coupling_ice2echam.functions index b2451607d..8688454d0 100644 --- a/couplings/echam/coupling_ice2echam.functions +++ b/couplings/echam/coupling_ice2echam.functions @@ -92,6 +92,18 @@ iterative_coupling_ice_echam6_update_glacial_mask() { 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 } @@ -426,7 +438,7 @@ update_orography_generate_target_orography_file() { 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}/ice_mask_before_T63grid.nc \ + -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" @@ -511,27 +523,13 @@ set_glacial_mask_echam_restart_file() { 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 - cp $ECHAM_RESTART_FILE echam_with_ice_mask.nc - - echo " - Preforming replacement" - - cdo -s -selvar,glac ${ECHAM_RESTART_FILE} just_glac.nc 2>> cdo_stderr_ice2echam - cdo -s -ifthenelse -setmisstoc,0 $dfile -chname,ice_mask,glac $ifile just_glac.nc just_glac_new.nc 2>> cdo_stderr_ice2echam - - ncdump just_glac_new.nc | sed s/int/double/g | ncgen -o tmp; mv tmp just_glac_new.nc - # This multipliation transforms "float" into "double" - ncap=${ncap:-ncap2} - $ncap -s "glac=glac*1.0000000000000;" \ - just_glac_new.nc \ - just_glac_new.tmp.nc - ncks -A just_glac_new.tmp.nc echam_with_ice_mask.nc - - #mv echam_with_ice_mask.nc ${ECHAM_RESTART_FILE} - ncatted -O -a history_of_appended_files,global,d,, echam_with_ice_mask.nc tmp.restart - ncatted -O -a history,global,d,, tmp.restart ${ECHAM_RESTART_FILE} + # 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 - rm just_glac_new.nc - rm just_glac_new.tmp.nc } ## @fn set_glacial_mask_echam_make_dummy_jan_surf() @@ -542,7 +540,15 @@ set_glacial_mask_echam_make_dummy_jan_surf() { # 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 - # Fix Lu: --------------------------- + # 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 @@ -550,9 +556,23 @@ set_glacial_mask_echam_make_dummy_jan_surf() { rm tmp fi # Fix Lu: --------------------------- - 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 + # 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} @@ -587,12 +607,20 @@ set_glacial_mask_echam_make_dummy_jan_surf() { 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 \ - -selvar,GLAC dummy_jan_surf.nc input_if_mask_false.nc + -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 \ @@ -604,6 +632,43 @@ set_glacial_mask_echam_make_dummy_jan_surf() { 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() @@ -772,6 +837,26 @@ EOF # 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 ############################################################################ @@ -786,6 +871,80 @@ EOF 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 @@ -795,11 +954,20 @@ set_glacial_mask_jsbach_update_restart_jsbach() { 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- Old JSBACH Init File: ${FORCING_DIR_jsbach}/jsbach.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 ${FORCING_DIR_jsbach}/jsbach.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 @@ -884,6 +1052,16 @@ set_glacial_mask_jsbach_update_restart_jsbach() { 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 # @@ -1052,11 +1230,20 @@ set_glacial_mask_jsbach_update_restart_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- Old JSBACH Init File: $(readlink ${FORCING_DIR_jsbach}/jsbach.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 \ - ${FORCING_DIR_jsbach}/jsbach.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 @@ -1124,83 +1311,25 @@ set_glacial_mask_jsbach_update_restart_veg() { # 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 - bare_fpc ) - # For bare fpc, we need the fraction of "bare" under glaciers to be 1 - cdo -s ifthenelse \ - ${COUPLE_DIR}/update_veg/new_glaciers.nc \ - ${COUPLE_DIR}/update_veg/new_glaciers.nc \ - 2d_restart_${EXP_ID}_veg_${var}.nc \ - 2d_restart_${EXP_ID}_veg_${var}.tmp - # Over newly deglaciated tiles, we put 95% bare soil - cdo -s ifthenelse \ - ${COUPLE_DIR}/update_veg/new_non_glaciers.nc \ - -mulc,0.95 ${COUPLE_DIR}/update_veg/new_non_glaciers.nc \ - 2d_restart_${EXP_ID}_veg_${var}.tmp \ - 2d_restart_${EXP_ID}_veg_${var}.nc - cdo -s -setvar,$var 2d_restart_${EXP_ID}_veg_${var}.nc ${var}.tmp - mv ${var}.tmp 2d_restart_${EXP_ID}_veg_${var}.nc - ;; - act_fpc ) - # NOTE: No correction is applied under - # new_glaciers, this is automatically taken - # care of by JSBACH internally. - # - # For active plant cover fraction, we put 0 on - # all tiles except tundra (tile 7) which gets - # 5%. - # - # Split up the levels - cdo -s splitlevel 2d_restart_${EXP_ID}_veg_${var}.nc 2d_restart_${EXP_ID}_veg_${var}_tile - mergelist="" - for tile in $(ls 2d_restart_${EXP_ID}_veg_${var}_tile0000??.nc); do - # NOTE: The period before the old - # dimension name, sfc, indicates - # optional presence of the dimension - # (ncrename -d renames dimensions). It - # seems as if the variable sometimes - # has the correct dimension name tiles - # (or maybe my test case was incorrect - # and had the dimension name sfc). To - # overcome this, we simply ensure that - # if sfc is the dimension name, it - # should be renamed to tiles, and if - # this is already the case, do nothing: - ncrename -d .sfc,tiles $tile - case ${tile} in - 2d_restart_${EXP_ID}_veg_${var}_tile000007.nc ) - echo -e "\t\t- Setting 5% initial cover of tundra for new deglacial cells" - # For the tundra tile - cdo -v -setvar,$var -ifthenelse \ - ${COUPLE_DIR}/update_veg/new_non_glaciers.nc \ - -mulc,0.05 ${COUPLE_DIR}/update_veg/new_non_glaciers.nc \ - 2d_restart_${EXP_ID}_veg_${var}_tile000007.nc \ - 2d_restart_${EXP_ID}_veg_${var}_tile000007.tmp - mv 2d_restart_${EXP_ID}_veg_${var}_tile000007.tmp 2d_restart_${EXP_ID}_veg_${var}_tile000007.nc - ;; - * ) - # For all other tiles: - cdo -s -setvar,$var -ifthenelse \ - ${COUPLE_DIR}/update_veg/new_non_glaciers.nc \ - -mulc,0.0 ${COUPLE_DIR}/update_veg/new_glaciers.nc \ - ${tile} ${tile}.tmp - mv ${tile}.tmp ${tile} - ;; - esac - mergelist="$mergelist $tile" - done - mkdir pre_merge - cp $mergelist pre_merge - echo -e "\t\t >>> PG DEBUG: Command for merging will be: cdo -v -O merge $mergelist ${var}.tmp" - echo -e "\t\t >>> PG DEBUG: The files in merglist ($mergelist) have been copied here:" - echo -e "\t\t >>> PG DEBUG: pre_merge -- directory contents:" - ls -ratlh pre_merge - cdo -v -O merge $mergelist ${var}.tmp - mv -v ${var}.tmp 2d_restart_${EXP_ID}_veg_${var}.nc - # NOTE: Same as above...not sure why this is needed twice?? - ncrename -d .sfc,tiles 2d_restart_${EXP_ID}_veg_${var}.nc - ;; - esac + 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="" @@ -1246,9 +1375,281 @@ set_glacial_mask_jsbach_update_restart_veg() { 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" @@ -1331,7 +1732,14 @@ update_land_runoff_set_namelist_modifications_next_echam_run() { if [ ! -s $echam_variable_modify_file ]; then :> $echam_variable_modify_file fi - echo "lgfw=1" > ${COUPLE_DIR}/oasis3mct_config_switches.dat + # 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 @@ -1345,11 +1753,11 @@ update_land_runoff_set_namelist_modifications_next_jsbach_run() { if [ ! -s $jsbach_variable_modify_file ]; then :> $jsbach_variable_modify_file fi - echo "submodelctl___lgfw___nml_entry=.TRUE." >> $jsbach_variable_modify_file + 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=.TRUE." >> $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 } diff --git a/couplings/echam/env_echam.py b/couplings/echam/env_echam.py index 637342e90..25a93b683 100644 --- a/couplings/echam/env_echam.py +++ b/couplings/echam/env_echam.py @@ -1,4 +1,30 @@ def prepare_environment(config): + # Get PISM domain from configuration + # Try multiple paths to find the domain: + # 1. From model2 config (iterative coupling) + # 2. From pism config directly + # 3. From echam config (user-specified override) + # 4. Default to "greenland" + pism_domain = None + + # Try model2 (iterative coupling configuration) + if "model2" in config: + model2_setup = config["model2"].get("setup_name", "pism") + if model2_setup in config: + pism_domain = config[model2_setup].get("domain") + + # Try pism directly + if not pism_domain and "pism" in config: + pism_domain = config["pism"].get("domain") + + # Try echam config (user override) + if not pism_domain: + pism_domain = config["echam"].get("pism_domain") + + # Default to greenland + if not pism_domain: + pism_domain = "greenland" + environment_dict = { "ICE_TO_ECHAM": int(config["general"]["first_run_in_chunk"]), "ECHAM_TO_ICE": int(config["general"]["last_run_in_chunk"]), @@ -7,7 +33,9 @@ def prepare_environment(config): "ISM_TO_ECHAM_update_glacial_mask": int(config["echam"].get("update_glacial_mask", True).__bool__()), "ISM_TO_ECHAM_update_land_runoff": 1, "COUPLE_DIR": config["general"]["experiment_couple_dir"], - "RES_echam": config["echam"]["resolution"], + "RES_echam": config["echam"]["resolution"], + # === FIX: Add DOMAIN_pism for Greenland GLAC update === + "DOMAIN_pism": pism_domain, "EXP_ID": config["general"]["command_line_config"]["expid"], "RESTART_DIR_echam": config["echam"]["experiment_restart_out_dir"], "DATA_DIR_echam": config["echam"]["experiment_outdata_dir"], diff --git a/couplings/pism/coupling_atmosphere2pism.functions b/couplings/pism/coupling_atmosphere2pism.functions index cc9472549..7524da295 100644 --- a/couplings/pism/coupling_atmosphere2pism.functions +++ b/couplings/pism/coupling_atmosphere2pism.functions @@ -85,7 +85,7 @@ atmosphere2pism() { which cdo echo -e "\t\t---> Preparing to couple generic atmosphere to specific ice sheet model PISM" ic_atm2pism_prepare - + #iterative_coupling_atmosphere_pism_regrid_method="INTERPOLATE" iterative_coupling_atmosphere_pism_regrid_method=${iterative_coupling_atmosphere_pism_regrid_method:-DOWNSCALE} bool_atm2pism_regrid_dyn_downscaling=${bool_atm2pism_regrid_dyn_downscaling:-0} echo -e "\t\t--> Regrid method >>> ${GREEN}${iterative_coupling_atmosphere_pism_regrid_method}${NOCOLOR} <<<" @@ -520,7 +520,7 @@ ic_atm2pism_regrid_dyn_downscale_generate_elevation_difference() { ncrename -v $atmosphere_name_elevation,usurf ${COUPLE_DIR}/elev_lo.nc ${COUPLE_DIR}/tmp.nc cdo settbounds,1mon -setreftime,1850-01-01,00:00:00 -setunit,m tmp.nc surface_altitude_coarse_resolution.nc - cleanup_list="$cleanup_list ${COUPLE_DIR}/elev_hi.nc ${COUPLE_DIR}/elev_lo.nc ${COUPLE_DIR}/elev_hi_minus_lo.nc" + #cleanup_list="$cleanup_list ${COUPLE_DIR}/elev_hi.nc ${COUPLE_DIR}/elev_lo.nc ${COUPLE_DIR}/elev_hi_minus_lo.nc" } ic_atm2pism_regrid_dyn_downscale_temperature_below_firn() { @@ -646,9 +646,9 @@ ic_atm2pism_regrid_downscale_split_names() { ${COUPLE_DIR}/atmosphere_file_for_ice_for_downscaling_ 2>> ${COUPLE_DIR}/cdo_stderr_atm2pism ##atmosphere_file_at_ice.${INTERPOL_TYPE_ATM}.nc # && rm atmosphere_file_for_ice.nc - for varname in $(cdo -s vardes ${COUPLE_DIR}/atmosphere_file_for_ice.nc | tr -s ' ' | awk -F ' ' '{print $2}'); do - cleanup_list="$cleanup_list ${COUPLE_DIR}/atmosphere_file_for_ice_for_downscaling_${varname}.nc" - done + #for varname in $(cdo -s vardes ${COUPLE_DIR}/atmosphere_file_for_ice.nc | tr -s ' ' | awk -F ' ' '{print $2}'); do + #cleanup_list="$cleanup_list ${COUPLE_DIR}/atmosphere_file_for_ice_for_downscaling_${varname}.nc" + #done } ## @fn ic_atm2pism_regrid_downscale_generate_elevation_difference() @@ -732,7 +732,9 @@ ic_atm2pism_regrid_downscale_temperature() { if [ "x${DOWNSCALE_TEMP}" == "x1" ]; then echo ""; echo " * downscaling temperature" DOWNSCALING_LAPSE_RATE=${DOWNSCALING_LAPSE_RATE:-"-0.005"} - if [ ${DOWNSCALING_LAPSE_RATE} -gt 0 ]; then + #if [ ${DOWNSCALING_LAPSE_RATE} -gt 0 ]; then + #if (( ${DOWNSCALING_LAPSE_RATE} > 0 )); then + if [ $(echo "${DOWNSCALING_LAPSE_RATE} > 0" | bc) -eq 1 ]; then echo "YOUR DOWNSCALING LAPSE RATE IS POSITIVE!" echo "IT WAS: $DOWNSCALING_LAPSE_RATE" exit 6 @@ -763,7 +765,9 @@ ic_atm2pism_regrid_downscale_temperature_below_firn() { if [ "x${DOWNSCALE_TEMP}" == "x1" ]; then echo ""; echo " * downscaling temperature_below_firn" DOWNSCALING_LAPSE_RATE=${DOWNSCALING_LAPSE_RATE:-"-0.005"} #old value = -0.007 - if [ ${DOWNSCALING_LAPSE_RATE} -gt 0 ]; then + #if (( ${DOWNSCALING_LAPSE_RATE} > 0 )); then + if [ $(echo "${DOWNSCALING_LAPSE_RATE} > 0" | bc) -eq 1 ]; then + #if [ ${DOWNSCALING_LAPSE_RATE} > 0 ]; then echo "YOUR DOWNSCALING LAPSE RATE IS POSITIVE!" echo "IT WAS: $DOWNSCALING_LAPSE_RATE" exit 6 diff --git a/couplings/pism/coupling_pism2atmosphere.functions b/couplings/pism/coupling_pism2atmosphere.functions index 7f86e3991..e45b6b1f3 100644 --- a/couplings/pism/coupling_pism2atmosphere.functions +++ b/couplings/pism/coupling_pism2atmosphere.functions @@ -108,9 +108,34 @@ function pism_atmosphere_get_newest_output { fi done - echo; echo " * Concatenating these files: ${PISM_OFILE_LIST} " + #echo; echo " * Concatenating these files: ${PISM_OFILE_LIST} " + #cdo cat ${PISM_OFILE_LIST} ${OUTPUT_DIR_pism}/latest_ex_file_pism.CHUNK.nc 2>> ${COUPLE_DIR}/cdo_stderr_pism2awiesm + #cdo timmean ${OUTPUT_DIR_pism}/latest_ex_file_pism.CHUNK.nc ${OUTPUT_DIR_pism}/latest_ex_file_pism.nc 2>> ${COUPLE_DIR}/cdo_stderr_pism2awiesm + echo; echo " * Concatenating these 100 files: ${PISM_OFILE_LIST} " cdo cat ${PISM_OFILE_LIST} ${OUTPUT_DIR_pism}/latest_ex_file_pism.CHUNK.nc 2>> ${COUPLE_DIR}/cdo_stderr_pism2awiesm - cdo timmean ${OUTPUT_DIR_pism}/latest_ex_file_pism.CHUNK.nc ${OUTPUT_DIR_pism}/latest_ex_file_pism.nc 2>> ${COUPLE_DIR}/cdo_stderr_pism2awiesm + + # --- 修改开始 --- + # 计算最后100年的起始年份 (例如: 结束年=2000, 2000-100+1 = 1901. 即 1901-2000) + AVG_START_YEAR_CALC=$((CHUNK_END_YEAR_pism - 100 + 1)) + + # 检查计算出的起始年是否早于该 chunk 的实际起始年 + # (防止 chunk 不足100年) + if [ ${AVG_START_YEAR_CALC} -lt ${CHUNK_START_YEAR_pism} ]; then + AVG_START_YEAR=${CHUNK_START_YEAR_pism} + echo " * Note: Chunk is shorter than 100 years. Averaging over available years." + else + AVG_START_YEAR=${AVG_START_YEAR_CALC} + fi + + echo " * Calculating time mean over the last 100 years (range: ${AVG_START_YEAR} to ${CHUNK_END_YEAR_pism})" + + # 在 timmean 之前使用 selyear 选择年份范围 + cdo timmean \ + -selyear,${AVG_START_YEAR}/${CHUNK_END_YEAR_pism} \ + ${OUTPUT_DIR_pism}/latest_ex_file_pism.CHUNK.nc \ + ${OUTPUT_DIR_pism}/latest_ex_file_pism.nc 2>> ${COUPLE_DIR}/cdo_stderr_pism2awiesm + # --- 修改结束 --- + #rm ${PISM_OFILE_LIST} ${PISM_cleanup_list} else diff --git a/couplings/utils/fix_desert_fpc_safe.py b/couplings/utils/fix_desert_fpc_safe.py new file mode 100644 index 000000000..ee1a7c6fc --- /dev/null +++ b/couplings/utils/fix_desert_fpc_safe.py @@ -0,0 +1,297 @@ +#!/usr/bin/env python3 +""" +Safe script to fix FPC issues in JSBACH veg restart file. + +This script directly modifies the veg restart file using netCDF4, +avoiding the pack/unpack process that causes precision issues. + +Two fixes are applied: +1. In deglaciated areas (where GLAC changed from 1 to 0): + - Set desert_fpc = 0.0 (allow vegetation to grow) + - Set bare_fpc = 1.0 - act_fpc_sum (to maintain sum = 1) + +2. For ALL landpoints where total FPC > 1: + - Reduce desert_fpc to bring total to 1 + - If not enough, also reduce bare_fpc + +Usage: + python fix_desert_fpc_safe.py [--dry-run] +""" + +import sys +import numpy as np + +try: + from netCDF4 import Dataset +except ImportError: + print("ERROR: netCDF4 module not found. Please load the appropriate module:") + print(" module load python3") + print(" or use: conda activate ") + sys.exit(1) + + +def get_landpoint_indices_from_2d_mask(mask_2d, landseamask): + """ + Convert 2D grid mask to 1D landpoint indices. + + Parameters: + ----------- + mask_2d : 2D array (lat, lon) with 1 where condition is true + landseamask : 2D array (lat, lon) with 1 for land points + + Returns: + -------- + Array of landpoint indices where mask_2d is 1 + """ + # Get all land points in order (row-major, i.e., lat varies slowest) + land_indices = np.where(landseamask.flatten() > 0.5)[0] + + # Get indices where mask is 1 + mask_flat = mask_2d.flatten() + + # Find which landpoints have mask = 1 + landpoint_indices = [] + for i, grid_idx in enumerate(land_indices): + if mask_flat[grid_idx] > 0.5: + landpoint_indices.append(i) + + return np.array(landpoint_indices) + + +def fix_all_fpc_issues(desert_fpc, bare_fpc, act_fpc_sum, tol=1e-10): + """ + Fix all landpoints where total FPC differs from 1 by more than tol. + + Strategy: + - If total > 1 + tol: + - First reduce desert_fpc to bring total to 1 + - If not enough, also reduce bare_fpc + - If total < 1 - tol: + - Increase bare_fpc by the deficit (bounded to [0, 1]) + """ + total_fpc = desert_fpc + bare_fpc + act_fpc_sum + problem_mask = np.abs(total_fpc - 1.0) > tol + problem_indices = np.where(problem_mask)[0] + + n_fixed = 0 + for idx in problem_indices: + diff = total_fpc[idx] - 1.0 + if diff > tol: + excess = diff + + # First try to reduce desert_fpc + if desert_fpc[idx] >= excess: + desert_fpc[idx] -= excess + n_fixed += 1 + else: + # Need to also reduce bare_fpc + remaining_excess = excess - desert_fpc[idx] + desert_fpc[idx] = 0.0 + + if bare_fpc[idx] >= remaining_excess: + bare_fpc[idx] -= remaining_excess + n_fixed += 1 + else: + # Fallback: enforce exact closure by setting bare to residual + bare_fpc[idx] = max(0.0, 1.0 - act_fpc_sum[idx]) + desert_fpc[idx] = 0.0 + n_fixed += 1 + elif diff < -tol: + deficit = -diff + # Prefer to fill deficit into bare_fpc (can't exceed 1 - act_fpc_sum) + max_bare = max(0.0, 1.0 - act_fpc_sum[idx]) + bare_fpc[idx] = min(max_bare, bare_fpc[idx] + deficit) + # If bare was already at max, adjust desert as residual (rare) + residual = 1.0 - (act_fpc_sum[idx] + bare_fpc[idx] + desert_fpc[idx]) + if abs(residual) > tol: + desert_fpc[idx] = min(1.0, max(0.0, desert_fpc[idx] + residual)) + n_fixed += 1 + + return desert_fpc, bare_fpc, n_fixed + + +def main(): + if len(sys.argv) < 4: + print(__doc__) + sys.exit(1) + + veg_restart_file = sys.argv[1] + jsbach_reference_file = sys.argv[2] + jsbach_new_file = sys.argv[3] + dry_run = '--dry-run' in sys.argv + + print(f"=== Fix FPC issues in veg restart ===") + print(f"Veg restart file: {veg_restart_file}") + print(f"Reference jsbach: {jsbach_reference_file}") + print(f"New jsbach: {jsbach_new_file}") + print(f"Dry run: {dry_run}") + print() + + # Load the glacier mask difference + print("Loading jsbach files to identify deglaciated areas...") + + with Dataset(jsbach_reference_file, 'r') as ref_nc: + glac_ref = ref_nc.variables['glac'][:] + if len(glac_ref.shape) == 3: + glac_ref = glac_ref[0, :, :] # Remove time dimension if present + # Get land sea mask from reference file + if 'slm' in ref_nc.variables: + landseamask = ref_nc.variables['slm'][:] + elif 'landseamask' in ref_nc.variables: + landseamask = ref_nc.variables['landseamask'][:] + else: + raise KeyError("Cannot find land sea mask (slm or landseamask) in reference file") + if len(landseamask.shape) == 3: + landseamask = landseamask[0, :, :] + + with Dataset(jsbach_new_file, 'r') as new_nc: + glac_new = new_nc.variables['glac'][:] + if len(glac_new.shape) == 3: + glac_new = glac_new[0, :, :] + + # Deglaciated areas: was glacier (ref=1), now land (new=0) + deglaciated_mask = (glac_ref > 0.5) & (glac_new < 0.5) + n_deglaciated_2d = np.sum(deglaciated_mask) + print(f"Number of deglaciated grid cells (2D): {n_deglaciated_2d}") + + # Load veg restart + print(f"\nLoading veg restart file: {veg_restart_file}") + + with Dataset(veg_restart_file, 'r+' if not dry_run else 'r') as veg_nc: + # Get current fpc values + desert_fpc = veg_nc.variables['desert_fpc'][:] + bare_fpc = veg_nc.variables['bare_fpc'][:] + + # Get act_fpc (has tiles dimension) + act_fpc = veg_nc.variables['act_fpc'][:] + # Get pot_fpc (has tiles dimension) if present + pot_fpc = veg_nc.variables['pot_fpc'][:] if 'pot_fpc' in veg_nc.variables else None + + # === Fix 0: Remove tiny values from act_fpc === + tiny_act = (act_fpc > 0) & (act_fpc < 1e-8) + n_tiny_act = np.sum(tiny_act) + if n_tiny_act > 0: + print(f"\n--- Fix 0: Tiny values in act_fpc ---") + print(f" Found {n_tiny_act} tiny values in act_fpc, setting to 0") + act_fpc[tiny_act] = 0.0 + if not dry_run: + veg_nc.variables['act_fpc'][:] = act_fpc + + # Also fix tiny values in desert_fpc + tiny_desert = (desert_fpc > 0) & (desert_fpc < 1e-8) + n_tiny_desert = np.sum(tiny_desert) + if n_tiny_desert > 0: + print(f" Found {n_tiny_desert} tiny values in desert_fpc, setting to 0") + desert_fpc[tiny_desert] = 0.0 + + act_fpc_sum = np.sum(act_fpc, axis=0) # Sum across tiles + + # === Fix 0b: Ensure pot_fpc is valid (sum across tiles must be 1) === + # JSBACH uses pot_fpc to derive cover_fract_pot and expects it to be a valid + # composition vector (sum=1) even when total vegetation fraction is small. + # Some restart files contain pot_fpc=0 for many landpoints, which triggers: + # FATAL ERROR in fpc_to_cover_fract_pot: sum of cover_fract_pot /= 1 + if pot_fpc is not None: + tol = 1e-10 + pot_sum = np.sum(pot_fpc, axis=0) + bad = np.where(np.abs(pot_sum - 1.0) > tol)[0] + if bad.size > 0: + print(f"\n--- Fix 0b: pot_fpc closure ---") + print(f" Found {bad.size} landpoints with |sum(pot_fpc)-1| > {tol:g}") + # Fix each bad landpoint: + for idx in bad: + s = pot_sum[idx] + if s > tol: + pot_fpc[:, idx] = pot_fpc[:, idx] / s + else: + # If pot_fpc is (near) all-zero, use act_fpc composition if available, + # otherwise put all potential vegetation into the first tile. + a = act_fpc_sum[idx] + if a > tol: + pot_fpc[:, idx] = act_fpc[:, idx] / a + else: + pot_fpc[:, idx] = 0.0 + pot_fpc[0, idx] = 1.0 + # Recompute and report + pot_sum2 = np.sum(pot_fpc, axis=0) + bad2 = int(np.sum(np.abs(pot_sum2 - 1.0) > tol)) + print(f" After fix: {bad2} landpoints still failing pot_fpc closure") + if not dry_run: + veg_nc.variables['pot_fpc'][:] = pot_fpc + + # Calculate total FPC + total_fpc = desert_fpc + bare_fpc + act_fpc_sum + tol = 1e-10 + n_problems_before = int(np.sum(np.abs(total_fpc - 1.0) > tol)) + print(f"\nBefore fix: {n_problems_before} landpoints with |total FPC - 1| > {tol:g}") + print(f" Total FPC: min={total_fpc.min():.12f}, max={total_fpc.max():.12f}") + + # === Fix 1: Deglaciated areas === + if n_deglaciated_2d > 0: + deglaciated_landpoints = get_landpoint_indices_from_2d_mask( + deglaciated_mask, landseamask + ) + n_deglaciated = len(deglaciated_landpoints) + print(f"\n--- Fix 1: Deglaciated areas ---") + print(f"Number of deglaciated landpoints: {n_deglaciated}") + + if n_deglaciated > 0: + # For deglaciated areas: + # - desert_fpc = 0 (allow vegetation to grow) + # - Initialize act_fpc to a small value (1% of pot_fpc) if currently 0 + # This is CRITICAL: JSBACH's fpc_to_cover_fract_pot function fails + # when act_fpc=0 for non-glacier points, producing 1e-10 values + # - bare_fpc = 1.0 - act_fpc_sum (so that total = 1) + + # First, initialize act_fpc for points where it's currently 0 + small_init = 0.01 # 1% initial vegetation + n_act_init = 0 + for idx in deglaciated_landpoints: + if act_fpc_sum[idx] < 1e-8 and pot_fpc is not None: + # Initialize act_fpc to a small fraction of pot_fpc + act_fpc[:, idx] = pot_fpc[:, idx] * small_init + n_act_init += 1 + + # Recalculate act_fpc_sum after initialization + act_fpc_sum = np.sum(act_fpc, axis=0) + + if n_act_init > 0: + print(f" Initialized act_fpc for {n_act_init} deglaciated points (1% of pot_fpc)") + if not dry_run: + veg_nc.variables['act_fpc'][:] = act_fpc + + new_desert_fpc = np.zeros(n_deglaciated) + new_bare_fpc = 1.0 - act_fpc_sum[deglaciated_landpoints] + # Do NOT force a positive floor (like 1e-10): it can make total FPC > 1 + new_bare_fpc = np.clip(new_bare_fpc, 0.0, 1.0) + + desert_fpc[deglaciated_landpoints] = new_desert_fpc + bare_fpc[deglaciated_landpoints] = new_bare_fpc + print(f" Set desert_fpc=0 and bare_fpc=1-act_fpc for {n_deglaciated} points") + + # === Fix 2: All remaining FPC issues === + print(f"\n--- Fix 2: All remaining FPC issues ---") + desert_fpc, bare_fpc, n_fixed = fix_all_fpc_issues( + desert_fpc, bare_fpc, act_fpc_sum, tol=tol + ) + print(f" Fixed {n_fixed} additional landpoints") + + # Verify + total_fpc_after = desert_fpc + bare_fpc + act_fpc_sum + n_problems_after = int(np.sum(np.abs(total_fpc_after - 1.0) > tol)) + print(f"\nAfter fix: {n_problems_after} landpoints with |total FPC - 1| > {tol:g}") + print(f" Total FPC: min={total_fpc_after.min():.12f}, max={total_fpc_after.max():.12f}") + + if dry_run: + print("\n[DRY RUN] No changes made to file.") + else: + print("\nApplying changes...") + veg_nc.variables['desert_fpc'][:] = desert_fpc + veg_nc.variables['bare_fpc'][:] = bare_fpc + print("Changes applied successfully!") + + print("\n=== Done ===") + + +if __name__ == '__main__': + main() diff --git a/couplings/utils/fix_jsbach_restart_cover_fract.py b/couplings/utils/fix_jsbach_restart_cover_fract.py new file mode 100644 index 000000000..89abcf705 --- /dev/null +++ b/couplings/utils/fix_jsbach_restart_cover_fract.py @@ -0,0 +1,130 @@ +#!/usr/bin/env python3 +""" +Fix cover_fract and cover_fract_pot in jsbach restart file for deglaciated areas. + +IMPORTANT: This script only updates points where the GLACIER MASK (glac) changed +from 1 to 0. It does NOT touch points that were already non-glacier (like tropical +evergreen forest which also uses Tile 1). + +Usage: + python fix_jsbach_restart_cover_fract.py restart_jsbach.nc jsbach.nc jsbach_reference.nc +""" + +import sys +import numpy as np +import netCDF4 as nc +import shutil +from datetime import datetime + +def fix_cover_fract(restart_file, jsbach_new_file, jsbach_ref_file=None, dry_run=False): + """ + Update cover_fract and cover_fract_pot in restart file for deglaciated areas only. + + Deglaciated areas are identified by comparing glac in jsbach_new with jsbach_reference: + - glac_ref > 0.5 AND glac_new < 0.5 => deglaciated + """ + print(f"=== Fix cover_fract in jsbach restart file ===") + print(f"Restart file: {restart_file}") + print(f"New JSBACH file: {jsbach_new_file}") + print(f"Reference JSBACH file: {jsbach_ref_file}") + print(f"Dry run: {dry_run}") + print() + + # Open new jsbach file + jsbach_new = nc.Dataset(jsbach_new_file, 'r') + cover_fract_new = jsbach_new.variables['cover_fract'][:] + glac_new = jsbach_new.variables['glac'][:] + slm = jsbach_new.variables['slm'][:] + jsbach_new.close() + + # Open reference jsbach file (to identify deglaciated areas) + if jsbach_ref_file: + jsbach_ref = nc.Dataset(jsbach_ref_file, 'r') + glac_ref = jsbach_ref.variables['glac'][:] + jsbach_ref.close() + else: + # If no reference, assume all glac=1 in restart were originally glaciated + glac_ref = np.ones_like(glac_new) + + # Identify truly deglaciated areas (was glacier, now not glacier) + deglaciated_2d = (glac_ref > 0.5) & (glac_new < 0.5) + n_deglaciated_2d = np.sum(deglaciated_2d) + print(f"Deglaciated grid cells (2D): {n_deglaciated_2d}") + + if n_deglaciated_2d == 0: + print("No deglaciated areas found. Nothing to do.") + return 0 + + # Open restart file + if not dry_run: + backup_file = f"{restart_file}.backup_{datetime.now().strftime('%Y%m%d%H%M%S')}" + shutil.copy2(restart_file, backup_file) + print(f"Backup created: {backup_file}") + + restart = nc.Dataset(restart_file, 'r+' if not dry_run else 'r') + + cover_fract_old = restart.variables['cover_fract'][:] + cover_fract_pot_old = restart.variables['cover_fract_pot'][:] + + # Build land mask mapping (2D -> 1D) + land_mask_2d = slm > 0.5 + n_land = int(np.sum(land_mask_2d)) + + print(f"Number of land points: {n_land}") + print(f"Restart cover_fract shape: {cover_fract_old.shape}") + print(f"JSBACH cover_fract shape: {cover_fract_new.shape}") + + # Map 2D to 1D + land_idx_map = {} + land_idx_1d = 0 + for i in range(glac_new.shape[0]): + for j in range(glac_new.shape[1]): + if land_mask_2d[i, j]: + land_idx_map[(i, j)] = land_idx_1d + land_idx_1d += 1 + + # Update only deglaciated points + updated_count = 0 + + for i in range(glac_new.shape[0]): + for j in range(glac_new.shape[1]): + if not land_mask_2d[i, j]: + continue + + if not deglaciated_2d[i, j]: + continue + + land_idx = land_idx_map[(i, j)] + + # Get new cover_fract from jsbach.nc + new_fract = cover_fract_new[:, i, j] + + if not dry_run: + # Update cover_fract + restart.variables['cover_fract'][:, land_idx] = new_fract + # Update cover_fract_pot (same as cover_fract for natural vegetation) + restart.variables['cover_fract_pot'][:, land_idx] = new_fract + + updated_count += 1 + + print(f"\nDeglaciated points updated: {updated_count}") + + if not dry_run: + restart.sync() + + restart.close() + + print("\nDone!") + return updated_count + +if __name__ == "__main__": + if len(sys.argv) < 3: + print("Usage: python fix_jsbach_restart_cover_fract.py restart_jsbach.nc jsbach_new.nc [jsbach_reference.nc] [--dry-run]") + sys.exit(1) + + restart_file = sys.argv[1] + jsbach_new_file = sys.argv[2] + jsbach_ref_file = sys.argv[3] if len(sys.argv) > 3 and not sys.argv[3].startswith('--') else None + dry_run = "--dry-run" in sys.argv + + fix_cover_fract(restart_file, jsbach_new_file, jsbach_ref_file, dry_run) diff --git a/couplings/utils/fix_jsbach_tiny_values.py b/couplings/utils/fix_jsbach_tiny_values.py new file mode 100644 index 000000000..31fa1164a --- /dev/null +++ b/couplings/utils/fix_jsbach_tiny_values.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 +""" +Fix tiny values (1e-10) in jsbach.nc and jsbach restart files. + +The jsbach_init_file program generates values of 1e-10 for tiles/variables +that should have 0 cover. While the sum may be 1.0, these tiny values can +cause JSBACH's validation to fail. + +This script: +1. Sets all values < 1e-8 to exactly 0.0 for fraction-like variables +2. Renormalizes cover_fract to ensure sum = 1.0 on land + +Usage: + python fix_jsbach_tiny_values.py [--restart] + + --restart: Also fix restart-specific variables (cover_fract_pot, veg_ratio, etc.) +""" + +import sys +import numpy as np + +try: + from netCDF4 import Dataset +except ImportError: + print("ERROR: netCDF4 module not found") + sys.exit(1) + + +def fix_variable(nc, var_name, normalize=False): + """Fix tiny values in a variable, optionally normalize to sum=1""" + if var_name not in nc.variables: + return 0 + + var = nc.variables[var_name][:] + tiny_mask = (var > 0) & (var < 1e-8) + n_tiny = np.sum(tiny_mask) + + if n_tiny == 0: + return 0 + + # Set tiny values to 0 + var[tiny_mask] = 0.0 + + if normalize: + # Renormalize so sum = 1.0 along first axis (tiles) + var_sum = var.sum(axis=0, keepdims=True) + # Only normalize where sum > 0 (land points) + var_sum_safe = np.where(var_sum > 0, var_sum, 1.0) + var = var / var_sum_safe + + nc.variables[var_name][:] = var + return n_tiny + + +def main(): + if len(sys.argv) < 2: + print(__doc__) + sys.exit(1) + + jsbach_file = sys.argv[1] + is_restart = '--restart' in sys.argv + + print(f"=== Fix tiny values in {jsbach_file} ===") + print(f"Mode: {'restart' if is_restart else 'init file'}") + + # Variables to fix + if is_restart: + # jsbach restart file variables + vars_to_fix = { + 'cover_fract': True, # normalize + 'cover_fract_pot': True, # normalize + 'forest_fract': False, + 'veg_ratio': False, + 'veg_ratio_max': False, + } + else: + # jsbach init file (jsbach.nc) variables + vars_to_fix = { + 'cover_fract': True, # normalize + 'natural_veg': False, + 'veg_fract': False, + } + + total_fixed = 0 + + with Dataset(jsbach_file, 'r+') as nc: + for var_name, normalize in vars_to_fix.items(): + n_fixed = fix_variable(nc, var_name, normalize) + if n_fixed > 0: + print(f" {var_name}: fixed {n_fixed} tiny values") + total_fixed += n_fixed + + if total_fixed == 0: + print("No tiny values found. Nothing to do.") + else: + print(f"Total: fixed {total_fixed} tiny values") + + print("=== Done ===") + + +if __name__ == '__main__': + main() diff --git a/couplings/utils/fix_veg_restart_act_fpc.py b/couplings/utils/fix_veg_restart_act_fpc.py new file mode 100755 index 000000000..60817f04d --- /dev/null +++ b/couplings/utils/fix_veg_restart_act_fpc.py @@ -0,0 +1,94 @@ +#!/usr/bin/env python3 +""" +Fix act_fpc in veg restart file for deglaciated areas. + +When ice melts (deglaciation), act_fpc at those grid points remains 0 because +the vegetation hasn't had time to grow. This causes fpc_to_cover_fract_pot +to produce near-zero cover_fract_pot values, triggering JSBACH errors. + +Solution: Initialize act_fpc at deglaciated points using pot_fpc values. + +Usage: + python fix_veg_restart_act_fpc.py veg_restart.nc +""" + +import sys +import numpy as np +import netCDF4 as nc +import shutil +from datetime import datetime + +def fix_act_fpc(veg_file, dry_run=False): + """ + Fix act_fpc in veg restart file by initializing deglaciated points + """ + print(f"=== Fix act_fpc in veg restart file ===") + print(f"Veg restart file: {veg_file}") + print(f"Dry run: {dry_run}") + print() + + # Open veg restart + if not dry_run: + backup_file = f"{veg_file}.backup_actfpc_{datetime.now().strftime('%Y%m%d%H%M%S')}" + shutil.copy2(veg_file, backup_file) + print(f"Backup created: {backup_file}") + + veg = nc.Dataset(veg_file, 'r+' if not dry_run else 'r') + + act_fpc = veg.variables['act_fpc'][:] + pot_fpc = veg.variables['pot_fpc'][:] + + print(f"act_fpc shape: {act_fpc.shape}") + print(f"pot_fpc shape: {pot_fpc.shape}") + + # Calculate sums + act_sum = np.sum(act_fpc, axis=0) + pot_sum = np.sum(pot_fpc, axis=0) + + print(f"\nBefore fix:") + print(f" act_fpc sum: min={act_sum.min():.6f}, max={act_sum.max():.6f}") + print(f" pot_fpc sum: min={pot_sum.min():.6f}, max={pot_sum.max():.6f}") + + # Find points where act_fpc sum is very small but pot_fpc sum is near 1 + # These are deglaciated points that need initialization + need_fix = (act_sum < 0.5) & (pot_sum > 0.5) + n_fix = np.sum(need_fix) + + print(f"\nPoints needing fix (act_fpc_sum < 0.5 and pot_fpc_sum > 0.5): {n_fix}") + + if n_fix > 0: + if not dry_run: + # Initialize act_fpc using pot_fpc values for deglaciated points + for i in np.where(need_fix)[0]: + act_fpc[:, i] = pot_fpc[:, i] + + veg.variables['act_fpc'][:] = act_fpc + veg.sync() + print(f"act_fpc updated at {n_fix} points") + else: + print(f"[DRY RUN] Would update act_fpc at {n_fix} points") + + # Verify + if not dry_run and n_fix > 0: + act_fpc_new = veg.variables['act_fpc'][:] + act_sum_new = np.sum(act_fpc_new, axis=0) + still_bad = np.sum(act_sum_new < 0.5) + + print(f"\nAfter fix:") + print(f" act_fpc sum: min={act_sum_new.min():.6f}, max={act_sum_new.max():.6f}") + print(f" Points still with act_fpc sum < 0.5: {still_bad}") + + veg.close() + + print("\nDone!") + return n_fix + +if __name__ == "__main__": + if len(sys.argv) < 2: + print("Usage: python fix_veg_restart_act_fpc.py veg_restart.nc [--dry-run]") + sys.exit(1) + + veg_file = sys.argv[1] + dry_run = "--dry-run" in sys.argv + + fix_act_fpc(veg_file, dry_run) diff --git a/couplings/utils/fix_veg_restart_pot_fpc.py b/couplings/utils/fix_veg_restart_pot_fpc.py new file mode 100644 index 000000000..8f3dcd002 --- /dev/null +++ b/couplings/utils/fix_veg_restart_pot_fpc.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python3 +""" +Fix pot_fpc in veg restart file for deglaciated areas. + +When cover_fract in jsbach restart is updated, pot_fpc in veg restart +must also be updated to match. Otherwise, JSBACH will crash with +memory errors or inconsistent state. + +Usage: + python fix_veg_restart_pot_fpc.py veg_restart.nc jsbach_restart.nc +""" + +import sys +import numpy as np +import netCDF4 as nc +import shutil +from datetime import datetime + +def fix_pot_fpc(veg_file, jsbach_file, dry_run=False): + """ + Update pot_fpc in veg restart file to match cover_fract in jsbach restart + """ + print(f"=== Fix pot_fpc in veg restart file ===") + print(f"Veg restart file: {veg_file}") + print(f"JSBACH restart file: {jsbach_file}") + print(f"Dry run: {dry_run}") + print() + + # Open jsbach restart to get cover_fract + jsbach = nc.Dataset(jsbach_file, 'r') + cover_fract = jsbach.variables['cover_fract'][:] + jsbach.close() + + # Open veg restart + if not dry_run: + backup_file = f"{veg_file}.backup_{datetime.now().strftime('%Y%m%d%H%M%S')}" + shutil.copy2(veg_file, backup_file) + print(f"Backup created: {backup_file}") + + veg = nc.Dataset(veg_file, 'r+' if not dry_run else 'r') + + pot_fpc_old = veg.variables['pot_fpc'][:] + + print(f"pot_fpc shape: {pot_fpc_old.shape}") + print(f"cover_fract shape: {cover_fract.shape}") + + if pot_fpc_old.shape != cover_fract.shape: + print("ERROR: Shape mismatch! Cannot proceed.") + veg.close() + return -1 + + # FIX 2025-01-25: Improved mismatch detection + # Find ANY points where pot_fpc and cover_fract differ significantly + diff = np.abs(pot_fpc_old - cover_fract) + mismatch_mask = np.any(diff > 1e-6, axis=0) # Check across all tiles + n_mismatch = np.sum(mismatch_mask) + + # Also specifically check for the common case: glacier in pot_fpc but not in cover_fract + glacier_mismatch = (pot_fpc_old[0] > 0.5) & (cover_fract[0] < 0.5) + n_glacier_mismatch = np.sum(glacier_mismatch) + + print(f"\nTotal mismatch points (any tile differs > 1e-6): {n_mismatch}") + print(f"Glacier mismatch points (pot_fpc[0]>0.5 but cover_fract[0]<0.5): {n_glacier_mismatch}") + + # Always sync pot_fpc to cover_fract to ensure consistency + # This is safe because cover_fract is the authoritative source from jsbach.nc + if not dry_run: + print(f"\nSyncing pot_fpc to match cover_fract...") + veg.variables['pot_fpc'][:] = cover_fract + veg.sync() + print(f"pot_fpc updated to match cover_fract") + + # Verify + if not dry_run: + pot_fpc_new = veg.variables['pot_fpc'][:] + diff_new = np.abs(pot_fpc_new - cover_fract) + new_mismatch = np.sum(np.any(diff_new > 1e-6, axis=0)) + print(f"\nAfter fix, mismatch points: {new_mismatch}") + + # Verify sum = 1 + pot_sum = np.sum(pot_fpc_new, axis=0) + bad_sum = np.sum(np.abs(pot_sum - 1.0) > 1e-6) + print(f"Points with pot_fpc sum != 1: {bad_sum}") + + veg.close() + + print("\nDone!") + return n_mismatch + +if __name__ == "__main__": + if len(sys.argv) < 3: + print("Usage: python fix_veg_restart_pot_fpc.py veg_restart.nc jsbach_restart.nc [--dry-run]") + sys.exit(1) + + veg_file = sys.argv[1] + jsbach_file = sys.argv[2] + dry_run = "--dry-run" in sys.argv + + fix_pot_fpc(veg_file, jsbach_file, dry_run)