mkpft.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

subroutine mkpft (fpft, ndiag,  noveg,  pctlnd_o, pft, pctpft)

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Make PFT data for vegetated patches (1 to maxpatch_pft)
! 
! Method: 
! This dataset consists of the %cover of the [numpft]+1 PFTs used by
! the model. The %cover pertains to the "vegetated" portion of the
! grid cell and sums to 100. The real portion of each grid cell
! covered by each PFT is the PFT cover times the fraction of the
! grid cell that is land. This is the quantity preserved when
! area-averaging from the input (1/2 degree) grid to the models grid.
! 
! Author: Mariana Vertenstein
! 
!-----------------------------------------------------------------------
! $Id: mkpft.F90,v 1.3.10.2 2002/04/27 15:38:54 erik Exp $
!-----------------------------------------------------------------------

  use precision
  use clm_varpar    !parameters
  use clm_varsur    !surface variables
  use clm_varctl    !run control variables
  use fileutils, only : getfil
  use areaMod      !area averaging routines 
  use shr_sys_mod, only : shr_sys_flush 
  implicit none

! ------------------------ arguments ------------------------------
  character(len=*), intent(in) :: fpft                        !input pft dataset file name
  integer , intent(in) :: ndiag                               !unit number for diagnostic output
  integer , intent(in) :: noveg                               !PFT number for no vegetation
  real(r8), intent(out):: pctlnd_o(lsmlon,lsmlat)             !output grid: % land per gridcell
  integer , intent(out):: pft(lsmlon,lsmlat,maxpatch_pft)     !output grid PFT (0 to numpft)
  real(r8), intent(out):: pctpft(lsmlon,lsmlat,maxpatch_pft)  !output grid PFT cover (% of vegetated area)
! -----------------------------------------------------------------

! ------------------------ local variables ------------------------
  character(len=256) locfn                    !local dataset file name

  integer :: nlon_i                 !input grid : longitude points (read in)
  integer :: nlat_i                 !input grid : latitude  points (read in)
  integer :: ncid,dimid,varid       !input netCDF id's
  integer :: ier                    !error status   

  integer  :: miss = 99999                    !missing data indicator
  real(r8) :: pctpft_o(lsmlon,lsmlat,0:numpft)!PFT percent on output grid
  real(r8) :: wst(0:numpft)                   !as pft_o at specific io, jo
  integer  :: wsti(maxpatch_pft)              !ranked indices of largest values in wst
  real(r8) :: wst_sum                         !sum of %pft
  real(r8) :: sumpct                          !sum of %pft over maxpatch_pft
  real(r8) :: diff                            !the difference (wst_sum - sumpct)
  real(r8) :: gpft_o(0:numpft)                !output grid: global area pfts
  real(r8) :: garea_o                         !output grid: global area
  real(r8) :: gpft_i(0:numpft)                !input grid: global area pfts
  real(r8) :: garea_i                         !input grid: global area
  integer  :: ii                              !longitude index for input grid
  integer  :: ji                              !latitude  index for input grid
  integer  :: io                              !longitude index for model grid
  integer  :: jo                              !latitude  index for model grid
  integer  :: k,n,m                           !indices

  integer numpft_i                            !number of plant types on input dataset
  real(r8) :: edge_i(4)                       !input grid: N,E,S,W edges (degrees)
  real(r8), allocatable :: pctpft_i(:,:,:)    !input grid: PFT percent 
  real(r8), allocatable :: landmask_i(:,:)    !input grid: fraction land (not ocn) per land gridcell
  real(r8), allocatable :: mask_i(:,:)        !input grid: mask (0, 1)
  real(r8), allocatable :: latixy_i(:,:)      !input grid: latitude (degrees)
  real(r8), allocatable :: longxy_i(:,:)      !input grid: longitude (degrees)
  integer , allocatable :: numlon_i(:)        !input grid: number longitude points by lat
  real(r8), allocatable :: lon_i(:,:)         !input grid: longitude, west edge (degrees)
  real(r8), allocatable :: lon_i_offset(:,:)!input grid: offset longitude, west edge (degrees)
  real(r8), allocatable :: lat_i(:)           !input grid: latitude, south edge (degrees)
  real(r8), allocatable :: area_i(:,:)        !input grid: cell area

  real(r8) :: mask_o                          !output grid: mask (0, 1)
  integer  :: novr_i2o                        !number of overlapping input cells
  integer  :: iovr_i2o(maxovr)                !lon index of overlap input cell
  integer  :: jovr_i2o(maxovr)                !lat index of overlap input cell
  real(r8) :: wovr_i2o(maxovr)                !weight    of overlap input cell
  real(r8) :: offset                          !used to shift x-grid 360 degrees

  real(r8) :: fld_o(lsmlon,lsmlat)            !output grid: dummy field 
  real(r8) :: 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                !max error: sum overlap weights ne 1

  character(len=35)  veg(0:numpft)            !vegetation types
  data veg( 0) /'not vegetated'                      /
  data veg( 1) /'needleleaf evergreen temperate tree'/
  data veg( 2) /'needleleaf evergreen boreal tree'   /
  data veg( 3) /'needleleaf deciduous boreal tree'   /
  data veg( 4) /'broadleaf evergreen tropical tree'  /
  data veg( 5) /'broadleaf evergreen temperate tree' /
  data veg( 6) /'broadleaf deciduous tropical tree'  /
  data veg( 7) /'broadleaf deciduous temperate tree' /
  data veg( 8) /'broadleaf deciduous boreal tree'    /
  data veg( 9) /'broadleaf evergreen shrub'          /
  data veg(10) /'broadleaf deciduous temperate shrub'/
  data veg(11) /'broadleaf deciduous boreal shrub'   /
  data veg(12) /'c3 arctic grass'                    /
  data veg(13) /'c3 non-arctic grass'                /
  data veg(14) /'c4 grass'                           / 
  data veg(15) /'corn'                               /
  data veg(16) /'wheat'                              /
