PBS job arrays let you submit many similar jobs with one script, where each task runs with a different index. This is ideal for looping over forecast hours, ensemble members, dates, or domains.

1. Basic PBS array syntax
#PBS -N gfs_apcp
#PBS -l select=1:ncpus=4:mem=16gb
#PBS -l walltime=01:00:00
#PBS -J 1-40          # array indices 1..40
#PBS -j oe
#PBS -r y             # MUST be rerunnable for arrays

cd $PBS_O_WORKDIR
echo "This is array task: $PBS_ARRAY_INDEX"


Submit with:

qsub gfs_apcp.pbs


Each task runs independently with a different
$PBS_ARRAY_INDEX.

2. Passing array index into Python
PBS script
python plot_precip.py ${PBS_ARRAY_INDEX}

Python
import sys

idx = int(sys.argv[1])     # 1..40
fhr = idx * 6              # example: 6-hourly forecast
print("Forecast hour:", fhr)

3. Mapping array index to real values
Example: forecast hours
FHR_LIST=(006 012 018 024 030 036)
IDX=$((PBS_ARRAY_INDEX-1))
FHR=${FHR_LIST[$IDX]}
python plot.py $FHR

Example: dates
START=20240101
DATE=$(date -d "$START +$((PBS_ARRAY_INDEX-1)) days" +%Y%m%d)
python run.py $DATE

4. Output files per array task

Use the index to avoid overwriting:

#PBS -o logs/job_${PBS_ARRAY_INDEX}.out
#PBS -e logs/job_${PBS_ARRAY_INDEX}.err

5. Common errors (you’ve likely seen this one)
❌ cannot submit non-rerunable Array Job

Cause: Arrays must be rerunnable.

Fix:

#PBS -r y


(or remove -r n if present)

6. Throttling array concurrency

Run only N tasks at once:

#PBS -J 1-100%10   # max 10 running simultaneously


Very useful on shared NOAA / HPC systems.

7. Environment variables available in arrays
Variable	Meaning
PBS_ARRAY_INDEX	Current task number
PBS_ARRAY_ID	Job array ID
PBS_JOBID	Unique job ID
PBS_O_WORKDIR	Submission directory
8. Typical NOAA-style use cases

GFS: one array index per forecast hour

GEFS: one index per ensemble member

NMME: one index per model or lead

Climatology: one index per month/year

If you want, I can:

Convert one of your existing PBS scripts to an array

Design a GFS/GEFS array layout (fhrs × domains)

Debug a specific PBS array error you’re hitting

yes Design a GFS array layout (fhrs × domains)

Below is a clean, production-ready GFS PBS array design that scales over forecast hours × domains, matches NOAA/HPC practices, and avoids job explosion.

GFS PBS Array Layout (Forecast Hours × Domains)
1. Design strategy (important)

Instead of nested arrays (PBS does not support them), we:

✅ Use ONE array index
✅ Map index → (forecast hour, domain)
✅ Throttle concurrency
✅ Keep memory/CPU predictable

2. Define the problem space
Example

Forecast hours: 0–240 every 6 h
→ 41 times: 000 006 ... 240

Domains:

global
africa
conus
asia

Total array size
41 fhrs × 4 domains = 164 tasks

3. PBS array script (recommended layout)
gfs_precip_array.pbs
#!/bin/bash
#PBS -N GFS_APCP
#PBS -l select=1:ncpus=4:mem=16gb
#PBS -l walltime=01:00:00
#PBS -J 1-164%12          # 12 tasks running at once
#PBS -j oe
#PBS -r y

cd $PBS_O_WORKDIR

# ----------------------------
# Configuration
# ----------------------------
FHR_LIST=( $(seq -f "%03g" 0 6 240) )
DOMAINS=(global africa conus asia)

