end do    end do!$OMP END PARALLEL DO! loop through overlapping cells on input grid to make area-average    do n = 1, nmax!$OMP PARALLEL DO PRIVATE (jo,io,ii,ji)       do jo = 1, nlat_o          do io =1, numlon_o(jo)             ii = i_ovr(io,jo,n)             ji = j_ovr(io,jo,n)             fld_o(io,jo) = fld_o(io,jo) + w_ovr(io,jo,n)*fld_i(ii,ji)          end do       end do!$OMP END PARALLEL DO    end do! set non-valid points for reduced grid to missing value    do jo = 1, nlat_o       do io = numlon_o(jo)+1, nlon_o          fld_o(io,jo) = 1.e36       end do    end do    return  end subroutine areaave!=======================================================================  subroutine areamap (nlon_i   , nlat_i   , nlon_o , nlat_o ,  &                      lon_i    , lat_i    , lon_o  , lat_o  ,  &                      numlon_i , numlon_o , mask_i , mx_ovr ,  &                      n_ovr    , i_ovr    , j_ovr  , w_ovr  ,  &                      fland_o  , area_o   )!----------------------------------------------------------------------- ! ! Purpose: ! weights and indices for area of overlap between grids!! Method: ! Get indices and weights for area-averaging between input and output grids.! For each output grid cell find:!!    o number of input grid cells that overlap with output grid cell (n_ovr)!    o longitude index (1 <= i_ovr <= nlon_i) of the overlapping input grid cell!    o latitude index  (1 <= j_ovr <= nlat_i) of the overlapping input grid cell!    o fractional overlap of input grid cell (w_ovr)!! so that for!! field values fld_i on an  input grid with dimensions nlon_i and nlat_i! field values fld_o on an output grid with dimensions nlon_o and nlat_o are!! fld_o(io,jo) = ! fld_i(i_ovr(io,jo,     1),j_ovr(io,jo,     1)) * w_ovr(io,jo,     1) +!                             ... + ... +! fld_i(i_ovr(io,jo,mx_ovr),j_ovr(io,jo,mx_ovr)) * w_ovr(io,jo,mx_ovr) !! Note: mx_ovr is some number greater than n_ovr. Weights of zero are ! used for the excess points! ! Author: Gordon Bonan! !-----------------------------------------------------------------------! ------------------------ arguments ---------------------------------    integer ,intent(in) :: nlon_i                     !input grid : max number of longitude points    integer ,intent(in) :: nlat_i                     !input grid : number of latitude points    integer ,intent(in) :: numlon_i(nlat_i)           !input grid : number longitude points for lat    real(r8),intent(inout) :: lon_i(nlon_i+1,nlat_i)  !input grid : cell longitude, west edge (deg)    real(r8),intent(in) :: lat_i(nlat_i+1)            !input grid : cell latitude,  south edge (deg)    real(r8),intent(in) :: mask_i(nlon_i,nlat_i)      !input grid : mask (0, 1)     integer ,intent(in) :: nlon_o                     !output grid: max number of longitude points    integer ,intent(in) :: nlat_o                     !output grid: number of latitude points    integer ,intent(in) :: numlon_o(nlat_o)           !output grid: number longitude points for lat    real(r8),intent(in) :: lon_o(nlon_o+1,nlat_o)     !output grid: cell longitude, west edge  (deg)    real(r8),intent(in) :: lat_o(nlat_o+1)            !output grid: cell latitude,  south edge (deg)    real(r8),intent(in) :: fland_o(nlon_o,nlat_o)     !output grid: fraction that is land    real(r8),intent(in) :: area_o(nlon_o,nlat_o)      !output grid: cell area    integer ,intent(in) :: mx_ovr                     !max num input cells that overlap output cell    integer ,intent(out):: n_ovr(nlon_o,nlat_o)       !number of overlapping input cells    integer ,intent(out):: i_ovr(nlon_o,nlat_o,mx_ovr)!lon index, overlapping input cell    integer ,intent(out):: j_ovr(nlon_o,nlat_o,mx_ovr)!lat index, overlapping input cell    real(r8),intent(out):: w_ovr(nlon_o,nlat_o,mx_ovr)!overlap weights for input cells! --------------------------------------------------------------------! ------------------------ local variables ---------------------------    integer :: io                   !output grid longitude loop index    integer :: ii                   !input  grid longitude loop index    integer :: jo                   !output grid latitude  loop index    integer :: ji                   !input  grid latitude  loop index    integer :: n                    !overlapping cell index    real(r8) :: offset              !used to shift x-grid 360 degrees    real(r8) :: f_ovr               !sum of overlap weights     real(r8) :: relerr = 0.00001    !max error: sum overlap weights ne 1    real(r8) :: dx_i                !input grid  longitudinal range    real(r8) :: dy_i                !input grid  latitudinal  range    real(r8) :: dx_o                !output grid longitudinal range    real(r8) :: dy_o                !output grid latitudinal  range! --------------------------------------------------------------------! --------------------------------------------------------------------! Initialize overlap weights on output grid to zero for maximum ! number of overlapping points. Set lat and lon indices of overlapping! input cells to dummy values. Set number of overlapping cells to zero! --------------------------------------------------------------------    do n = 1, mx_ovr        do jo = 1, nlat_o                do io = 1, numlon_o(jo)               i_ovr(io,jo,n) = 1             j_ovr(io,jo,n) = 1             w_ovr(io,jo,n) = 0.           end do       end do    end do    do jo = 1, nlat_o             do io = 1, numlon_o(jo)            n_ovr(io,jo) = 0       end do    end do! --------------------------------------------------------------------! First pass to find cells that overlap and area of overlap! --------------------------------------------------------------------    call areaovr (nlon_i , nlat_i , numlon_i , lon_i  , lat_i  , &                  nlon_o , nlat_o , numlon_o , lon_o  , lat_o  , &                  mx_ovr , n_ovr  , i_ovr    , j_ovr  , w_ovr  )! --------------------------------------------------------------------! Second pass to find cells that overlap and area of overlap! --------------------------------------------------------------------! Shift x-grid to locate periodic grid intersections. This! assumes that all lon_i(1,j) have the same value for all! latitudes j and that the same holds for lon_o(1,j)    if (lon_i(1,1) < lon_o(1,1)) then       offset = 360.0    else       offset = -360.0    end if    do ji = 1, nlat_i       do ii = 1, numlon_i(ji) + 1          lon_i(ii,ji) = lon_i(ii,ji) + offset       end do    end do! find overlap    call areaovr (nlon_i , nlat_i , numlon_i , lon_i  , lat_i  , &                  nlon_o , nlat_o , numlon_o , lon_o  , lat_o  , &                  mx_ovr , n_ovr  , i_ovr    , j_ovr  , w_ovr  )! restore x-grid (un-shift x-grid)    do ji = 1, nlat_i       do ii = 1, numlon_i(ji) + 1          lon_i(ii,ji) = lon_i(ii,ji) - offset       end do    end do! --------------------------------------------------------------------! Normalize areas of overlap to get fractional contribution of each! overlapping grid cell (input grid) to grid cell average on output grid.! Normally, do this by dividing area of overlap by area of output grid cell.! But, only have data for land cells on input grid. So if output grid cell! overlaps with land and non-land cells (input grid), do not have valid! non-land data for area-average. Instead, weight by area of land using! [mask_i], which has a value of one for land and zero for ocean. If! [mask_i] = 1, input grid cell contributes to output grid cell average.! If [mask_i] = 0, input grid cell does not contribute to output grid cell ! average.! --------------------------------------------------------------------    do jo = 1, nlat_o       do io = 1, numlon_o(jo)! find total land area of overlapping input cells          f_ovr = 0.          do n = 1, n_ovr(io,jo)             ii = i_ovr(io,jo,n)             ji = j_ovr(io,jo,n)             f_ovr = f_ovr + w_ovr(io,jo,n)*mask_i(ii,ji)          end do! make sure area of overlap is less than or equal to output grid cell area          if ((f_ovr-area_o(io,jo))/area_o(io,jo) > relerr) then             write (6,*) 'AREAMAP error: area not conserved for lon,lat = ',io,jo             write (6,'(a30,e20.10)') 'sum of overlap area = ',f_ovr             write (6,'(a30,e20.10)') 'area of output grid = ',area_o(io,jo)             call endrun          end if! make weights          do n = 1, n_ovr(io,jo)             ii = i_ovr(io,jo,n)             ji = j_ovr(io,jo,n)             if (f_ovr > 0.) then                w_ovr(io,jo,n) = w_ovr(io,jo,n)*mask_i(ii,ji) / f_ovr             else                w_ovr(io,jo,n) = 0.             end if          end do       end do    end do! --------------------------------------------------------------------! Error check: overlap weights for input grid cells must sum to 1. This! is always true if both grids span the same domain. However, if one! grid is a subset of the other grid, this is only true when mapping! from the full grid to the subset. When input grid covers a smaller! domain than the output grid, this test is not valid.! --------------------------------------------------------------------    dx_i = lon_i(nlon_i+1,1) - lon_i(1,1)    dx_o = lon_o(nlon_o+1,1) - lon_o(1,1)    if (lat_i(nlat_i+1) > lat_i(1)) then      !South to North grid       dy_i = lat_i(nlat_i+1) - lat_i(1)    else                                         !North to South grid       dy_i = lat_i(1) - lat_i(nlat_i+1)    end if    if (lat_o(nlat_o+1) > lat_o(1)) then      !South to North grid       dy_o = lat_o(nlat_o+1) - lat_o(1)    else                                         !North to South grid       dy_o = lat_o(1) - lat_o(nlat_o+1)    end if    if (abs(dx_i-dx_o)>relerr .or. abs(dy_i-dy_o)>relerr) then       if (dx_i<dx_o .or. dy_i<dy_o) then          write (6,*) 'AREAMAP warning: area-average not valid for '          write (6,*) '   input  grid of ',nlon_i,' x ',nlat_i          write (6,*) '   output grid of ',nlon_o,' x ',nlat_o          return            end if    end if    do jo = 1, nlat_o       do io = 1, numlon_o(jo)          f_ovr = 0.          do n = 1, mx_ovr             f_ovr = f_ovr + w_ovr(io,jo,n)          end do! error check only valid if output grid cell has land. non-land cells! will have weights equal to zero          if (fland_o(io,jo) > 0.) then             if (abs(f_ovr-1.) > relerr) then                write (6,*) 'AREAMAP error: area not conserved for lon,lat = ',io,jo                write (6,'(a30,e20.10)') 'sum of overlap weights = ',f_ovr                call endrun             end if          end if       end do    end do    return  end subroutine areamap!=======================================================================  subroutine areaovr (nlon_i , nlat_i , numlon_i , lon_i , lat_i  , &                      nlon_o , nlat_o , numlon_o , lon_o , lat_o  , &                      mx_ovr , n_ovr  , i_ovr    , j_ovr , w_ovr  )!----------------------------------------------------------------------- ! ! Purpose: ! Find area of overlap between grid cells!! Method: ! For each output grid cell: find overlapping input grid cell and area of! input grid cell that overlaps with output grid cell. Cells overlap if:!! southern edge of input grid < northern edge of output grid AND! northern edge of input grid > southern edge of output grid!! western edge of input grid < eastern edge of output grid AND! eastern edge of input grid > western edge of output grid!