! Method: ! ! Author: Mariana Vertenstein! !-----------------------------------------------------------------------! ------------------------ arguments ---------------------------------    integer , intent(in) :: nlat              !dimension: number of latitude points    integer , intent(in) :: nlon              !dimension: number of longitude points    integer , intent(in) :: numlon(nlat)      !number of grid cells per latitude strip    real(r8), intent(in) :: longxy(nlon,nlat) !longitude at center of grid cell    real(r8), intent(in) :: latixy(nlon,nlat) !latitude at center of grid cell    real(r8), intent(in) :: edgen             !northern edge of grid (degrees)    real(r8), intent(in) :: edgee             !eastern edge of grid (degrees)    real(r8), intent(in) :: edges             !southern edge of grid (degrees)    real(r8), intent(in) :: edgew             !western edge of grid (degrees)    real(r8), intent(out):: lats(nlat+1)      !grid cell latitude, southern edge (degrees)    real(r8), intent(out):: lonw(nlon+1,nlat) !grid cell longitude, western edge (degrees)! --------------------------------------------------------------------! ------------------- local variables -----------------------------    integer i,j             !indices    real(r8) dx             !cell width! --------------------------------------------------------------------! --------------------------------------------------------------------! Latitudes -- southern/northern edges for each latitude strip. ! For grids oriented South to North, the southern! and northern edges of latitude strip [j] are:!        southern = lats(j  )!        northern = lats(j+1)! For grids oriented North to South: the southern! and northern edges of latitude strip [j] are:!        northern = lats(j  )!        southern = lats(j+1)! In both cases, [lats] must be dimensioned lats(lat+1)! --------------------------------------------------------------------    if (nlat == 1) then       lats(1) = edges        lats(nlat+1) = edgen     else       if (latixy(1,2) > latixy(1,1)) then    !South to North grid          lats(1) = edges           lats(nlat+1) = edgen        else                                   !North to South grid          lats(1) = edgen           lats(nlat+1) = edges        end if    end if    do j = 2, nlat       lats(j) = (latixy(1,j-1) + latixy(1,j)) / 2.    end do! --------------------------------------------------------------------! Longitudes -- western edges. Longitudes for the western edge of the ! cells must increase continuously and span 360 degrees. Assume that! grid starts at Dateline with western edge on Dateline Western edges ! correspond to [longxy] (longitude at center of cell) and range from ! -180 to 180 with negative longitudes west of Greenwich. ! Partial grids that do not span 360 degrees are allowed so long as they! have the convention of Grid 1 with !      western edge of grid: >= -180 and < 180!      eastern edge of grid: > western edge  and <= 180! [lonw] must be dimensioned lonw(lon+1,lat) because each latitude! strip can have variable longitudinal resolution! --------------------------------------------------------------------! Western edge of first grid cell -- since grid starts with western! edge on Dateline, lonw(1,j)=-180. This is the same as [edgew].! Remaining grid cells. On a global grid lonw(numlon+1,j)=lonw(1,j)+360. ! This is the same as [edgee].  Set unused longitudes to non-valid number    do j = 1, nlat       dx = (edgee - edgew) / numlon(j)       lonw(1,j) = edgew       do i = 2, numlon(j)+1          lonw(i,j) = lonw(1,j) + (i-1)*dx       end do       do i = numlon(j)+2, nlon          lonw(i,j) = -999.       end do    end do    return  end subroutine celledge_regional!=======================================================================  subroutine celledge_global (nlat, nlon, numlon, longxy, latixy, &                              lats, lonw )!----------------------------------------------------------------------- ! ! Purpose: ! Southern and western edges of grid cells - global grid!! Method: ! ! Author: Mariana Vertenstein! !-----------------------------------------------------------------------! ------------------------ arguments ---------------------------------    integer , intent(in) :: nlat               !dimension: number of latitude points    integer , intent(in) :: nlon               !dimension: number of longitude points    integer , intent(in) :: numlon(nlat)       !number of grid cells per latitude strip    real(r8), intent(in) :: longxy(nlon,nlat)  !longitude at center of grid cell    real(r8), intent(in) :: latixy(nlon,nlat)  !latitude at center of grid cell    real(r8), intent(out):: lats(nlat+1)       !grid cell latitude, southern edge (degrees)    real(r8), intent(out):: lonw(nlon+1,nlat)  !grid cell longitude, western edge (degrees)! --------------------------------------------------------------------! ------------------- local variables -----------------------------    integer i,j             !indices    real(r8) dx             !cell width! --------------------------------------------------------------------! --------------------------------------------------------------------! Latitudes -- southern/northern edges for each latitude strip. ! For grids oriented South to North, the southern! and northern edges of latitude strip [j] are:!        southern = lats(j  )!        northern = lats(j+1)! For grids oriented North to South: the southern! and northern edges of latitude strip [j] are:!        northern = lats(j  )!        southern = lats(j+1)! In both cases, [lats] must be dimensioned lats(lat+1)! --------------------------------------------------------------------    do j = 1, nlat+1              !southern edges       if (j == 1) then           !south pole          lats(j) = -90.       else if (j == nlat+1) then !north pole          lats(j) = 90.       else                       !edge = average latitude          lats(j) = (latixy(1,j-1) + latixy(1,j)) / 2.       end if    end do! --------------------------------------------------------------------! Longitudes -- western edges. Longitudes for the western edge of the ! cells must increase continuously and span 360 degrees.! --------------------------------------------------------------------    if (longxy(1,1) >= 0.) then       do j = 1, nlat          dx = 360./(numlon(j))          do i = 1, numlon(j)+1             lonw(i,j) = -dx/2. + (i-1)*dx          end do          do i = numlon(j)+2, nlon             lonw(i,j) = -999.          end do       end do    else       write(6,*)'global non-regional grids currently only supported ', &            'for grids starting at greenwich and centered on Greenwich'       call endrun    endif    return  end subroutine celledge_global!=======================================================================  subroutine areaini_point (io      , jo          , nlon_i  , nlat_i  , numlon_i, &                            lon_i   , lon_i_offset, lat_i   , area_i  , mask_i  , &                            nlon_o  , nlat_o      , numlon_o, lon_o   , lat_o   , &                            area_o  , fland_o     , novr_i2o, iovr_i2o, jovr_i2o, &                            wovr_i2o)!----------------------------------------------------------------------- ! ! Purpose: ! area averaging initialization !! Method: ! This subroutine is used in conjunction with areaave.F for area-average! mapping of a field from one grid to another. !!    areaini_point  - initializes indices and weights for area-averaging from !                     input grid to output grid!    areamap_point  - called by areaini: finds indices and weights!    areaovr_onit   - called by areamap: finds if cells overlap and area of overlap! ! To map from one grid to another, must first call areaini to build! the indices and weights (iovr_i2o, jovr_i2o, wovr_i2o). Then must! call areaave to get new field on output grid.!! Not all grid cells on the input grid will be used in the area-averaging! of a field to the output grid. Only input grid cells with [mask_i] = 1! contribute to output grid cell average. If [mask_i] = 0, input grid cell ! does not contribute to output grid cell. This distinction is not usually! required for atm -> land mapping, because all cells on the atm grid have! data. But when going from land -> atm, only land grid cells have data.! Non-land grid cells on surface grid do not have data. So if output grid cell! overlaps with land and non-land cells (input grid), can only use land! grid cells when computing area-average. !! o Input and output grids can be ANY resolution BUT:!!   a. Grid orientation -- Grids can be oriented south to north!      (i.e. cell(lat+1) is north of cell(lat)) or from north to !      south (i.e. cell(lat+1) is south of cell(lat)). Both grids must be !      oriented from west to east, i.e., cell(lon+1) must be east of cell(lon)!!   b. Grid domain -- Grids do not have to be global. Both grids are defined!      by their north, east, south, and west edges (edge_i and edge_o in!      this order, i.e., edge_i(1) is north and edge_i(4) is west).!      !      For partial grids, northern and southern edges are any latitude!      between 90 (North Pole) and -90 (South Pole). Western and eastern!      edges are any longitude between -180 and 180, with longitudes !      west of Greenwich negative.!!      For global grids, northern and southern edges are 90 (North Pole) !      and -90 (South Pole). The grids do not have to start at the!      same longitude, i.e., one grid can start at Dateline and go east;!      the other grid can start at Greenwich and go east. Longitudes for!      the western edge of the cells must increase continuously and span!      360 degrees. Examples!!                              West edge    East edge!                            ---------------------------------------------------!      Dateline            :        -180 to 180        (negative W of Greenwich)!      Greenwich (centered):    0 - dx/2 to 360 - dx/2 !!   c. Both grids can have variable number of longitude points for each!      latitude strip. However, the western edge of the first point in each!      latitude must be the same for all latitudes. Likewise, for the!      eastern edge of the last point. That is, each latitude strip must span !      the same longitudes, but the number of points to do this can be different!!   d. One grid can be a sub-set (i.e., smaller domain) than the other grid.!      In this way, an atmospheric dataset for the entire globe can be!      used in a simulation for a region 30N to 50N and 130W to 70W -- the !      code will extract the appropriate data. The two grids do not have to!      be the same resolution. Area-averaging will work for full => partial!      grid but obviously will not work for partial => full grid.!! o 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 as!!   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,maxovr),j_ovr(io,jo,maxovr)) * w_ovr(io,jo,maxovr)!! o Error checks:!!   Overlap weights of input cells sum to 1 for each output cell.!   Global sum of dummy field is conserved for input => output area-average.! ! Author: Gordon Bonan! !-----------------------------------------------------------------------! ------------------------ arguments --------------------------------    integer , intent(in)    :: io                     !output grid longitude index    integer , intent(in)    :: jo                     !output grid latitude index     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 lon points at each lat    real(r8), intent(in)    :: lon_i(nlon_i+1,nlat_i) !input grid: longitude, west edge (degrees)    real(r8), intent(in)    :: lon_i_offset(nlon_i+1,nlat_i) !input grid : cell lons, west edge (deg)    real(r8), intent(in)    :: lat_i(nlat_i+1)        !input grid: latitude, south edge (degrees)    real(r8), intent(in)    :: area_i(nlon_i,nlat_i)  !input grid: cell area    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 lon points at each lat    real(r8), intent(in)    :: lon_o(nlon_o+1,nlat_o) !output grid: longitude, west edge  (degrees)    real(r8), intent(in)    :: lat_o(nlat_o+1)        !output grid: latitude, south edge (degrees)    real(r8), intent(in)    :: area_o                 !output grid: cell area    real(r8), intent(in)    :: fland_o                !output grid: fraction that is land    integer , intent(out)   :: novr_i2o               !number of overlapping input cells    integer , intent(out)   :: iovr_i2o(maxovr)        !lon index of overlap input cell    integer , intent(out)   :: jovr_i2o(maxovr)        !lat index of overlap input cell    real(r8), intent(out)   :: wovr_i2o(maxovr)        !weight    of overlap input cell! --------------------------------------------------------------------! ------------------------ local variables ---------------------------    real(r8) :: relerr = 0.00001       !relative error for error checks    integer  :: ii                     !input  grid longitude loop index    integer  :: ji                     !input  grid latitude  loop index    integer  :: n                      !overlap index    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! --------------------------------------------------------------------! -----------------------------------------------------------------