diff --git a/Inst_eddies/Tests/test_set_max_speed.f90 b/Inst_eddies/Tests/test_set_max_speed.f90 index b3512d21c604de77616326b592224ee44e7cdcc5..39c3437b5e7e9d84182384f99626b44fb25c782b 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 8cb6c176debbc611adc7d68bf881e0600cc6eb32..b3d0115bb9f602ae8f3dbfa83ce374f8642691a0 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 aa0bcdccdfe5af6647a351dc19e62998fe7969d2..6945de8fa61b8c5de53bfa31d2784045b460139f 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 400ebaabcbb07caf2b9b7e6c57fe21988a16b1ba..d45f94cdcd775df276d625b4535f21dc48d3ec7b 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