From 0656885d7cf084920d1166b864eeea8fcc99bc85 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ipsl.fr>
Date: Thu, 12 Sep 2024 15:58:15 +0200
Subject: [PATCH] Manage non-uniform longitude-latitude grid
In the case of non-uniform longitude-latitude grid: we read the whole
coordinate arrays from the NetCDF file; we have a different way of
inverting the projection, in `invert_proj`; the program `eddy_graph`
does not use `corner_whole` and step, which are not defined, it uses
instead the shapefile `extr_proj`.
---
Inst_eddies/Tests/test_set_max_speed.f90 | 3 +-
Inst_eddies/input_ssh.f90 | 65 +++++++++++++++++-------
Inst_eddies/invert_proj.f90 | 27 ++++++++--
Overlap/read_grid.F90 | 18 ++++---
4 files changed, 85 insertions(+), 28 deletions(-)
diff --git a/Inst_eddies/Tests/test_set_max_speed.f90 b/Inst_eddies/Tests/test_set_max_speed.f90
index b3512d21..39c3437b 100644
--- a/Inst_eddies/Tests/test_set_max_speed.f90
+++ b/Inst_eddies/Tests/test_set_max_speed.f90
@@ -79,7 +79,8 @@ program test_set_max_speed
CALL shpc_close(hshpc)
if (e%out_cont%closed) then
- e%extr%coord_proj = nint((e%extr%coord - corner_whole) / step) + 1
+ if (uniform_lon_lat) &
+ e%extr%coord_proj = nint((e%extr%coord - corner_whole) / step) + 1
call get_command_arg_dyn(4, path, &
"Required 4th argument: cont-list-shapefile")
diff --git a/Inst_eddies/input_ssh.f90 b/Inst_eddies/input_ssh.f90
index 8cb6c176..b3d0115b 100644
--- a/Inst_eddies/input_ssh.f90
+++ b/Inst_eddies/input_ssh.f90
@@ -14,6 +14,18 @@ module input_ssh_m
logical, protected:: uniform_lon_lat = .true.
! grid is cartesian and uniform in longitude and latitude
+ real, allocatable, save, protected:: longitude(:) ! (nlon)
+ ! allocated only for non-uniform grid, coordinate for the whole grid, in rad
+
+ real, allocatable, save, protected:: latitude(:) ! (nlat)
+ ! allocated only for non-uniform grid, coordinate for the whole grid, in rad
+
+ real, allocatable, save, protected:: lon_steps(:) ! (nlon - 1)
+ ! allocated only for non-uniform grid, longitude steps, in rad
+
+ real, allocatable, save, protected:: lat_steps(:) ! (nlat - 1)
+ ! allocated only for non-uniform grid, longitude steps, in rad
+
private get_var
contains
@@ -23,9 +35,9 @@ contains
use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
! Libraries:
- use jumble, only: new_unit, assert, deg_to_rad
+ use jumble, only: new_unit, assert, deg_to_rad, ediff1d
use netcdf95, only: nf95_open, nf95_find_coord, nf95_inquire_dimension, &
- nf95_get_var, nf95_close, nf95_nowrite
+ nf95_get_var, nf95_close, nf95_nowrite, nf95_gw_var
use config_m, only: max_radius
@@ -65,30 +77,47 @@ contains
print *, "input_ssh: Enter input_ssh_nml:"
read(unit = *, nml = input_ssh_nml)
call nf95_open(ssh_fname, nf95_nowrite, ncid)
- call nf95_find_coord(ncid, dimid = dimid, varid = varid, &
- std_name = "longitude")
- call nf95_inquire_dimension(ncid, dimid, nclen = nlon)
- call nf95_get_var(ncid, varid, corner_deg(1))
- call nf95_get_var(ncid, varid, lon_max, start = [nlon])
- call nf95_find_coord(ncid, dimid = dimid, varid = varid, &
- std_name = "latitude")
- call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
- call nf95_get_var(ncid, varid, corner_deg(2))
- call nf95_get_var(ncid, varid, lat_max, start = [nlat])
- call assert(lon_max > corner_deg(1) .and. lat_max > corner_deg(2), &
- "input_ssh: coordinate order")
- step_deg = [(lon_max - corner_deg(1)) / (nlon - 1), &
- (lat_max - corner_deg(2)) / (nlat - 1)]
+
+ if (uniform_lon_lat) then
+ call nf95_find_coord(ncid, dimid = dimid, varid = varid, &
+ std_name = "longitude")
+ call nf95_inquire_dimension(ncid, dimid, nclen = nlon)
+ call nf95_get_var(ncid, varid, corner_deg(1))
+ call nf95_get_var(ncid, varid, lon_max, start = [nlon])
+ call nf95_find_coord(ncid, dimid = dimid, varid = varid, &
+ std_name = "latitude")
+ call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
+ call nf95_get_var(ncid, varid, corner_deg(2))
+ call nf95_get_var(ncid, varid, lat_max, start = [nlat])
+ call assert(lon_max > corner_deg(1) .and. lat_max > corner_deg(2), &
+ "input_ssh: coordinate order")
+ step_deg = [(lon_max - corner_deg(1)) / (nlon - 1), &
+ (lat_max - corner_deg(2)) / (nlat - 1)]
+ step = step_deg * deg_to_rad
+ corner_whole = corner_deg * deg_to_rad
+ else
+ call nf95_find_coord(ncid, varid = varid, std_name = "longitude")
+ call nf95_gw_var(ncid, varid, longitude)
+ longitude = longitude * deg_to_rad
+ call nf95_find_coord(ncid, varid = varid, std_name = "latitude")
+ call nf95_gw_var(ncid, varid, latitude)
+ latitude = latitude * deg_to_rad
+ call assert(longitude(2) > longitude(1) &
+ .and. latitude(2) > latitude(1), "input_ssh: coordinate order")
+ nlon = size(longitude)
+ nlat = size(latitude)
+ lon_steps = ediff1d(longitude)
+ lat_steps = ediff1d(latitude)
+ end if
+
call new_unit(unit)
open(unit, file = "grid_nml.txt", status = "replace", action = "write")
write(unit, nml = grid_nml)
close(unit)
- step = step_deg * deg_to_rad
if (periodic) call assert(4 * max_radius(1) < nlon, "input_ssh: max_radius")
copy = merge(max_radius(1), 0, periodic)
allocate(ssh(1 - copy:nlon + copy, nlat), u(1 - copy:nlon + copy, nlat), &
v(1 - copy:nlon + copy, nlat))
- corner_whole = corner_deg * deg_to_rad
! Read ssh, u and v:
diff --git a/Inst_eddies/invert_proj.f90 b/Inst_eddies/invert_proj.f90
index aa0bcdcc..6945de8f 100644
--- a/Inst_eddies/invert_proj.f90
+++ b/Inst_eddies/invert_proj.f90
@@ -1,6 +1,7 @@
module invert_proj_m
- use input_ssh_m, only: corner_whole, step
+ use input_ssh_m, only: uniform_lon_lat, corner_whole, step, longitude, &
+ latitude
implicit none
@@ -23,7 +24,11 @@ contains
!-------------------------------------------------------
- invert_proj_int = corner_whole + (ind - 1.) * step
+ if (uniform_lon_lat) then
+ invert_proj_int = corner_whole + (ind - 1.) * step
+ else
+ invert_proj_int = [longitude(ind(1)), latitude(ind(2))]
+ end if
end function invert_proj_int
@@ -32,12 +37,28 @@ contains
pure type(polyline) function invert_proj_pol(ind)
use contour_531, only: polyline, convert_to_reg_coord
+ use input_ssh_m, only: lon_steps, lat_steps
type(polyline), intent(in):: ind
+ ! Local
+ integer i, k(2)
+
!-------------------------------------------------------
- invert_proj_pol = convert_to_reg_coord(ind, corner_whole, step)
+ if (uniform_lon_lat) then
+ invert_proj_pol = convert_to_reg_coord(ind, corner_whole, step)
+ else
+ invert_proj_pol%n_points = ind%n_points
+ invert_proj_pol%closed = ind%closed
+ allocate(invert_proj_pol%points(2, ind%n_points))
+
+ do i = 1, ind%n_points
+ k = floor(ind%points(:, i))
+ invert_proj_pol%points(:, i) = [longitude(k(1)), latitude(k(2))] &
+ + (ind%points(:, i) - k) * [lon_steps(k(1)), lat_steps(k(2))]
+ end do
+ end if
end function invert_proj_pol
diff --git a/Overlap/read_grid.F90 b/Overlap/read_grid.F90
index 400ebaab..d45f94cd 100644
--- a/Overlap/read_grid.F90
+++ b/Overlap/read_grid.F90
@@ -8,11 +8,11 @@ module read_grid_m
logical, protected:: periodic = .false. ! grid is periodic in longitude
- real, protected:: corner_whole(2)
+ real, protected:: corner_whole(2) = huge(0.)
! longitude and latitude of the corner of the whole grid in input
! NetCDF file, in rad
- real, protected:: step(2) ! longitude and latitude steps, in rad
+ real, protected:: step(2) = huge(0.) ! longitude and latitude steps, in rad
logical, protected:: uniform_lon_lat = .true.
! grid is cartesian and uniform in longitude and latitude
@@ -52,17 +52,23 @@ contains
action = "read", position = "rewind")
read(unit, nml = grid_nml)
close(unit)
- corner_whole = corner_deg * deg_to_rad
- step = step_deg * deg_to_rad
+
+ if (uniform_lon_lat) then
+ corner_whole = corner_deg * deg_to_rad
+ step = step_deg * deg_to_rad
+ end if
end if
#ifndef CPP_SEQUENTIAL
call ezmpi_bcast(nlon, root = 0)
call ezmpi_bcast(nlat, root = 0)
call ezmpi_bcast(periodic, root = 0)
- call ezmpi_bcast(corner_whole, root = 0)
- call ezmpi_bcast(step, root = 0)
call ezmpi_bcast(uniform_lon_lat, root = 0)
+
+ if (uniform_lon_lat) then
+ call ezmpi_bcast(corner_whole, root = 0)
+ call ezmpi_bcast(step, root = 0)
+ end if
#endif
end subroutine read_grid
--
GitLab