Updated I have modified my solutions based on what your have requested in comments. This time I first found the common months for each Plot and then went on to calculate the mean for every year.

library(dplyr)
library(purrr)

# I created a character vector of Plot names in order to ease iteration
Plots <- paste("Plot", c(1:3), sep = "")

# Then we iterate over each name selecting only the corresponding column 
# from the original data set and find the common months for each iteration.
Plots %>%
  map(~ df %>% 
        select(.x, months, year) %>%
        drop_na() %>%
        group_split(year) %>%
        map(~ .x %>% 
              select(months) %>%
              distinct()) %>% 
        reduce(~ intersect(.x, .y))) %>%
  set_names(nm = Plots) -> commons

Then we use the commons' months to apply a filter on the entire data set:

# This time we iterate over both commons and plot names at the same time

map2(Plots, commons, ~ 
       df %>% 
       select(.x, months, year) %>%
       filter(months %in% .y[["months"]]) %>%
       group_by(year) %>%
       summarise(across(starts_with("Plot"), ~ mean(.x, na.rm = TRUE)))) %>%
  exec(cbind, !!!.) %>%
  select(-c(3, 5))


    year Plot1 Plot2 Plot3
   <dbl> <dbl> <dbl> <dbl>
 1  2000 0.749 0.771 0.566
 2  2001 0.78  0.724 0.512
 3  2002 0.741 0.707 0.528
 4  2003 0.703 0.806 0.559
 5  2004 0.734 0.757 0.549
 6  2005 0.725 0.758 0.547
 7  2006 0.732 0.796 0.542
 8  2007 0.725 0.763 0.585
 9  2008 0.761 0.736 0.569
10  2009 0.766 0.776 0.575
11  2010 0.768 0.78  0.622
12  2011 0.781 0.787 0.549
13  2012 0.776 0.784 0.553
14  2013 0.774 0.787 0.598
15  2014 0.748 0.809 0.577
16  2015 0.731 0.797 0.587
17  2016 0.758 0.67  0.653
18  2017 0.799 0.76  0.599
19  2018 0.794 0.81  0.575

You have three Plot columns and you want to do the same thing with each one. In cases like this you can make your life easier by 'pivoting' the data into 'long' (i.e. 'tidy') format. Here is an approach using tidyr and dplyr:

library(dplyr)
library(tidyr)

tidied_df <- df %>%
  
  # First pivot into 'long' or 'tidy' format
  pivot_longer(Plot1:Plot3, names_to = "Plot", 
               values_to = "Value", names_prefix = "Plot") %>% 
  
  # First remove all missing values of 'Value'
  filter(!is.na(Value)) %>% 
  
  # Now turn 'implicit' missing months to explicit
  complete(Plot, year, months) %>% 
  
  # Can now add a column indicating whether the whole month is 'missing'
  # for a given year
  group_by(Plot, year, months) %>% 
  mutate(month_missing = all(is.na(Value))) %>% 
  
  # Now use the helper column to remove a month from all years if any
  # one year has only missing values for that month
  group_by(Plot, months) %>% 
  filter(!any(month_missing)) %>% 
  ungroup()
  
# We can take a look to see which months are present across all years
# for each plot like so:
tidied_df %>% 
  group_by(Plot) %>% 
  summarise(
    months = paste(sort(unique(months)), collapse = ", "), 
    .groups = "drop"
  )
#> # A tibble: 3 x 2
#>   Plot  months                         
#>   <chr> <chr>                          
#> 1 1     2, 3, 4, 5, 6, 7, 9, 10, 11, 12
#> 2 2     4, 5, 6, 7, 8, 9, 10, 11       
#> 3 3     3, 4, 5, 6, 7, 8, 9, 10
  

# Finally, we can compute the mean across each plot:
tidied_df %>% 
  group_by(Plot, year) %>% 
  summarise(Value = mean(Value, na.rn = TRUE), .groups = "drop")
#> # A tibble: 57 x 3
#>    Plot   year Value
#>    <chr> <dbl> <dbl>
#>  1 1      2000 0.749
#>  2 1      2001 0.78 
#>  3 1      2002 0.741
#>  4 1      2003 0.703
#>  5 1      2004 0.734
#>  6 1      2005 0.725
#>  7 1      2006 0.732
#>  8 1      2007 0.725
#>  9 1      2008 0.761
#> 10 1      2009 0.766
#> # … with 47 more rows

Edit

As requested, here is a method to filter for values of Plot where all years have at least 4 months in common:

tidied_df %>% 
  group_by(Plot) %>% 
  filter(n_distinct(months) > 4) %>% 
  ungroup()