Skip to content
Snippets Groups Projects
Commit 4399e353 authored by Lionel GUEZ's avatar Lionel GUEZ Committed by GUEZ Lionel
Browse files

Compute good contour in projection coordinates

Procedure `good_contour` now takes input in projection coordinates and
returns a contour in projection coordinates. The dummy argument step
of `good_contour` is thus removed.
parent b4355bef
No related branches found
No related tags found
No related merge requests found
&MAIN_NML
INSIDE_POINT=-0.7,0
INSIDE_POINT=-8.75,0.
LEVEL=5
STEP=2*0.08
CORNER=-2,-2 /
CORNER=-25.,-25. /
&MAIN_NML INSIDE_POINT=0.3,0,
&MAIN_NML INSIDE_POINT=3.75,0,
LEVEL=5
corner = -2., -2., step = 0.08, 0.08 /
corner = -25., -25. /
&MAIN_NML LEVEL= 5, corner = -2., -2., step = 0.08, 0.08 /
&MAIN_NML LEVEL= 5, corner = -25., -25. /
-0.5, -1.
0.3, 0.
-6.25, -12.5
3.75, 0.
-0.5, -1.
-0.7, 0.
-6.25, -12.5
-8.75, 0.
......@@ -21,7 +21,6 @@ program test_good_contour
real, allocatable:: Z(:, :)
real:: inside_point(2) = [0., 0.]
real, allocatable:: outside_points(:, :) ! (2, n_out)
real:: step(2) = [1., 1.]
real:: corner(2) = [1., 1.]
TYPE(shpfileobject) shphandle
integer field_number, shape_number
......@@ -29,7 +28,7 @@ program test_good_contour
integer, allocatable:: dimids(:)
character(len = 10):: varname = "z"
namelist /main_nml/ inside_point, level, step, corner, varname
namelist /main_nml/ inside_point, level, corner, varname
!-------------------------------------------------------------------
......@@ -70,7 +69,7 @@ program test_good_contour
print *, outside_points(:, i)
end do
c = good_contour(corner, step, z, level, inside_point, outside_points)
c = good_contour(corner, z, level, inside_point, outside_points)
if (c%n_points /= 0) then
call shp_create_03("test_good_contour", shpt_polygon, shphandle)
......
......@@ -41,13 +41,14 @@ contains
! 360 degrees.
! Libraries:
use contour_531, only: polyline
use contour_531, only: polyline, convert_to_reg_coord, convert_to_ind
use jumble, only: assert
use ccw_orient_m, only: ccw_orient
use config_m, only: min_area
use derived_types, only: ssh_contour, null_ssh_contour
use good_contour_m, only: good_contour
use input_ssh_m, only: corner_whole => corner
use spher_polyline_area_m, only: spher_polyline_area
type(ssh_contour), intent(out):: out_cont
......@@ -88,13 +89,20 @@ contains
real, parameter:: precision = 5e-4 ! in m
logical mask(size(ssh, 1), size(ssh, 2))
integer n_max_cont ! >= 3
integer i
real corner_proj(2), coord_extr_proj(2)
real outside_points_proj(2, size(outside_points, 2))
!-----------------------------------------------------------------
corner_proj = (corner - corner_whole) / step + 1.
coord_extr_proj = (coord_extr - corner_whole) / step + 1.
outside_points_proj = convert_to_ind(outside_points, corner_whole, step)
n_max_cont = size(cont_list)
n_cont = 0
cont_list(1)%polyline = good_contour(corner, step, ssh, innermost_level_2, &
coord_extr, outside_points)
cont_list(1)%polyline = convert_to_reg_coord(good_contour(corner_proj, &
ssh, innermost_level_2, coord_extr_proj, outside_points_proj), &
corner_whole, step)
test_n_points: if (cont_list(1)%n_points == 0) then
out_cont = null_ssh_contour()
......@@ -106,8 +114,9 @@ contains
level_try = merge(maxval(ssh, mask), minval(ssh, mask), cyclone)
call assert(merge(level_try > level_good, level_try < level_good, &
cyclone), "get_1_outerm level_try")
cont_list(2)%polyline = good_contour(corner, step, ssh, level_try, &
coord_extr, outside_points)
cont_list(2)%polyline = convert_to_reg_coord(good_contour(corner_proj, &
ssh, level_try, coord_extr_proj, outside_points_proj), &
corner_whole, step)
if (cont_list(2)%n_points /= 0) then
! This can only happen if there are degenerate extrema
......@@ -121,8 +130,10 @@ contains
do while (abs(level_bad - level_good) > precision)
level_try = (level_good + level_bad) / 2.
cont_list(n_cont + 1)%polyline = good_contour(corner, step, ssh, &
level_try, coord_extr, outside_points)
cont_list(n_cont + 1)%polyline &
= convert_to_reg_coord(good_contour(corner_proj, ssh, &
level_try, coord_extr_proj, outside_points_proj), &
corner_whole, step)
if (cont_list(n_cont + 1)%n_points /= 0) then
cont_list(n_cont + 1)%ssh = level_try
......
......@@ -4,7 +4,7 @@ module good_contour_m
contains
type(polyline) function good_contour(corner, step, z, level, inside_point, &
type(polyline) function good_contour(corner, z, level, inside_point, &
outside_points)
! Finds a contour in the plane, if it exists, which is closed,
......@@ -20,8 +20,6 @@ contains
use geometry, only: polygon_contains_point
real, intent(in):: corner(:) ! (2) [x, y] corresponding to z(1, 1)
real, intent(in):: step(:) ! (2)
real, intent(in):: z(:, :)
! We are assuming that significant values of z cannot be greater
! than zmax and that the missing value flag in z is greater than
......@@ -48,7 +46,7 @@ contains
n_out = 0
end if
call find_contours_reg_grid(corner, step, z, level, contours, zmax)
call find_contours_reg_grid(corner, [1., 1.], z, level, contours, zmax)
n_cont = size(contours)
found_good_contour = .false.
i = 1
......
......@@ -22,11 +22,13 @@ contains
use jumble, only: rad_to_deg
use shapelib, only: shpt_polygon
use shapelib_03, only: shp_append_object_03, dbf_write_attribute_03
use contour_531, only: convert_to_reg_coord
use ccw_orient_m, only: ccw_orient
use complete_ssh_m, only: complete_ssh
use derived_types, only: eddy, null_ssh_contour, ssh_contour, missing_speed
use good_contour_m, only: good_contour
use input_ssh_m, only: corner_whole => corner
use mean_speed_m, only: mean_speed
use spher_polyline_area_m, only: spher_polyline_area
use cont_list_m, only: write_cont_list, hshp_cont_list, ifield_ssh, &
......@@ -71,16 +73,19 @@ contains
real, allocatable:: speed(:) ! (n_cont) speed on the contour
integer i, i_outer, ishape
real corner_proj(2), extr_coord_proj(2)
!---------------------------------------------------------------
corner_proj = (corner - corner_whole) / step + 1.
extr_coord_proj = (extr_coord - corner_whole) / step + 1.
i_outer = n_cont
if (n_cont >= 2) call complete_ssh(cont_list%ssh, n_cont)
! Now find the contours associated to the new values of SSH:
do i = i_outer + 1, n_cont
cont_list(i)%polyline = good_contour(corner, step, ssh, &
cont_list(i)%ssh, extr_coord)
cont_list(i)%polyline = convert_to_reg_coord(good_contour(corner_proj, &
ssh, cont_list(i)%ssh, extr_coord_proj), corner_whole, step)
end do
allocate(speed(n_cont))
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment