!           lon_o(io,jo)      lon_o(io+1,jo)!!              |                   |!              --------------------- lat_o(jo+1)!              |                   |!              |                   |!    xxxxxxxxxxxxxxx lat_i(ji+1)   |!    x         |   x               |!    x  input  |   x   output      |!    x  cell   |   x    cell       |!    x  ii,ji  |   x   io,jo       |!    x         |   x               |!    x         ----x---------------- lat_o(jo  )!    x             x!    xxxxxxxxxxxxxxx lat_i(ji  )!    x             x! lon_i(ii,ji) lon_i(ii+1,ji)!!! The above diagram assumes both grids are oriented South to North. Other! combinations of North to South and South to North grids are possible:!!     Input Grid    Output Grid!     -------------------------! (1)   S to N        S to N! (2)   N to S        N to S! (3)   S to N        N to S! (4)   N to S        S to N!! The code has been modified to allow for North to South grids. Verification! that these changes work are: !    o (1) and (4) give same results for output grid!    o (2) and (3) give same results for output grid!    o (2) and (4) give same results for output grid when output grid inverted!! WARNING: this code does not vectorize but is only called during start-up! ! 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 of lon points for lat    real(r8), intent(in) :: lon_i(nlon_i+1,nlat_i) !input grid : cell longitude, W edge (deg)    real(r8), intent(in) :: lat_i(nlat_i+1)        !input grid : cell latitude, S edge (deg)    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 of lon points for lat    real(r8), intent(in) :: lon_o(nlon_o+1,nlat_o) !output grid: cell longitude, W edge (deg)    real(r8), intent(in) :: lat_o(nlat_o+1)        !output grid: cell latitude, S edge (deg)    integer , intent(in) :: mx_ovr                 !maximum number of overlapping input cells    integer , intent(inout) :: n_ovr(nlon_o,nlat_o       ) !number of overlapping input cells    integer , intent(inout) :: i_ovr(nlon_o,nlat_o,mx_ovr) !lon index, overlapping input cell    integer , intent(inout) :: j_ovr(nlon_o,nlat_o,mx_ovr) !lat index, overlapping input cell    real(r8), intent(inout) :: w_ovr(nlon_o,nlat_o,mx_ovr) !area of overlap for input cells! --------------------------------------------------------------------! ------------------------ local variables ---------------------------    integer io             !output grid longitude loop index    integer jo             !output grid latitude  loop index    integer indexo1        !output grid lat. index according to orientn    integer indexo2        !output grid lat. index according to orientn    integer indexo3        !output grid lat. index according to orientn    integer ii             !input  grid longitude loop index    integer ji             !input  grid latitude  loop index    integer indexi1        !input grid lat. index according to orientn    integer indexi2        !input grid lat. index according to orientn    integer indexi3        !input grid lat. index according to orientn    real(r8) lonw          !west longitudes of overlap     real(r8) lone          !east longitudes of overlap     real(r8) dx            !difference in longitudes    real(r8) lats          !south latitudes of overlap    real(r8) latn          !north latitudes of overlap    real(r8) dy            !difference in latitudes     real(r8) deg2rad       !pi/180    real(r8) a_ovr         !area of overlap! --------------------------------------------------------------------    deg2rad = SHR_CONST_PI / 180.    do jo = 1, nlat_o! choose the right index according to the orientation of the data       if (lat_o(nlat_o+1) > lat_o(1)) then          indexo1 = jo+1        !south to north along the edges           indexo2 = jo          !south to north along the edges          indexo3 = jo          !south to north at the center of cell       else          indexo1 = nlat_o+1-jo !north to south along the edges          indexo2 = nlat_o+2-jo !north to south along the edges          indexo3 = nlat_o+1-jo !north to south at the center of cell       end if       do io = 1, numlon_o(indexo3)! loop through all input grid cells to find overlap with output grid          do ji = 1, nlat_i                            ! choose the right index according to the orientation of the data             if (lat_i(nlat_i+1) > lat_i(1)) then                indexi1 = ji          !south to north along the edges                indexi2 = ji+1        !south to north along the edges                indexi3 = ji          !south to north at the center of cell             else                indexi1 = nlat_i+2-ji !north to south along the edges                indexi2 = nlat_i+1-ji !north to south along the edges                indexi3 = nlat_i+1-ji !north to south at the center of cell             end if! lat okay             if ( lat_i(indexi1)<lat_o(indexo1) .and. &                  lat_i(indexi2)>lat_o(indexo2) ) then                do ii = 1, numlon_i(indexi3)! lon okay                   if (lon_i(ii,indexi3)<lon_o(io+1,indexo3) .and. &                       lon_i(ii+1,indexi3)>lon_o(io,indexo3)) then! increment number of overlapping cells. make sure 0 < n_ovr < mx_ovr                      n_ovr(io,indexo3) = n_ovr(io,indexo3) + 1                      if (n_ovr(io,indexo3) > mx_ovr) then                         write (6,*) 'AREAOVR error: n_ovr= ', &                              n_ovr(io,indexo3),' exceeded mx_ovr = ', &                              mx_ovr,' for output lon,lat = ',io,indexo3                         call endrun                      end if! determine area of overlap                      lone = min(lon_o(io+1,indexo3),lon_i(ii+1,indexi3))*deg2rad !e edge                      lonw = max(lon_o(io  ,indexo3),lon_i(ii  ,indexi3))*deg2rad !w edge                      dx = max(0.0_r8,(lone-lonw))                      latn = min(lat_o(indexo1),lat_i(indexi2))*deg2rad !n edge                      lats = max(lat_o(indexo2),lat_i(indexi1))*deg2rad !s edge                      dy = max(0.0_r8,(sin(latn)-sin(lats)))                      a_ovr = dx*dy*re*re! save lat, lon indices of overlapping cell and area of overlap                      i_ovr(io,indexo3,n_ovr(io,indexo3)) = ii                      j_ovr(io,indexo3,n_ovr(io,indexo3)) = indexi3                      w_ovr(io,indexo3,n_ovr(io,indexo3)) = a_ovr                   end if                end do             end if          end do       end do    end do    return  end subroutine areaovr!=======================================================================  subroutine cellarea_regional (nlat  ,nlon  ,numlon, lats ,lonw , &                                edgen, edgee, edges , edgew, area)!----------------------------------------------------------------------- ! ! Purpose: ! Area of grid cells (square kilometers) - regional grid! (can become global as special case)!! 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) :: edgen              !northern 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(in) :: edgee              !eastern edge of grid (degrees)    real(r8), intent(in) :: lats(nlat+1)       !grid cell latitude, southern edge (degrees)    real(r8), intent(in) :: lonw(nlon+1,nlat)  !grid cell longitude, western edge (degrees)    real(r8), intent(out):: area(nlon,nlat)    !cell area (km**2)! --------------------------------------------------------------------! ------------------------ local variables ---------------------------    integer i,j                 !indices    real(r8) deg2rad            !pi/180    real(r8) global             !summed area    real(r8) dx                 !cell width: E-W    real(r8) dy                 !cell width: N-S    real(r8) error              !true area for error check! --------------------------------------------------------------------    deg2rad = (SHR_CONST_PI) / 180.    global = 0.    do j = 1, nlat       do i = 1, numlon(j)          dx = (lonw(i+1,j) - lonw(i,j)) * deg2rad          if (lats(j+1) > lats(j)) then        !South to North grid             dy = sin(lats(j+1)*deg2rad) - sin(lats(j)*deg2rad)          else                                 !North to South grid             dy = sin(lats(j)*deg2rad) - sin(lats(j+1)*deg2rad)          end if          area(i,j) = dx*dy*re*re          global = global + area(i,j)       end do    end do! make sure total area from grid cells is same as area of grid! as defined by its edges    dx = (edgee - edgew) * deg2rad    dy = sin(edgen*deg2rad) - sin(edges*deg2rad)    error =  dx*dy*re*re    if (abs(global-error)/error > 0.00001) then       write (6,*) 'CELLAREA error: correct area is ',error, &            ' but summed area of grid cells is ',global       call endrun    end if    return  end subroutine cellarea_regional!=======================================================================  subroutine cellarea_global (nlat , nlon, numlon, lats, lonw, area)!----------------------------------------------------------------------- ! ! Purpose: ! Area of grid cells (square kilometers)- 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) :: lats(nlat+1)     !grid cell latitude, southern edge (degrees)    real(r8), intent(in) :: lonw(nlon+1,nlat)!grid cell longitude, western edge (degrees)    real(r8), intent(out):: area(nlon,nlat)  !cell area (km**2)! --------------------------------------------------------------------! ------------------------ local variables ---------------------------    integer i,j                 !indices    real(r8) deg2rad            !pi/180    real(r8) dx                 !cell width: E-W    real(r8) dy                 !cell width: N-S! --------------------------------------------------------------------! Note: assume that cam latitudes go S->N! Note: cannot make sure total area from grid cells is same as area of grid! as defined by its edges as in offline case since the edges are not all at! the same longitudes for every latitude    deg2rad = (SHR_CONST_PI) / 180.    do j = 1, nlat       do i = 1, numlon(j)          dx = (lonw(i+1,j) - lonw(i,j)) * deg2rad          dy = sin(lats(j+1)*deg2rad) - sin(lats(j)*deg2rad)  !s->n latiatudes          area(i,j) = dx*dy*re*re       end do    end do    return  end subroutine cellarea_global!=======================================================================  subroutine celledge_regional (nlat   , nlon  , numlon , longxy ,  &                                latixy , edgen , edgee  , edges  ,  &                                edgew  , lats  , lonw   )!----------------------------------------------------------------------- ! ! Purpose: ! Southern and western edges of grid cells - regional grid ! (can become global as special case)!