convert_soitex.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

program convert_soitex

  implicit none
  include 'netcdf.inc'

!-----------------------------------------------------------------
! make surface type netcdf file
!-----------------------------------------------------------------

  integer, parameter :: r8 = selected_real_kind(12)

! File specific settings

  integer, parameter :: nlon = 4320      !input grid : longitude points
  integer, parameter :: nlat = 2160      !input grid : latitude  points
  integer, parameter :: nlay = 10        !input grid : number of soil layers
  integer, parameter :: nmapunits = 4931 !input grid : # of igbp soil 'mapunits'
  integer, parameter :: mapunitmax = 6998!input grid : max value of 'mapunits'

  real(r8) :: dzsoi(10), zsoi(10)        !soil layer thickness and depth

  real(r8) :: lon(nlon)                   !longitude dimension array (1d)
  real(r8) :: lat(nlat)                   !latitude dimension array (1d) 
  real(r8) :: longxy(nlon,nlat)           !longitude dimension array (2d)
  real(r8) :: latixy(nlon,nlat)           !longitude dimension array (2d)
  real(r8) :: edge(4)                     !N,E,S,W edges of grid
  real(r8) :: dx,dy                       !grid increments

  integer  :: mu                          !current mapunit
  real(r8) :: pct_sand(mapunitmax,nlay)   !pct sand 
  real(r8) :: pct_clay(mapunitmax,nlay)   !pct clay
  real(r8) :: landmask(nlon,nlat)         !land mask
  real(r8) :: mapunit(nlon,nlat)          !global map of igbp soil mapunits
  real(r8) :: temp(nlon,nlat)             !same; used as temporary buffer

  integer :: dimlon_id                    !netCDF dimension id
  integer :: dimlat_id                    !netCDF dimension id
  integer :: dimlay_id                    !netCDF dimension id
  integer :: dimmapunits_id               !netCDF dimension id
  integer :: dimmapunitmax_id             !netCDF dimension id

  integer :: dzsoi_id                     !soil thickness by layer
  integer :: zsoi_id                      !soil depth by layer
  integer :: lon_id                       !1d longitude array id
  integer :: lat_id                       !1d latitude array id
  integer :: lay_id                       !1d layer array id
  integer :: mapunit_id                   !2d mapunits array id
  integer :: longxy_id                    !2d longitude array id
  integer :: latixy_id                    !2d latitude array id
  integer :: edgen_id                     !northern edge of grid (edge(1)) id
  integer :: edgee_id                     !eastern  edge of grid (edge(2)) id
  integer :: edges_id                     !southern edge of grid (edge(3)) id
  integer :: edgew_id                     !western  edge of grid (edge(4)) id
  integer :: pct_sand_id                  !sand id
  integer :: pct_clay_id                  !clay id
  integer :: landmask_id                  !landmask id

  integer :: i,j,k                        !indices
  integer :: ndata = 1                    !input unit
  integer :: ncid                         !netCDF file id
  integer :: dim1_id(1)                   !netCDF dimension id for 1-d variables
  integer :: dim2_id(2)                   !netCDF dimension id for 2-d variables
  integer :: status                       !status

  character(len=256) ::  filei,   fileo   !input,output filenames
  character(len=256) ::  filei1,  filei2,  filei3,  filei4,  filei5,  filei6
  character(len=256) ::  filei7,  filei8,  filei9, filei10, filei11, filei12
  character(len=256) :: filei13, filei14, filei15, filei16, filei17, filei18
  character(len=256) :: filei19, filei20
  character(len=256) :: name,unit         !netCDF attributes

!-----------------------------------------------------------------

