Allow user to configure source coordinate reference system

elmerice/Solvers/yac2elmer.F90

@@ -24,7 +24,7 @@ SUBROUTINE YAC2Elmer( Model,Solver,dt,TransientSimulation )
   ! parameters to be read in from this solvers section in the sif 
   LOGICAL :: couple_to_ebfm, couple_to_icon         ! define which component is coupled to Elmer 
 
-  CHARACTER(LEN=1024) ::  config_file, model_tstep
+  CHARACTER(LEN=1024) ::  config_file, model_tstep, grid_crs
   INTEGER :: I, t, ierr, dt_hours
   INTEGER, POINTER :: t_icePerm(:), smbPerm(:), runoffPerm(:)
   ! INTEGER, POINTER :: cltPerm(:), prPerm(:)  ! ICON is not supported at the moment
@@ -60,6 +60,15 @@ SUBROUTINE YAC2Elmer( Model,Solver,dt,TransientSimulation )
      CALL FATAL(SolverName,'No keyword >Couple To ICON< found in yac2elmer solver')
   END IF
 
+  ! read grid CRS (coordinate reference system)
+  grid_crs = GetString(SolverParams, 'Grid CRS', Found)
+  IF (.NOT. Found) THEN
+     grid_crs = "EPSG:3413"  ! default value
+     CALL INFO(SolverName, 'No keyword >Grid CRS< found, using default: EPSG:3413', Level=3)
+  ELSE
+     CALL INFO(SolverName, 'Using grid projection (CRS): ' // TRIM(grid_crs), Level=3)
+  END IF
+
   IF (.NOT. (couple_to_ebfm .OR. couple_to_icon)) THEN
     CALL FATAL(SolverName,'At least one of >Couple To EBFM< or >Couple To ICON< must be TRUE')
   END IF
@@ -91,7 +100,7 @@ SUBROUTINE YAC2Elmer( Model,Solver,dt,TransientSimulation )
         'Running with ' // I2S(ParEnv % PEs) // ' partitions', Level=3)
     END IF
 
-    CALL coupling_setup(ThisMesh, TRIM(model_tstep), couple_to_ebfm, couple_to_icon)
+    CALL coupling_setup(ThisMesh, TRIM(grid_crs), TRIM(model_tstep), couple_to_ebfm, couple_to_icon)
 
     IF (couple_to_ebfm) THEN
       ! setting up Elmer-side variables for receiving YAC variables

fem/src/elmer_coupling.F90

@@ -729,9 +729,16 @@ SUBROUTINE coupling_init(coupling_config_file, elmer_rank, yac_comm)
 
   END SUBROUTINE coupling_init
 
-  SUBROUTINE coupling_setup(grid, timestepstring, couple_to_ebfm_in, couple_to_icon_in)
+  !> Set up the coupling configuration for YAC
+  !>
+  !> @param grid The Elmer mesh/grid to be coupled
+  !> @param grid_crs Coordinate reference system (CRS) of the grid (e.g., "EPSG:3413")
+  !> @param timestepstring Timestep configuration string for YAC
+  !> @param couple_to_ebfm_in Enable coupling to EBFM
+  !> @param couple_to_icon_in Enable coupling to ICON
+  SUBROUTINE coupling_setup(grid, grid_crs, timestepstring, couple_to_ebfm_in, couple_to_icon_in)
 
-    USE, INTRINSIC :: iso_c_binding, ONLY: C_INT, C_DOUBLE
+    USE, INTRINSIC :: iso_c_binding, ONLY: C_INT, C_DOUBLE, C_CHAR, C_NULL_CHAR
 
     ! need to use Types_ without Messages because of circular dependency
     ! originating from coupling_finalize in Messages Fatal
@@ -743,6 +750,7 @@ SUBROUTINE coupling_setup(grid, timestepstring, couple_to_ebfm_in, couple_to_ico
     INTEGER, POINTER :: this_cell_ids(:)
 
     TYPE(Mesh_t), POINTER, INTENT(IN) :: grid
+    CHARACTER(LEN=*), INTENT(IN) :: grid_crs
     TYPE(Element_t), POINTER :: element
     CHARACTER(LEN=*), INTENT(IN) :: timestepstring
     LOGICAL, INTENT(IN) :: couple_to_ebfm_in, couple_to_icon_in
@@ -764,16 +772,17 @@ SUBROUTINE coupling_setup(grid, timestepstring, couple_to_ebfm_in, couple_to_ico
 
     INTERFACE
 
-      SUBROUTINE convert_epsg3413_to_lonlat_c(x_vertices, y_vertices, nbr_vertices) &
-        bind ( C, name='convert_epsg3413_to_lonlat' )
+      SUBROUTINE convert_to_lonlat_c(x_vertices, y_vertices, nbr_vertices, crs) &
+        bind ( C, name='convert_to_lonlat' )
 
-        USE, INTRINSIC :: iso_c_binding, ONLY: C_INT, C_DOUBLE
+        USE, INTRINSIC :: iso_c_binding, ONLY: C_INT, C_DOUBLE, C_CHAR
 
         INTEGER(KIND=C_INT), VALUE, INTENT(IN)    :: nbr_vertices
         REAL(C_DOUBLE),             INTENT(INOUT) :: x_vertices(*)
         REAL(C_DOUBLE),             INTENT(INOUT) :: y_vertices(*)
+        CHARACTER(KIND=C_CHAR),     INTENT(IN)    :: crs(*)
 
-      END SUBROUTINE convert_epsg3413_to_lonlat_c
+      END SUBROUTINE convert_to_lonlat_c
 
     END INTERFACE
 
@@ -816,8 +825,8 @@ END SUBROUTINE convert_epsg3413_to_lonlat_c
       y_cells(i) = SUM(y_vertices(this_cell_ids(1:n))) / n
     END DO
 
-    CALL convert_epsg3413_to_lonlat_c(x_vertices, y_vertices, nbr_vertices)
-    CALL convert_epsg3413_to_lonlat_c(x_cells, y_cells, nbr_cells)
+    CALL convert_to_lonlat_c(x_vertices, y_vertices, nbr_vertices, TRIM(grid_crs)//C_NULL_CHAR)
+    CALL convert_to_lonlat_c(x_cells, y_cells, nbr_cells, TRIM(grid_crs)//C_NULL_CHAR)
 
     ! register Elmer grid in YAC
     ! * is defined as an unstructured grid

fem/src/project_to_lonlat.c

@@ -37,11 +37,11 @@
 
 #include "project_to_lonlat.h"
 
-void convert_epsg3413_to_lonlat(double * x, double * y, int n) {
+void convert_to_lonlat(double * x, double * y, int n, const char * crs) {
   // define transformation
   PJ * P =
     proj_create_crs_to_crs(
-      PJ_DEFAULT_CTX, "EPSG:3413", "+proj=longlat +datum=WGS84", NULL);
+      PJ_DEFAULT_CTX, crs, "+proj=longlat +datum=WGS84", NULL);
 
   if (!P) {
     fputs("failed to create transformation", stderr);

fem/src/project_to_lonlat.h

@@ -34,12 +34,13 @@
 #define PROJECT_TO_LONLAT_H
 
 /**
- * Convert coordinates from EPSG:3413 to longitude/latitude in radians.
+ * Convert coordinates from a specified coordinate reference system (CRS) to longitude/latitude in radians.
  * 
- * @param x Array of x coordinates (EPSG:3413), modified in place to longitude in radians.
- * @param y Array of y coordinates (EPSG:3413), modified in place to latitude in radians.
+ * @param x Array of x coordinates (in source CRS), modified in place to longitude in radians.
+ * @param y Array of y coordinates (in source CRS), modified in place to latitude in radians.
  * @param n Number of coordinates in the arrays.
+ * @param crs Source CRS (e.g., "EPSG:3413").
  */
-void convert_epsg3413_to_lonlat(double * x, double * y, const int n);
+void convert_to_lonlat(double * x, double * y, const int n, const char * crs);
 
 #endif // PROJECT_TO_LONLAT_H