diff --git a/ADASMELIAE.py b/ADASMELIAE.py
index b028941a772b77421473f0facb5999e30c951c98..ba3fbbee963c398971e541a73a5360ff2837d69d 100755
--- a/ADASMELIAE.py
+++ b/ADASMELIAE.py
@@ -45,16 +45,16 @@ DEBUG = (
 model_id = os.getenv("MODEL_ID")
 home_dir = os.getenv("HOME_DIR")
 weather_data_dir = os.getenv("WEATHER_DATA_DIR")
-result_tif_dir = os.getenv("RESULT_TIF_DIR")
-result_mapfile_dir = os.getenv("RESULT_MAPFILE_DIR")
+result_tif_base_dir = os.getenv("RESULT_TIF_DIR")
+result_mapfile_base_dir = os.getenv("RESULT_MAPFILE_DIR")
 tmp_dir = os.getenv("TMP_DIR")
 template_dir = f"{home_dir}mapfile/"
 
+
 weather_data_filename_pattern = os.getenv("WEATHER_DATA_FILENAME_DATEFORMAT")
 start_MM_DD = os.getenv("START_DATE_MM_DD")
 end_MM_DD = os.getenv("END_DATE_MM_DD")
 
-
 config_file_path = f"{home_dir}ADASMELIAE.cfg"
 mask_file_path = (
     f"{home_dir}{os.getenv('MASK_FILE')}" if os.getenv("MASK_FILE") else None
@@ -126,12 +126,12 @@ def remove_temporary_files():
 
 
 # Remove previously calculated results
-def remove_old_results():
+def remove_old_results(tif_dir, mapfile_dir):
     logging.info("Remove previously calculated results")
-    if glob.glob(f"{result_tif_dir}*.tif"):
-        run_command(command=f"rm {result_tif_dir}*.tif")
-    if glob.glob(f"rm {result_mapfile_dir}*.map"):
-        run_command(command=f"rm {result_mapfile_dir}*.map")
+    if glob.glob(f"{tif_dir}*.tif"):
+        run_command(command=f"rm {tif_dir}*.tif")
+    if glob.glob(f"rm {mapfile_dir}*.map"):
+        run_command(command=f"rm {mapfile_dir}*.map")
 
 # Get unique dates in given file
 def get_unique_dates(input_file):
@@ -158,19 +158,39 @@ if __name__ == "__main__":
     start_time = time.time()  # For logging execution time
     today = datetime.now().date()
 
-    year = today.year
+    # Check if a year argument is provided
+    if len(sys.argv) > 1:
+        year = int(sys.argv[1])
+    else:
+        today = datetime.now().date()
+        year = today.year
+
+    # Ensure that model is not run for future year
+    if year > today.year:
+        logging.error(f"Cannot run model for future year ({year}). Quit.")
+        sys.exit()
+
+    # Ensure that output year dirs exist
+    result_tif_year_dir = f"{result_tif_base_dir}/{year}"
+    result_mapfile_year_dir = f"{result_mapfile_base_dir}/{year}"
+    os.makedirs(result_tif_year_dir, exist_ok=True)
+    os.makedirs(result_mapfile_year_dir, exist_ok=True)
+
     model_start_date = datetime.strptime(f"{year}-{start_MM_DD}", "%Y-%m-%d").date()
     model_end_date = datetime.strptime(f"{year}-{end_MM_DD}", "%Y-%m-%d").date()
 
     start_date = model_start_date
-    end_date = today + timedelta(days=2)
+    if year == today.year: # Run for current year
+        end_date = today + timedelta(days=2)
+    else: # Run for previous year
+        end_date = model_end_date + timedelta(days=2)
 
     logging.info(
         f"Model start date {model_start_date}. Model end date {model_end_date}."
     )
 
     logging.info(
-        f"Attempt to run model {model_id} for start date {start_date} and end date {end_date} (today + 2 days)"
+        f"Attempt to run model {model_id} for year {year}, start date {start_date} and end date {end_date}"
     )
 
     # Ensure that model is not run after model run period (include 2 extra days at the end to replace forecast with reanalysed)
@@ -183,7 +203,7 @@ if __name__ == "__main__":
         logging.error("Model period not started. Quit.")
         sys.exit()
     
-    weather_data_file = weather_data_dir + today.strftime(weather_data_filename_pattern)
+    weather_data_file = weather_data_dir + weather_data_filename_pattern.replace('%Y', str(year))
     max_temp_in_period_file = f"{tmp_dir}{year}_max_temp_for_period.nc"
     unmasked_result_file = f"{tmp_dir}{year}_unmasked_result.nc"
     final_result_file = f"{tmp_dir}{year}_result.nc"
@@ -233,7 +253,7 @@ if __name__ == "__main__":
     ]
     logging.info(f"Split into {len(timestep_filenames)} timestep files")
 
-    remove_old_results()
+    remove_old_results(result_tif_year_dir, result_mapfile_year_dir)
 
     for file in timestep_files:
         filename = os.path.basename(file)
@@ -242,8 +262,8 @@ if __name__ == "__main__":
         with open("/dev/null", "w") as devnull:
             warning_status_lcc = f"{tmp_dir}result_WARNING_STATUS_{current_file_date_str}_lcc.tif"
             result_lcc = f"{tmp_dir}result_{current_file_date_str}_lcc.tif"
-            warning_status = f"{result_tif_dir}result_WARNING_STATUS_{current_file_date_str}.tif"
-            result = f"{result_tif_dir}result_{current_file_date_str}.tif"
+            warning_status = f"{result_tif_year_dir}/result_WARNING_STATUS_{current_file_date_str}.tif"
+            result = f"{result_tif_year_dir}/result_{current_file_date_str}.tif"
 
             logging.debug(f"Calculate result for {current_file_date_str}")
 
@@ -265,6 +285,7 @@ if __name__ == "__main__":
                 command=f"gdal_translate -ot Float32 -of GTiff  NETCDF:'{file}':{tm_max} {result_lcc}",
                 stdout=devnull,
             )