! Determine input and output file names

  filei = '/ptmp/slevis/lsminput/old/world.ascii'
  filei1= '/ptmp/slevis/lsminput/old/clm/Clay2.map'
  filei2= '/ptmp/slevis/lsminput/old/clm/Clay5.map'
  filei3= '/ptmp/slevis/lsminput/old/clm/Clay9.map'
  filei4= '/ptmp/slevis/lsminput/old/clm/Clay17.map'
  filei5= '/ptmp/slevis/lsminput/old/clm/Clay29.map'
  filei6= '/ptmp/slevis/lsminput/old/clm/Clay49.map'
  filei7= '/ptmp/slevis/lsminput/old/clm/Clay83.map'
  filei8= '/ptmp/slevis/lsminput/old/clm/Clay138.map'
  filei9= '/ptmp/slevis/lsminput/old/clm/Clay230.map'
  filei10= '/ptmp/slevis/lsminput/old/clm/Clay343.map'
  filei11= '/ptmp/slevis/lsminput/old/clm/Sand2.map'
  filei12= '/ptmp/slevis/lsminput/old/clm/Sand5.map'
  filei13= '/ptmp/slevis/lsminput/old/clm/Sand9.map'
  filei14= '/ptmp/slevis/lsminput/old/clm/Sand17.map'
  filei15= '/ptmp/slevis/lsminput/old/clm/Sand29.map'
  filei16= '/ptmp/slevis/lsminput/old/clm/Sand49.map'
  filei17= '/ptmp/slevis/lsminput/old/clm/Sand83.map'
  filei18= '/ptmp/slevis/lsminput/old/clm/Sand138.map'
  filei19= '/ptmp/slevis/lsminput/old/clm/Sand230.map'
  filei20= '/ptmp/slevis/lsminput/old/clm/Sand343.map'
  fileo = '/ptmp/slevis/lsminput/new/clm/5minx5min/mksrf_soitex.nc'

! -----------------------------------------------------------------
! Determine grid for input data 
!
! IGBP input data
!
! Data are igbp soil mapunits. Each mapunit can be thought of as a
! unique soil profile. Map resolution is 5min x 5min, stored in
! latitude bands, from north to south. In a given latitude band,
! data begin at the dateline (180W) and proceed eastward.

! Also from igbp we read the %sand and %clay that correspond to each
! mapunit by LSM soil layer. The raw data set created here, assembles all
! this data in one netcdf file. The final product goes from south to north.
! -----------------------------------------------------------------

! Define soil thicknesses and soil depths (model dependent)

  zsoi(1) = 0.0175
  zsoi(2) = 0.0451
  zsoi(3) = 0.0906
  zsoi(4) = 0.1656
  zsoi(5) = 0.2892
  zsoi(6) = 0.4930
  zsoi(7) = 0.8290
  zsoi(8) = 1.3829
  zsoi(9) = 2.2962
  zsoi(10) = 3.4332

  dzsoi(1) = zsoi(1) - 0.
  dzsoi(2) = zsoi(2) - zsoi(1)
  dzsoi(3) = zsoi(3) - zsoi(2)
  dzsoi(4) = zsoi(4) - zsoi(3)
  dzsoi(5) = zsoi(5) - zsoi(4)
  dzsoi(6) = zsoi(6) - zsoi(5)
  dzsoi(7) = zsoi(7) - zsoi(6)
  dzsoi(8) = zsoi(8) - zsoi(7)
  dzsoi(9) = zsoi(9) - zsoi(8)
  dzsoi(10) = zsoi(10) - zsoi(9)

! Define North, East, South, West edges of grid

  edge(1) =   90.
  edge(2) =  180.
  edge(3) =  -90.
  edge(4) = -180.

! Make latitudes and longitudes at center of grid cell

  dx = (edge(2)-edge(4)) / nlon
  dy = (edge(1)-edge(3)) / nlat

  do j = 1, nlat
     do i = 1, nlon
        latixy(i,j) = (edge(1)-dy/2.) - (j-1)*dy
        longxy(i,j) = (edge(4)+dx/2.) + (i-1)*dx
       end do
  end do

  lat(:) = latixy(1,:)
  lon(:) = longxy(:,1)

! -----------------------------------------------------------------
! create netcdf file
! -----------------------------------------------------------------

  call wrap_create (fileo, nf_clobber, ncid)
  call wrap_put_att_text (ncid, nf_global, 'data_type', 'igbp_soil_texture_data')

! Define dimensions

  call wrap_def_dim (ncid, 'lon' , nlon, dimlon_id)
  call wrap_def_dim (ncid, 'lat' , nlat, dimlat_id)
  call wrap_def_dim (ncid, 'number_of_layers'   , nlay      , dimlay_id)
  call wrap_def_dim (ncid, 'number_of_mapunits' , nmapunits , dimmapunits_id)
  call wrap_def_dim (ncid, 'max_value_mapunit'  , mapunitmax, dimmapunitmax_id)

