get_snapshot.f

module get_snapshot_m

  implicit none

contains

  subroutine get_snapshot(s, min_amp, min_area, m, n_proc, k_end, max_delta, &
       nlon, nlat, k, max_radius, corner, step, periodic)

    use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN

    ! Libraries:
    use contour_531, only: convert_to_ind, null_polyline
    use netcdf, only: nf90_nowrite
    use netcdf95, only: nf95_open, nf95_close
    use nr_util, only: twopi

    use derived_types, only: snapshot, null_ssh_contour, missing_speed
    use get_var_m, only: get_var
    use nearby_extr_m, only: nearby_extr
    use receive_snapshot_m, only: receive_snapshot
    use set_max_speed_m, only: set_max_speed
    use set_all_outerm_m, only: set_all_outerm

    type(snapshot), intent(out):: s

    real, intent(in):: min_amp
    ! minimum amplitude of ssh, between extremum and outermost
    ! contour, in m

    real, intent(in):: min_area
    ! minimum area of an outermost contour, in m2
    
    integer, intent(in):: m ! rank of MPI process
    integer, intent(in):: n_proc ! number of MPI processes
    integer, intent(in):: k_end ! last date index analyzed by this MPI process

    integer, intent(in):: max_delta
    ! maximum interval of date indices between which we look for
    ! overlapping of eddies

    integer, intent(in):: nlon, nlat
    integer, intent(in):: k ! date index

    integer, intent(in):: max_radius(:) ! (2) maximum radius of an
    ! eddy in longitude and latitude, in number of grid points

    real, intent(in):: corner(:) ! (2) longitude and latitude of the
    ! corner of the whole grid, in rad
    
    real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
    logical, intent(in):: periodic ! grid is periodic in longitude

    ! Local:

    integer ncid
    real ssh(1 - merge(max_radius(1), 0, periodic):nlon + merge(max_radius(1), &
         0, periodic), nlat) ! sea-surface height, in m
    real, dimension(1 - merge(max_radius(1), 0, periodic):nlon &
         + merge(max_radius(1), 0, periodic), nlat):: u, v ! wind, in m s-1
    integer i

    ! Window around each extremum:
    
    integer llc(2) ! indices in global grid of lower left corner
    integer urc(2) ! indices in global grid of upper right corner

    real corner_window(2) ! longitude and latitude of the window around each
    ! extremum, in rad

    real, allocatable:: outside_points(:, :) ! (2, :) longitude and
    ! latitude, in rad, of all the significant extrema, except the
    ! target extremum

    !--------------------------------------------------------------

    if (m == n_proc .or. k <= k_end - max_delta) then
       ! Read ssh, u and v at date k:

       call nf95_open("h.nc", nf90_nowrite, ncid)
       call get_var(periodic, max_radius(1), ssh, ncid, nlon, k, name = "adt", &
            new_fill_value = huge(0.))
       ! (We cannot keep the original fill value because Contour_531
       ! works with an upper limit for valid values.)
       call nf95_close(ncid)

       call nf95_open("uv.nc", nf90_nowrite, ncid)
       call get_var(periodic, max_radius(1), u, ncid, nlon, k, name = "u", &
            new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
       call get_var(periodic, max_radius(1), v, ncid, nlon, k, name = "v", &
            new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
       ! (We will need quiet NaNs rather the original fill values for
       ! u and v when we compute the max-speed contours and when we
       ! search the ssh of max-speed contours.)
       call nf95_close(ncid)

       call set_all_outerm(s, min_amp, max_radius, step, periodic, ssh, &
            corner, min_area)

       ! Done with outermost contours, now let us take care of the
       ! max-speed contours.

       do i = 1, s%number_vis_extr
          if (s%list_vis(i)%valid) then
             ! Restrict the field to the outermost contour:

             llc = floor(convert_to_ind(minval(s%list_vis(i)%out_cont%points, &
                  dim = 2), corner, step))

             urc = ceiling(convert_to_ind(maxval( &
                  s%list_vis(i)%out_cont%points, dim = 2), corner, step))

             ! Should have no effect except because of roundup error:
             urc(2) = min(urc(2), nlat)
             if (.not. periodic) urc(1) = min(urc(1), nlon)

             corner_window = corner + (llc - 1) * step
             
             outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
                  llc(2):urc(2)), s%list_vis, i)

             if (periodic) outside_points(1, :) = outside_points(1, :) &
                  + ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
                  * twopi
             ! (Shift the longitude of each point to a value close to
             ! the longitude of the target extremum.)
             
             call set_max_speed(s%list_vis(i), s%ind_extr(:, i) - llc + 1, &
                  outside_points, ssh(llc(1):urc(1), llc(2):urc(2)), &
                  u(llc(1):urc(1), llc(2):urc(2)), &
                  v(llc(1):urc(1), llc(2):urc(2)), corner_window, step)
          else
             s%list_vis(i)%speed_cont = null_ssh_contour()
             s%list_vis(i)%max_speed = missing_speed
             s%list_vis(i)%radius4 = 0
          end if
       end do

       s%number_eddies = s%number_vis_extr
    else
       call receive_snapshot(s, source = m + 1, tag = k)
    end if

  end subroutine get_snapshot

end module get_snapshot_m