! -----------------------------------------------------------------

  write (6,*) 'Attempting to make PFTs .....'
  call shr_sys_flush(6)

! -----------------------------------------------------------------
! Read input PFT file
! -----------------------------------------------------------------

! Obtain input grid info

  call getfil (fpft, locfn, 0)
  call wrap_open(locfn, 0, ncid)

  call wrap_inq_dimid  (ncid, 'lon', dimid)
  call wrap_inq_dimlen (ncid, dimid, nlon_i)

  call wrap_inq_dimid  (ncid, 'lat', dimid)
  call wrap_inq_dimlen (ncid, dimid, nlat_i)

  call wrap_inq_dimid  (ncid, 'pft', dimid)
  call wrap_inq_dimlen (ncid, dimid, numpft_i)

  if (numpft_i .ne. numpft+1) then
     write(6,*)'MKPFT: parameter numpft+1= ',numpft+1, &
          'does not equal input dataset numpft= ',numpft_i
     call endrun
  endif

  allocate (latixy_i(nlon_i,nlat_i), stat=ier)   
  if (ier/=0) call endrun
  allocate (longxy_i(nlon_i,nlat_i), stat=ier)   
  if (ier/=0) call endrun
  allocate (numlon_i(nlat_i), stat=ier)
  if (ier/=0) call endrun
  allocate (lon_i(nlon_i+1,nlat_i), stat=ier)
  if (ier/=0) call endrun
  allocate (lon_i_offset(nlon_i+1,nlat_i), stat=ier)
  if (ier/=0) call endrun
  allocate (lat_i(nlat_i+1), stat=ier)        
  if (ier/=0) call endrun
  allocate (landmask_i(nlon_i,nlat_i), stat=ier)
  if (ier/=0) call endrun
  allocate (mask_i(nlon_i,nlat_i), stat=ier)     
  if (ier/=0) call endrun
  allocate (area_i(nlon_i,nlat_i), stat=ier)  
  if (ier/=0) call endrun
  allocate (pctpft_i(nlon_i,nlat_i,0:numpft), stat=ier)     
  if (ier/=0) call endrun

  call wrap_inq_varid (ncid, 'LATIXY', varid)
  call wrap_get_var_realx (ncid, varid, latixy_i)

  call wrap_inq_varid (ncid, 'LONGXY', varid)
  call wrap_get_var_realx (ncid, varid, longxy_i)

  call wrap_inq_varid (ncid, 'EDGEN', varid)
  call wrap_get_var_realx (ncid, varid, edge_i(1))

  call wrap_inq_varid (ncid, 'EDGEE', varid)
  call wrap_get_var_realx (ncid, varid, edge_i(2))

  call wrap_inq_varid (ncid, 'EDGES', varid)
  call wrap_get_var_realx (ncid, varid, edge_i(3))

  call wrap_inq_varid (ncid, 'EDGEW', varid)
  call wrap_get_var_realx (ncid, varid, edge_i(4))

