areamod.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

!           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)
!