mklai.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

subroutine mklai (flai, ndiag, pft, mlai, msai, mhgtt, mhgtb)

!----------------------------------------------------------------------- 
! 
! Purpose: 
! make LAI/SAI/height data for vegetated patches (1 to maxpatch_pft)
! 
! Method: 
! Portions of this code could be moved out of the month loop
! for improved efficiency
! 
! Author: Sam Levis
! 
!-----------------------------------------------------------------------
! $Id: mklai.F90,v 1.2.10.2 2002/04/27 15:38:54 erik Exp $
!-----------------------------------------------------------------------

  use precision
  use clm_varpar    !parameters
  use clm_varctl    !run control variables
  use clm_varsur    !surface 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) :: flai                         !input lai-sai-hgt dataset
  integer , intent(in) :: ndiag                                !unit number for diagnostic output
  integer , intent(in) :: pft(lsmlon,lsmlat,maxpatch_pft)      !PFT (0 to numpft)
  real(r8), intent(out):: mlai(lsmlon,lsmlat,maxpatch_pft,12)  !monthly lai
  real(r8), intent(out):: msai(lsmlon,lsmlat,maxpatch_pft,12)  !monthly sai
  real(r8), intent(out):: mhgtt(lsmlon,lsmlat,maxpatch_pft,12) !monthly height (top)
  real(r8), intent(out):: mhgtb(lsmlon,lsmlat,maxpatch_pft,12) !monthly height (bottom)
! -----------------------------------------------------------------

! ------------------------ 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 :: beg4d(4),len4d(4)      !netCDF variable edges
  integer :: ntim                   !number of input time samples
  integer :: ier                    !error status   

  real(r8) :: glai_o(0:numpft)      !output grid: global area pfts
  real(r8) :: gsai_o(0:numpft)      !output grid: global area pfts
  real(r8) :: ghgtt_o(0:numpft)     !output grid: global area pfts
  real(r8) :: ghgtb_o(0:numpft)     !output grid: global area pfts
  real(r8) :: garea_o               !output grid: global area
  real(r8) :: glai_i(0:numpft)      !input grid: global area pfts
  real(r8) :: gsai_i(0:numpft)      !input grid: global area pfts
  real(r8) :: ghgtt_i(0:numpft)     !input grid: global area pfts
  real(r8) :: ghgtb_i(0:numpft)     !input grid: global area pfts
  real(r8) :: garea_i               !input grid: global area

  integer  :: ii                    !longitude index for input grid
  integer  :: io                    !longitude index for model grid
  integer  :: ji                    !latitude  index for input grid
  integer  :: jo                    !latitude  index for model grid
  integer  :: k,l,m,n               !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 :: mlai_i(:,:,:)      !monthly lai in
  real(r8), allocatable :: msai_i(:,:,:)      !monthly sai in
  real(r8), allocatable :: mhgtt_i(:,:,:)     !monthly height (top) in
  real(r8), allocatable :: mhgtb_i(:,:,:)     !monthly height (bottom) in
  real(r8), allocatable :: landmask_i(:,:)    !input grid: fraction land
  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), allocatable :: mask_i(:,:)        !input grid: mask (0, 1)

  real(r8) :: mlai_o(lsmlon,lsmlat,0:numpft)  !monthly lai out
  real(r8) :: msai_o(lsmlon,lsmlat,0:numpft)  !monthly sai out
  real(r8) :: mhgtt_o(lsmlon,lsmlat,0:numpft) !monthly height (top) out
  real(r8) :: mhgtb_o(lsmlon,lsmlat,0:numpft) !monthly height (bottom) out

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

  write (6,*) 'Attempting to make LAIs/SAIs/heights .....'
  call shr_sys_flush(6)

! -----------------------------------------------------------------
! Determine input grid info
! -----------------------------------------------------------------

  call getfil (flai, 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,*)'MKLAI: parameter numpft+1= ',numpft+1, &
          'does not equal input dataset numpft= ',numpft_i
     call endrun
  endif

  call wrap_inq_dimid  (ncid, 'time', dimid)
  call wrap_inq_dimlen (ncid, dimid, ntim)
  if (ntim .ne. 12) then
     write(6,*)'MKLAI: must have 12 time samples on input data'
     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 (area_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 (mlai_i(nlon_i,nlat_i,0:numpft), stat=ier) 
  if (ier/=0) call endrun
  allocate (msai_i(nlon_i,nlat_i,0:numpft), stat=ier) 
  if (ier/=0) call endrun
  allocate (mhgtt_i(nlon_i,nlat_i,0:numpft), stat=ier)
  if (ier/=0) call endrun
  allocate (mhgtb_i(nlon_i,nlat_i,0:numpft), stat=ier)
  if (ier/=0) call endrun
  allocate (landmask_i(nlon_i,nlat_i), 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 time independent input data

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

! -----------------------------------------------------------------
! Determine input grid cell and cell areas
! -----------------------------------------------------------------


  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

! -----------------------------------------------------------------
! Loop over input months
! -----------------------------------------------------------------

  do m = 1, ntim

! Get input data for the month

     beg4d(1) = 1 ; len4d(1) = nlon_i
     beg4d(2) = 1 ; len4d(2) = nlat_i
     beg4d(3) = 1 ; len4d(3) = numpft_i
     beg4d(4) = m ; len4d(4) = 1

     call wrap_inq_varid (ncid, 'MONTHLY_LAI', varid)