! Obtain input data

  call wrap_inq_varid (ncid, 'LANDMASK', varid)
  call wrap_get_var_realx (ncid, varid, landmask_i)

  call wrap_inq_varid (ncid, 'PCT_PFT', varid)
  call wrap_get_var_realx (ncid, varid, pctpft_i)

  call wrap_close(ncid)

! -----------------------------------------------------------------
! Map data from input grid to land model grid. Get:
! -----------------------------------------------------------------

  numlon_i(:) = nlon_i

  call celledge (nlat_i    , nlon_i    , numlon_i  , longxy_i  ,  &
                 latixy_i  , edge_i(1) , edge_i(2) , edge_i(3) ,  &
                 edge_i(4) , lat_i     , lon_i     )

  call cellarea (nlat_i    , nlon_i    , numlon_i  , lat_i     ,  &
                 lon_i     , edge_i(1) , edge_i(2) , edge_i(3) ,  &
                 edge_i(4) , area_i    )

  do ji = 1, nlat_i
     do ii = 1, numlon_i(ji)
        mask_i(ii,ji) = 1.
     end do
  end do

! Shift x-grid to locate periodic grid intersections. This
! assumes that all lon_i(1,j) have the same value for all
! latitudes j and that the same holds for lon_o(1,j)

  if (lon_i(1,1) < lonw(1,1)) then
     offset = 360.0
  else
     offset = -360.0
  end if
  
  do ji = 1, nlat_i
     do ii = 1, numlon_i(ji) + 1
        lon_i_offset(ii,ji) = lon_i(ii,ji) + offset
     end do
  end do

! Process each cell on land model grid
! novr_i2o - number of input grid cells that overlap each land grid cell
! iovr_i2o - longitude index of overlapping input grid cell
! jovr_i2o - latitude  index of overlapping input grid cell
! wovr_i2o - fraction of land grid cell overlapped by input grid cell

!$OMP PARALLEL DO PRIVATE (io,jo,ii,ji,n,m,mask_o,novr_i2o,iovr_i2o,jovr_i2o,wovr_i2o,fld_i, &
!$OMP &  wst, wst_sum, sumpct, wsti, diff)
  do jo = 1, lsmlat
     do io = 1, numlon(jo)

! Determine areas of overlap and indices

        mask_o = 1.

        call areaini_point (io        , jo          , nlon_i  , nlat_i  , numlon_i , &
                           lon_i      , lon_i_offset, lat_i   , area_i  , mask_i   , &
                           lsmlon     , lsmlat      , numlon  , lonw    , lats     , &
                           area(io,jo), mask_o      , novr_i2o, iovr_i2o, jovr_i2o , &
                           wovr_i2o)                             

! Make percent of grid cell that is land [pctlnd_o] and map data 
! from input grid to land model grid again, this time with the real mask_i

        pctlnd_o(io,jo) = 0.   
        do n = 1, novr_i2o  !overlap cell index
           ii = iovr_i2o(n) !lon index (input grid) of overlap cell
           ji = jovr_i2o(n) !lat index (input grid) of overlap cell