NFHR=${#FHR_LIST[@]}
NDOM=${#DOMAINS[@]}

IDX=$((PBS_ARRAY_INDEX-1))

# ----------------------------
# Index mapping
# ----------------------------
FHR_INDEX=$(( IDX / NDOM ))
DOM_INDEX=$(( IDX % NDOM ))

FHR=${FHR_LIST[$FHR_INDEX]}
DOMAIN=${DOMAINS[$DOM_INDEX]}

echo "Array index : $PBS_ARRAY_INDEX"
echo "Forecast hr : $FHR"
echo "Domain      : $DOMAIN"

# ----------------------------
# Run plotting
# ----------------------------
python gfs_maps_precip_6hrly.py \
       --fhr $FHR \
       --domain $DOMAIN

4. Python side (clean & explicit)
gfs_maps_precip_6hrly.py
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("--fhr", required=True)
parser.add_argument("--domain", required=True)
args = parser.parse_args()

fhr = args.fhr
domain = args.domain

print(f"Processing F{fhr} for domain {domain}")

5. Why this layout is optimal

✅ Single submission

qsub gfs_precip_array.pbs


✅ No duplicate code
✅ Easy to add domains or hours
✅ Scheduler-friendly
✅ Natural log separation by index

6. Output & logging (best practice)
#PBS -o logs/gfs_${PBS_ARRAY_INDEX}.out
#PBS -e logs/gfs_${PBS_ARRAY_INDEX}.err


Inside Python:

outfile = f"gfs_apcp_{domain}_f{fhr}.png"

7. Alternative layout (domain-major order)

If you prefer all fhrs of one domain run together:

DOM_INDEX=$(( IDX / NFHR ))
FHR_INDEX=$(( IDX % NFHR ))


Use this if:

Domain-specific data cached locally

You want better I/O locality

8. Scaling rules (NOAA systems)
Situation	Recommendation
Heavy GRIB I/O	%8–12 throttle
Cartopy plots	2–4 CPUs
Global plots	≥16 GB RAM
Many domains	Domain-major ordering
9. Common failure points (and fixes)
❌ Only one FHR runs

✔ Forgot #PBS -J

❌ non-rerunable Array Job

✔ Remove -r n or add #PBS -r y

❌ Index out of range

✔ Check:

#PBS -J 1-$((NFHR*NDOM))

10. Want next?

I can:

Convert this to GEFS (member × fhr × domain)

Add automatic PDY/cycle handling

Build restart-safe logic (skip if PNG exists)

Optimize for WCOSS / Hera / Jet

Just tell me which system you’re on and what domains you want.

Build restart-safe logic (skip if PNG exists)

Below is a restart-safe, production-grade pattern that NOAA ops commonly use:
if the PNG already exists → the array task exits cleanly.
This lets you re-qsub the same PBS array after a crash with zero wasted CPU.

Restart-Safe PBS Array Logic (Skip if PNG Exists)
1. Principle

Each array task is idempotent:

Same (fhr, domain) → same output filename
If file exists and is valid → skip

2. PBS script (restart-safe)
gfs_precip_array.pbs
#!/bin/bash
#PBS -N GFS_APCP
#PBS -l select=1:ncpus=4:mem=16gb
#PBS -l walltime=01:00:00
#PBS -J 1-164%12
#PBS -j oe
#PBS -r y

cd $PBS_O_WORKDIR || exit 1

# ----------------------------
# Configuration
# ----------------------------
FHR_LIST=( $(seq -f "%03g" 0 6 240) )
DOMAINS=(global africa conus asia)

NFHR=${#FHR_LIST[@]}
NDOM=${#DOMAINS[@]}

IDX=$((PBS_ARRAY_INDEX-1))
FHR_INDEX=$(( IDX / NDOM ))
DOM_INDEX=$(( IDX % NDOM ))

FHR=${FHR_LIST[$FHR_INDEX]}
DOMAIN=${DOMAINS[$DOM_INDEX]}

OUTDIR=png/${DOMAIN}
OUTFILE=${OUTDIR}/gfs_apcp_${DOMAIN}_f${FHR}.png

mkdir -p ${OUTDIR}

# ----------------------------
# Restart-safe check
# ----------------------------
if [[ -s ${OUTFILE} ]]; then
    echo "SKIP: ${OUTFILE} already exists"
    exit 0
fi

echo "RUN : F${FHR} ${DOMAIN}"

python gfs_maps_precip_6hrly.py \
       --fhr $FHR \
       --domain $DOMAIN \
       --outfile ${OUTFILE}

STATUS=$?

if [[ $STATUS -ne 0 ]]; then
    echo "FAILED: F${FHR} ${DOMAIN}"
    exit $STATUS
fi

echo "DONE: ${OUTFILE}"

Why -s?
[[ -s file.png ]]


✔ File exists
✔ File size > 0 (guards against crashed zero-byte files)

3. Python side (atomic write – VERY important)

Never write directly to the final PNG.
Write to a temp file, then rename.

gfs_maps_precip_6hrly.py
import argparse
import os
import matplotlib.pyplot as plt

parser = argparse.ArgumentParser()
parser.add_argument("--fhr", required=True)
parser.add_argument("--domain", required=True)
parser.add_argument("--outfile", required=True)
args = parser.parse_args()

tmpfile = args.outfile + ".tmp"

# ---- plotting ----
plt.figure()
plt.title(f"GFS APCP F{args.fhr} {args.domain}")
plt.plot([0, 1], [0, 1])  # placeholder
plt.savefig(tmpfile, dpi=150, bbox_inches="tight")
plt.close()

# ---- atomic rename ----
os.replace(tmpfile, args.outfile)

Why this matters

If the job crashes:

No partial PNG

No false “exists” detection

Safe restarts every time

4. Optional: checksum validation (extra safety)

If you want strong validation:

if [[ -s ${OUTFILE} ]] && python - <<EOF
from PIL import Image
Image.open("${OUTFILE}").verify()
EOF
then
    echo "VALID PNG → SKIP"
    exit 0
fi


(Usually not needed unless filesystem is flaky.)

5. Restart workflow (ops-style)
# First run
qsub gfs_precip_array.pbs

# Node crash / walltime hit / queue issue

# Just resubmit — completed plots will be skipped
qsub gfs_precip_array.pbs


No tracking, no state files, no bookkeeping.

6. Bonus: skip logic inside Python (alternative)

If PBS logic is locked down:

if os.path.exists(outfile) and os.path.getsize(outfile) > 0:
    print("SKIP:", outfile)
    sys.exit(0)


PBS-side skipping is preferred (cheaper).

7. NOAA best-practice summary
Feature	Recommendation
Output naming	Deterministic
Skip test	-s outfile
Write method	temp → atomic rename
Arrays	Rerunnable
Recovery	Re-qsub same job

If you want, I can next:

Add walltime-aware early exit

Build GEFS restart-safe member × fhr × domain

Add “touch .done” sentinel files

Integrate with PDY/cycle auto-discovery

Just say which one.

Integrate with PDY/cycle auto-discovery

Below is a fully integrated, ops-style solution that automatically discovers PDY and cycle, works with restart-safe PBS arrays, and requires zero manual edits per run.

This pattern is widely used for GFS / GEFS on NOAA systems.

GFS PBS Array with PDY & Cycle Auto-Discovery
1. Assumptions (adjust paths if needed)

GFS GRIB location follows:

/path/to/gfs/prod/gfs.${PDY}/${CYC}/


Filenames:

gfs.t${CYC}z.pgrb2.0p25.f006

2. PDY & cycle discovery logic (robust)
Rules

Use latest complete cycle

Prefer 00 → 06 → 12 → 18

Fall back automatically if latest not fully available

3. PBS script (production-ready)
gfs_precip_array.pbs
#!/bin/bash
#PBS -N GFS_APCP
#PBS -l select=1:ncpus=4:mem=16gb
#PBS -l walltime=01:00:00
#PBS -J 1-164%12
#PBS -j oe
#PBS -r y

cd $PBS_O_WORKDIR || exit 1

# ----------------------------
# Discover PDY & cycle
# ----------------------------
GFSROOT=/path/to/gfs/prod

for DAY in $(date -u +%Y%m%d) $(date -u -d "yesterday" +%Y%m%d); do
  for CYC in 18 12 06 00; do
    TESTDIR=${GFSROOT}/gfs.${DAY}/${CYC}
    TESTFILE=${TESTDIR}/gfs.t${CYC}z.pgrb2.0p25.f006
    if [[ -s $TESTFILE ]]; then
        PDY=$DAY
        cycle=$CYC
        GFS_DIR=$TESTDIR
        break 2
    fi
  done
done

if [[ -z $PDY ]]; then
    echo "ERROR: No valid GFS cycle found"
    exit 2
fi

echo "Using GFS: PDY=${PDY} CYC=${cycle}"

# ----------------------------
# Forecast hours & domains
# ----------------------------
FHR_LIST=( $(seq -f "%03g" 0 6 240) )
DOMAINS=(global africa conus asia)

NFHR=${#FHR_LIST[@]}
NDOM=${#DOMAINS[@]}

IDX=$((PBS_ARRAY_INDEX-1))
FHR_INDEX=$(( IDX / NDOM ))
DOM_INDEX=$(( IDX % NDOM ))

FHR=${FHR_LIST[$FHR_INDEX]}
DOMAIN=${DOMAINS[$DOM_INDEX]}

# ----------------------------
# Output
# ----------------------------
OUTDIR=png/${PDY}/${cycle}/${DOMAIN}
OUTFILE=${OUTDIR}/gfs_apcp_${PDY}_${cycle}_f${FHR}.png
mkdir -p ${OUTDIR}

# ----------------------------
# Restart-safe skip
# ----------------------------
if [[ -s ${OUTFILE} ]]; then
    echo "SKIP: ${OUTFILE}"
    exit 0
fi

# ----------------------------
# Run plot
# ----------------------------
python gfs_maps_precip_6hrly.py \
       --pdy ${PDY} \
       --cycle ${cycle} \
       --fhr ${FHR} \
       --domain ${DOMAIN} \
       --gfsdir ${GFS_DIR} \
       --outfile ${OUTFILE}

STATUS=$?
[[ $STATUS -ne 0 ]] && exit $STATUS

echo "DONE: ${OUTFILE}"

4. Python side (clean interface)
gfs_maps_precip_6hrly.py
import argparse
import os
import sys

parser = argparse.ArgumentParser()
parser.add_argument("--pdy", required=True)
parser.add_argument("--cycle", required=True)
parser.add_argument("--fhr", required=True)
parser.add_argument("--domain", required=True)
parser.add_argument("--gfsdir", required=True)
parser.add_argument("--outfile", required=True)
args = parser.parse_args()

grib = f"{args.gfsdir}/gfs.t{args.cycle}z.pgrb2.0p25.f{args.fhr}"

if not os.path.exists(grib):
    print("MISSING:", grib)
    sys.exit(3)

tmp = args.outfile + ".tmp"

# ---- plotting here ----
# plt.savefig(tmp)

os.replace(tmp, args.outfile)

5. Why this is ops-grade

✅ No hard-coded PDY
✅ No hard-coded cycle
✅ Safe re-runs
✅ Handles late data
✅ One submission per day

6. Optional upgrades (easy add-ons)
A. Require full forecast availability
[[ ! -s ${GFS_DIR}/gfs.t${cycle}z.pgrb2.0p25.f240 ]] && continue

B. Environment override (manual reruns)
PDY=${PDY_OVERRIDE:-$PDY}
cycle=${CYCLE_OVERRIDE:-$cycle}

C. Log metadata
echo "${PDY} ${cycle} F${FHR} ${DOMAIN}" >> run.log

7. Typical NOAA directory outcome
png/
 └── 20260111/
     └── 00/
         ├── africa/
         │   ├── gfs_apcp_20260111_00_f006.png
         │   └── ...
         ├── conus/
         └── global/

8. Next logical steps

I can:

Add GEFS auto-discovery (member dirs)

Add cycle-age thresholds (e.g., ≥2h old)

Integrate DBNet / HPSS archiving

Build multi-resolution (0p25 / 0p50) logic

Tell me which HPC (WCOSS2 / Hera / Jet / local) and I’ll tune it exactly.

does the pbs array uses individual cpu or its cores ?

Short answer:

👉 A PBS array uses separate jobs, not shared cores.
Each array element is an independent job that gets its own CPU cores as requested.

How PBS arrays actually work

Think of a PBS array as:

“Submit N independent jobs with the same script but different indices.”

Example
#PBS -J 1-20
#PBS -l select=1:ncpus=4


This creates 20 separate jobs.

Each job gets:

1 node

4 CPU cores

They do NOT share cores with other array elements.

What resources are allocated?
Per array element:
select=1:ncpus=4:mem=16gb


Means each task gets:

4 cores

16 GB RAM

Total possible usage (if all run at once):
20 × 4 cores = 80 cores
20 × 16 GB  = 320 GB RAM


This is why throttling matters.

Throttling controls total core usage
#PBS -J 1-164%12


At most:

12 jobs × 4 cores = 48 cores total


This is the correct way to limit load.

Important distinction: cores vs threads
PBS allocates cores

PBS does not know about threads.

If your Python code uses:

NumPy / MKL

OpenMP

BLAS

They may spawn multiple threads per job, unless controlled.

CRITICAL best practice (you should do this)

Add this to your PBS script:

export OMP_NUM_THREADS=1
export MKL_NUM_THREADS=1
export OPENBLAS_NUM_THREADS=1
export NUMEXPR_NUM_THREADS=1


Otherwise:

Each array job may use all cores on the node

You accidentally oversubscribe CPUs

Jobs run slower and annoy schedulers 😄

Cartopy / matplotlib case (your workload)

For your GFS plotting jobs:

Scenario	Recommended
One plot per task	ncpus=1–2
Heavy Cartopy projections	ncpus=2–4
Pure Python	OMP_NUM_THREADS=1

Most of the time:

#PBS -l select=1:ncpus=2:mem=8gb


is enough.

Common misunderstanding (very common)

❌ “Array jobs split cores among themselves”
✔ False

❌ “Array jobs run on the same node sharing CPUs”
✔ Not guaranteed; usually false

Each element is scheduled like a normal job.

Summary (one-glance)
Question	Answer
Does array share CPUs?	❌ No
Does each element get cores?	✅ Yes
Is %N a core limit?	❌ No (job count limit)
How to limit cores overall?	%N × ncpus
Should threads be capped?	✅ Always

If you want, I can:

Tune optimal ncpus for your Cartopy plots

Show node-sharing vs exclusive (place=shared/excl)

Help you reduce memory footprint per array job

Just say the word.

yes

Great — here is a practical, tested tuning guide specifically for GFS Cartopy / matplotlib plotting in PBS arrays, based on NOAA-style workloads like yours.

Optimal ncpus Tuning for Cartopy Plot Jobs
1. The key truth (important)

👉 Cartopy is NOT CPU-parallel
👉 More cores ≠ faster plots

Most of the time:

Cartopy + matplotlib run single-threaded

Extra CPUs just sit idle unless NumPy/OpenMP misbehaves

2. Recommended configurations (use this table)
A. Typical GFS precipitation maps (your case)
Plot type	Recommended
Single map per job	ncpus=1
Heavy coastlines/borders	ncpus=2
Africa / CONUS regional	ncpus=1
Global 0.25°	ncpus=2
Global 0.1° (rare)	ncpus=3–4

✅ Start with ncpus=1
Increase only if runtime is clearly CPU-bound.

3. What actually slows Cartopy (not CPU)
Bottleneck	Solution
Shapefile loading	Cache features
High-res coastlines	Use 110m or 50m
Too many contour levels	Reduce levels
Large grids	Decimate before plotting
Disk I/O (GRIB)	Local scratch
4. Best-practice PBS resource line (recommended)

For your GFS array:

#PBS -l select=1:ncpus=1:mem=6gb


And throttle:

#PBS -J 1-164%24


This gives:

24 jobs × 1 core = 24 cores total


Much better cluster citizenship 👍

5. CRITICAL: force single-threaded math libs

Put this near the top of your PBS script:

export OMP_NUM_THREADS=1
export MKL_NUM_THREADS=1
export OPENBLAS_NUM_THREADS=1
export NUMEXPR_NUM_THREADS=1


Without this:

Each job may spawn 8–32 threads

Even with ncpus=1

Massive oversubscription

6. Cartopy-specific optimizations (huge gains)
A. Cache Cartopy features (DO THIS)
export CARTOPY_USER_BACKGROUNDS=/scratch/$USER/cartopy
export CARTOPY_CACHE_DIR=/scratch/$USER/cartopy


First run is slower; after that, plots are much faster.

B. Use lower-resolution features
import cartopy.feature as cfeature

ax.coastlines(resolution="50m", linewidth=0.6)
ax.add_feature(cfeature.BORDERS.with_scale("50m"), linewidth=0.4)


Avoid 10m unless absolutely needed.

C. Decimate BEFORE plotting (you already do this 👍)
precip = precip[::2, ::2]
lats   = lats[::2, ::2]
lons   = lons[::2, ::2]


This often halves runtime.

7. When should you use ncpus > 1?

Only if ONE job does multiple plots:

for domain in domains:
    make_plot(domain)


In that case:

ncpus=2–4


But for array = one plot per task → ncpus=1 is ideal.

8. Recommended final PBS template (drop-in)
#PBS -N GFS_APCP
#PBS -l select=1:ncpus=1:mem=6gb
#PBS -l walltime=00:20:00
#PBS -J 1-164%24
#PBS -j oe
#PBS -r y

export OMP_NUM_THREADS=1
export MKL_NUM_THREADS=1
export OPENBLAS_NUM_THREADS=1
export NUMEXPR_NUM_THREADS=1

9. How to verify you chose correctly

Inside Python:

import os
print("OMP_NUM_THREADS =", os.getenv("OMP_NUM_THREADS"))


On compute node:

top


You should see:

~100% of one CPU

Not 400%, 800%, etc.

10. Rule of thumb (remember this)

Arrays scale by job count, not by CPU count per job

For plotting:

Many jobs × 1 core  >>>  Few jobs × many cores

Next optimizations I can do for you

Profile exact runtime per FHR

Build node-local GRIB staging

Optimize Cartopy transform calls

Convert to batch multi-plot per job (even faster)

Just tell me which one you want next. 