+
             # Need to reproject the files, to ensure we have the projection given in the generated mapfile. We always use EPSG:4326 for this
             run_command(
                 command=f"gdalwarp -t_srs EPSG:4326 {warning_status_lcc} {warning_status}",
@@ -305,7 +326,7 @@ if __name__ == "__main__":
             "language_codes": language_codes,
         }
     )
-    with open(f"{result_mapfile_dir}/{model_id}.map", "w") as f:
+    with open(f"{result_mapfile_year_dir}/{model_id}.map", "w") as f:
         f.write(output)
     # remove_temporary_files()
 
diff --git a/env-sample b/env-sample
index a86993f6a23ce435004129246c9fa35f49655e3c..8e246512cb4992424652ec12f71c4186b123f550 100644
--- a/env-sample
+++ b/env-sample
@@ -2,6 +2,8 @@
 
 # Default value is False
 DEBUG=False
+# If not given, the model is run for current year
+YEAR=2024
 # This is used to auto generate some variables and file names
 MODEL_ID="ADASMELIAE"
 # The start date MM-DD of the model
@@ -31,4 +33,4 @@ MAPSERVER_IMAGE_PATH=/foo/mapserver/image/
 # The value of the EXTENT parameter in Mapserver's mapfile. Units are DD (Decimal degrees)
 MAPSERVER_EXTENT="-23.5 29.5 62.5 70.5"
 # Path to optional CSV file with polygons for masking result.
-MASK_FILE=europe_coastline.csv
\ No newline at end of file
+MASK_FILE=europe_coastline.csv
diff --git a/run_ADASMELIAE.sh b/run_ADASMELIAE.sh
index caa9d0c309a7934b11c6fa6f33822b48841947cf..4e9547c55145a88230543660e49f8d88d89fa0a7 100755
--- a/run_ADASMELIAE.sh
+++ b/run_ADASMELIAE.sh
@@ -42,6 +42,12 @@ then
     exit
 fi
 
+if [ -z "${YEAR}" ]
+then
+    YEAR=$(date +%Y)
+    echo "YEAR is not set. Using the current year: ${YEAR}"
+fi
+
 # Paths to scripts and requirements
 LOG_FILE=${HOME_DIR}log/ADASMELIAE.log
 REQUIREMENTS=${HOME_DIR}requirements.txt
@@ -60,7 +66,7 @@ if [ -n "$VIRTUAL_ENV" ]; then
     pip install -q -r $REQUIREMENTS
 
     # Run the model
-    python3 ${HOME_DIR}ADASMELIAE.py >> "$LOG_FILE" 2>&1
+    python3 ${HOME_DIR}ADASMELIAE.py "$YEAR" >> "$LOG_FILE" 2>&1
 
     # Deactivate the virtual environment
     deactivate