areamod.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <preproc.h>

module areaMod

!----------------------------------------------------------------------- 
! 
! Purpose: 
! area averaging routines
!
! Method: 
! This subroutine is used in conjunction with areaave.F for area-average
!
!-----------------------------------------------------------------------
! $Id: areaMod.F90,v 1.9.2.6 2002/05/13 18:00:01 erik Exp $
!-----------------------------------------------------------------------

  use precision
  use clm_varcon, only : re
  use shr_const_mod, only : SHR_CONST_PI
  implicit none

  integer, parameter :: maxovr = 100000

  INTERFACE celledge
     MODULE procedure celledge_regional
     MODULE procedure celledge_global
  END INTERFACE

  INTERFACE cellarea
     MODULE procedure cellarea_regional
     MODULE procedure cellarea_global
  END INTERFACE

!=======================================================================
contains
!=======================================================================

  subroutine areaini (nlon_i , nlat_i, numlon_i, lon_i, lat_i, area_i, mask_i , &
                      nlon_o , nlat_o, numlon_o, lon_o, lat_o, area_o, fland_o, &
                      mx_ovr , 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  - initializes indices and weights for area-averaging from 
!               input grid to output grid
!    areamap  - called by areaini: finds indices and weights
!    areaovr  - called by areamap: finds if cells overlap and area of overlap
!    areaave  - does area-averaging from input grid to output grid
! 
! 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,mx_ovr),j_ovr(io,jo,mx_ovr)) * w_ovr(io,jo,mx_ovr)
!
! 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) :: 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(inout) :: lon_i(nlon_i+1,nlat_i)    !input grid: longitude, west edge (degrees)
    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(nlon_o,nlat_o)        !output grid: cell area
    real(r8), intent(in) :: fland_o(nlon_o,nlat_o)       !output grid: fraction that is land
    integer , intent(in) :: mx_ovr                       !maximum number of overlapping cells 
    integer , intent(out):: novr_i2o(nlon_o,nlat_o)      !number of overlapping input cells
    integer , intent(out):: iovr_i2o(nlon_o,nlat_o,mx_ovr)!lon index of overlap input cell
    integer , intent(out):: jovr_i2o(nlon_o,nlat_o,mx_ovr)!lat index of overlap input cell
    real(r8), intent(out):: wovr_i2o(nlon_o,nlat_o,mx_ovr)!weight    of overlap input cell
! --------------------------------------------------------------------

! ------------------------ local variables ---------------------------
    real(r8),allocatable :: fld_o(:,:) !output grid: dummy field
    real(r8),allocatable :: fld_i(:,:) !input grid: dummy field
    real(r8) :: sum_fldo               !global sum of dummy output field
    real(r8) :: sum_fldi               !global sum of dummy input field
    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  :: io                     !output grid longitude loop index
    integer  :: jo                     !output grid latitude  loop 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
! --------------------------------------------------------------------

! Dynamically allocate memory

    allocate (fld_o(nlon_o,nlat_o))
    allocate (fld_i(nlon_i,nlat_i))

! -----------------------------------------------------------------
! Get indices and weights for mapping from input grid to output grid
! -----------------------------------------------------------------

    call 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   , &
                  novr_i2o , iovr_i2o , jovr_i2o , wovr_i2o , &
                  fland_o  , area_o   ) 

! -----------------------------------------------------------------
! Error check: global sum fld_o = global sum fld_i. 
! This true only if both grids span the same domain. 
! -----------------------------------------------------------------

    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
       write (6,*) 'AREAINI warning: conservation check 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

! make dummy input field and sum globally

    sum_fldi = 0.
    do ji = 1, nlat_i      
       do ii = 1, numlon_i(ji)
          fld_i(ii,ji) = ((ji-1)*nlon_i + ii) * mask_i(ii,ji)
          sum_fldi = sum_fldi + area_i(ii,ji)*fld_i(ii,ji)
       end do
    end do

! area-average output field from input field

    call areaave (nlat_i   , nlon_i   , numlon_i , fld_i , &
                  nlat_o   , nlon_o   , numlon_o , fld_o , &
                  iovr_i2o , jovr_i2o , wovr_i2o , mx_ovr )

! global sum of output field -- must multiply by fraction of output
! grid that is land as determined by input grid

    sum_fldo = 0.
    do jo = 1, nlat_o
       do io = 1, numlon_o(jo)
          sum_fldo = sum_fldo + area_o(io,jo)*fld_o(io,jo) * fland_o(io,jo)
       end do
    end do

! check for conservation

    if ( abs(sum_fldo/sum_fldi-1.) > relerr ) then
       write (6,*) 'AREAINI error: input field not conserved'
       write (6,'(a30,e20.10)') 'global sum output field = ',sum_fldo
       write (6,'(a30,e20.10)') 'global sum input  field = ',sum_fldi
       call endrun
    end if

    deallocate (fld_o)
    deallocate (fld_i)

    return
  end subroutine areaini

!=======================================================================

  subroutine areaave (nlat_i , nlon_i , numlon_i, fld_i , &
                      nlat_o , nlon_o , numlon_o, fld_o , &
                      i_ovr  , j_ovr  , w_ovr   , nmax  )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! area averaging of field from input to output grids
!
! Method: 
! 
! Author: Gordon Bonan
! 
!-----------------------------------------------------------------------

! ------------------------ arguments ---------------------------------
    integer ,intent(in) :: nlat_i                    !input grid : number of latitude points      
    integer ,intent(in) :: nlon_i                    !input grid : max number longitude points    
    integer ,intent(in) :: numlon_i(nlat_i)          !input grid : number of lon points at each lat
    real(r8),intent(in) :: fld_i(nlon_i,nlat_i)      !input grid : field
    integer ,intent(in) :: nlat_o                    !output grid: number of latitude points      
    integer ,intent(in) :: nlon_o                    !output grid: max number of longitude points 
    integer ,intent(in) :: numlon_o(nlat_o)          !output grid: number of lon points at each lat
    real(r8),intent(out):: fld_o(nlon_o,nlat_o)      !field for output grid
    integer ,intent(in) :: nmax                      !input grid : max number of overlapping cells
    integer ,intent(in) :: i_ovr(nlon_o,nlat_o,nmax) !lon index, overlapping input cell
    integer ,intent(in) :: j_ovr(nlon_o,nlat_o,nmax) !lat index, overlapping input cell
    real(r8),intent(in) :: w_ovr(nlon_o,nlat_o,nmax) !overlap weights for input cells
! --------------------------------------------------------------------

! ------------------- local variables -----------------------------
    integer jo                !latitude index for output grid
    integer io                !longitude index for output grid
    integer ji                !latitude index for input grid
    integer ii                !longitude index for input grid
    integer n                 !overlapping cell index
! -----------------------------------------------------------------

! initialize field on output grid to zero everywhere

!$OMP PARALLEL DO PRIVATE (jo,io)
    do jo = 1, nlat_o
       do io = 1, numlon_o(jo)
          fld_o(io,jo) = 0.