mksoitex.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

subroutine mksoitex (fsoitex, ndiag, pctgla_o, sand_o, clay_o)

!----------------------------------------------------------------------- 
! 
! Purpose: 
! make %sand and %clay from IGBP soil data, which includes
! igbp soil 'mapunits' and their corresponding textures
! 
! Method: 
! 
! Authors: Gordon Bonan and Sam Levis
! 
!-----------------------------------------------------------------------
! $Id: mksoitex.F90,v 1.4.6.2 2002/04/27 15:38:55 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) :: fsoitex        !soil texture dataset file name
  integer , intent(in) :: ndiag                  !unit # for diagnostic output
  real(r8), intent(in) :: pctgla_o(lsmlon,lsmlat)      !% glacier (output grid)
  real(r8), intent(out):: sand_o(lsmlon,lsmlat,nlevsoi) !% sand (output grid)
  real(r8), intent(out):: clay_o(lsmlon,lsmlat,nlevsoi) !% clay (output grid)
! -----------------------------------------------------------------

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

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

  integer :: nlay                        !number of soil layers
  integer :: mapunitmax                  !maximum value of igbp soil mapunits
  integer :: mapunittemp                 !temporary igbp soil mapunit
  integer :: ii                          !longitude index for IGBP grid
  integer :: ji                          !latitude  index for IGBP grid
  integer :: io                          !longitude index for land grid
  integer :: jo                          !latitude  index for land grid
  integer :: l,m,n                       !loop indices
  integer :: miss = 99999                !missing data indicator
  integer :: k                           !igbp soil mapunit index
  integer, parameter :: num=2            !get 1st and 2nd largest areas of overlap
  integer :: wsti(num)                   !index to 1st and 2nd largest values in wst
  integer, parameter :: nwstmax = 10000  !maximum size of overlap weights, by soil mapunit
  real(r8):: wst(0:nwstmax)              !overlap weights, by soil mapunit
  integer, parameter :: nlsm=4           !number of soil textures (sand, silt, clay)
  character(len=38) :: soil(0:nlsm)      !name of each soil texture
  character(len=38) :: typ               !soil texture based on %sand, silt, clay
  real(r8) :: gast_i(0:nlsm)             !input grid : global area, by texture type
  real(r8) :: gast_o(0:nlsm)             !output grid: global area, by texture type

  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) :: mask_o                         !output grid: mask (0, 1)
  real(r8) :: offset                         !used to shift x-grid 360 degrees
  real(r8) :: mapunit_o(lsmlon,lsmlat)       !not needed except as diagnostic

  real(r8) :: edge_i(4)                      !input grid: N,E,S,W edges (degrees)
  real(r8), allocatable :: sand_i(:,:)       !input grid: percent sand
  real(r8), allocatable :: clay_i(:,:)       !input grid: percent clay
  real(r8), allocatable :: landmask_i(:,:)   !input grid: land=1, no_soil_data=0
  real(r8), allocatable :: mapunit_i(:,:)    !input grid: igbp soil mapunits
  real(r8), allocatable :: latixy_i(:,:)     !input grid: latitude (degrees)
  real(r8), allocatable :: longxy_i(:,:)     !input grid: longitude (degrees)
  integer , allocatable :: numlon_i(:)       !input grid: # longitude pts by lat
  real(r8), allocatable :: lon_i(:,:)        !input grid: longitude, W edge (deg)
  real(r8), allocatable :: lon_i_offset(:,:) !input grid: offset longitude, west edge (degrees)
  real(r8), allocatable :: lat_i(:)          !input grid: latitude,  S edge (deg)
  real(r8), allocatable :: area_i(:,:)       !input grid: cell area
  real(r8), allocatable :: mask_i(:,:)   !input grid: mask (0, 1)

  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 %sand and %clay .....'
  call shr_sys_flush(6)

! -----------------------------------------------------------------
! Define the model surface types: 0 to nlsm
! -----------------------------------------------------------------

  soil(0) = 'no soil: ocean, glacier, lake, no data'
  soil(1) = 'clays                                 '
  soil(2) = 'sands                                 '
  soil(3) = 'loams                                 '
  soil(4) = 'silts                                 '

! -----------------------------------------------------------------
! Read input data
! -----------------------------------------------------------------

! Obtain input grid info

  call getfil (fsoitex, 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, 'number_of_layers', dimid)
  call wrap_inq_dimlen (ncid, dimid, nlay)

  call wrap_inq_dimid  (ncid, 'max_value_mapunit', dimid)
  call wrap_inq_dimlen (ncid, dimid, mapunitmax)

  !NOTE: wst must not be allocatable if it is declared private in a thread                                         
  if (nwstmax < mapunitmax) then
     write(6,*)'MKSOITEX: parameter nwstmax must be increased to ',mapunitmax
     call endrun
  endif

  allocate (sand_i(mapunitmax,nlay), stat=ier)
  if (ier/=0) call endrun
  allocate (clay_i(mapunitmax,nlay), stat=ier)
  if (ier/=0) call endrun     
  allocate (landmask_i(nlon_i,nlat_i), stat=ier)
  if (ier/=0) call endrun
  allocate (mapunit_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 (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  

  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, 'MAPUNITS', varid)
  call wrap_get_var_realx (ncid, varid, mapunit_i)

  call wrap_inq_varid (ncid, 'PCT_SAND', varid)
  call wrap_get_var_realx (ncid, varid, sand_i)

  call wrap_inq_varid (ncid, 'PCT_CLAY', varid)
  call wrap_get_var_realx (ncid, varid, clay_i)

  call wrap_close(ncid)

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

! Determine input grid cell edges 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
  
! 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,l,k,mask_o,novr_i2o,iovr_i2o,jovr_i2o,wovr_i2o,fld_i,wst,wsti)
  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    , &