! Define input file independent variables 

  name = 'lon'
  unit = 'degrees east'
  dim1_id(1) = dimlon_id
  call wrap_def_var (ncid,'LON', nf_float, 1, dim1_id, lon_id)
  call wrap_put_att_text (ncid, lon_id, 'long_name', name)
  call wrap_put_att_text (ncid, lon_id, 'units'    , unit)

  name = 'lat'
  unit = 'degrees north'
  dim1_id(1) = dimlat_id
  call wrap_def_var (ncid,'LAT', nf_float, 1, dim1_id, lat_id)
  call wrap_put_att_text (ncid, lat_id, 'long_name', name)
  call wrap_put_att_text (ncid, lat_id, 'units'    , unit)

  name = 'longitude-2d'
  unit = 'degrees east'
  dim2_id(1) = dimlon_id
  dim2_id(2) = dimlat_id
  call wrap_def_var (ncid, 'LONGXY', nf_float, 2, dim2_id, longxy_id)
  call wrap_put_att_text (ncid, longxy_id, 'long_name', name)
  call wrap_put_att_text (ncid, longxy_id, 'units'    , unit)

  name = 'latitude-2d'
  unit = 'degrees north'
  call wrap_def_var (ncid, 'LATIXY', nf_float, 2, dim2_id, latixy_id)
  call wrap_put_att_text (ncid, latixy_id, 'long_name', name)
  call wrap_put_att_text (ncid, latixy_id, 'units'    , unit)

  name = 'land mask'
  unit = 'unitless'
  call wrap_def_var (ncid ,'LANDMASK' ,nf_float, 2, dim2_id, landmask_id)
  call wrap_put_att_text (ncid, landmask_id, 'long_name', name)
  call wrap_put_att_text (ncid, landmask_id, 'units'    , unit)

! to possibly replace the next two variables
! find out about dimensioned variables
! (eg, see how pressure levels are treated)

  name = 'soil layer thickness'
  unit = 'm'
  dim1_id(1) = dimlay_id
  call wrap_def_var (ncid,'DZSOI', nf_float, 1, dim1_id, dzsoi_id)
  call wrap_put_att_text (ncid, dzsoi_id, 'long_name', name)
  call wrap_put_att_text (ncid, dzsoi_id, 'units'    , unit)

  name = 'soil layer depth'
  unit = 'm'
  dim1_id(1) = dimlay_id
  call wrap_def_var (ncid,'ZSOI', nf_float, 1, dim1_id, zsoi_id)
  call wrap_put_att_text (ncid, zsoi_id, 'long_name', name)
  call wrap_put_att_text (ncid, zsoi_id, 'units'    , unit)

  name = 'northern edge of surface grid'
  unit = 'degrees north'
  call wrap_def_var (ncid, 'EDGEN', nf_float, 0, 0, edgen_id)
  call wrap_put_att_text (ncid, edgen_id, 'long_name', name)
  call wrap_put_att_text (ncid, edgen_id, 'units'    , unit)

  name = 'eastern edge of surface grid'
  unit = 'degrees east'
  call wrap_def_var (ncid, 'EDGEE', nf_float, 0, 0, edgee_id)
  call wrap_put_att_text (ncid, edgee_id, 'long_name', name)
  call wrap_put_att_text (ncid, edgee_id, 'units'    , unit)

  name = 'southern edge of surface grid'
  unit = 'degrees north'
  call wrap_def_var (ncid, 'EDGES', nf_float, 0, 0, edges_id)
  call wrap_put_att_text (ncid, edges_id, 'long_name', name)
  call wrap_put_att_text (ncid, edges_id, 'units'    , unit)

  name = 'western edge of surface grid'
  unit = 'degrees east'
  call wrap_def_var (ncid, 'EDGEW', nf_float, 0, 0, edgew_id)
  call wrap_put_att_text (ncid, edgew_id, 'long_name', name)
  call wrap_put_att_text (ncid, edgew_id, 'units'    , unit)

! Define soil type and soil texture variables

  name = 'igbp soil mapunit'
  unit = 'unitless'
  dim2_id(1) = dimlon_id
  dim2_id(2) = dimlat_id
  call wrap_def_var (ncid ,'MAPUNITS' ,nf_float, 2, dim2_id, mapunit_id)
  call wrap_put_att_text (ncid, mapunit_id, 'long_name', name)
  call wrap_put_att_text (ncid, mapunit_id, 'units'    , unit)

  name = 'percent sand'
  unit = 'unitless'
  dim2_id(1) = dimmapunitmax_id
  dim2_id(2) = dimlay_id
  call wrap_def_var (ncid ,'PCT_SAND' ,nf_float, 2, dim2_id, pct_sand_id)
  call wrap_put_att_text (ncid, pct_sand_id, 'long_name', name)
  call wrap_put_att_text (ncid, pct_sand_id, 'units'    , unit)