#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.