convert_soicol.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

program convert_soicol

  implicit none
  include 'netcdf.inc'

!-----------------------------------------------------------------
! make soil colr netcdf file
!-----------------------------------------------------------------

  integer, parameter :: r8 = selected_real_kind(12)

! File specific settings

  integer, parameter :: nlon = 128        !input grid : longitude points
  integer, parameter :: nlat =  64        !input grid : latitude  points

  real(r8) :: gaulat(nlat)                !input grid: Gaussian latitudes
  real(r8) :: gauwt(nlat)                 !input grid: Gaussian weights

  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

  real(r8) :: soil_color(nlon,nlat)       !lsm soil color
  real(r8) :: landmask(nlon,nlat)         !land mask derived from soil color

  integer :: dimlon_id                    !netCDF dimension id
  integer :: dimlat_id                    !netCDF dimension id

  integer :: lon_id                       !1d longitude array id
  integer :: lat_id                       !1d latitude 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 :: soil_color_id                !soil color id
  integer :: landmask_id                  !landmask id

  integer :: i,j                          !indicis
  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) :: name,unit         !netCDF attributes

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

! Determine input and output file names

  filei = '/ptmp/slevis/lsmv2_2/input/0.5x0.5/bats.dat'
  fileo = '/ptmp/slevis/lsmv2_2/input/0.5x0.5/mksrf_soicol.nc'

! -----------------------------------------------------------------
! Determine grid for input data 
!
! BATS data are on T42 Gaussian grid, approximately 2.8 x 2.8 degrees,
! stored in latitude bands, from south to north. In a given latitude band, 
! data begin at Greenwich, centered on Greenwich, and proceed eastward. 
! -----------------------------------------------------------------

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

  edge(1) =  90.
  edge(2) = 358.59375
  edge(3) = -90.
  edge(4) =  -1.40625

! Make latitudes and longitudes at center of grid cell

  dx = (edge(2)-edge(4)) / nlon
  call mkgaulat (gaulat, gauwt, nlat)

  do j = 1, nlat
     do i = 1, nlon
        latixy(i,j) = -asin(gaulat(j)) * 180./(4.*atan(1.))
        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', 'soil_color_data')

! Define dimensions

  call wrap_def_dim (ncid, 'lon' , nlon, dimlon_id)
  call wrap_def_dim (ncid, 'lat' , nlat, dimlat_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'
  dim2_id(1) = dimlon_id
  dim2_id(2) = dimlat_id
  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 = '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 input file specific variables

  name = 'soil color'
  unit = 'unitless'
  dim2_id(1) = lon_id
  dim2_id(2) = lat_id
  call wrap_def_var (ncid ,'SOIL_COLOR' ,nf_float, 2, dim2_id, soil_color_id)
  call wrap_put_att_text (ncid, soil_color_id, 'long_name', name)
  call wrap_put_att_text (ncid, soil_color_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)

! End of definition

  status = nf_enddef(ncid)

! Read in formatted surface data

  open (unit=ndata,file=trim(filei),status='unknown',form='formatted',iostat=status)
  if (status .ne. 0) then
     write (6,*)'failed to open ',trim(filei),' on unit ',ndata,' ierr=',status
     stop
  end if
  do j = 1, nlat
     do i = 1, nlon
        read (ndata,*) soil_color(i,j)
        if (soil_color(i,j)<0 .or. soil_color(i,j)>8) then
           write (6,*) 'ERROR: BATS soil color = ',soil_color(i,j), &
                ' is not valid for lon,lat = ',i,j
           stop
        end if
        if (soil_color(i,j)==0) then
           landmask(i,j) = 0.
        else
           landmask(i,j) = 1.
        end if
     end do
  end do
  close(ndata)

! Create output file

  call wrap_put_var_realx (ncid, lon_id       , lon)
  call wrap_put_var_realx (ncid, lat_id       , lat)
  call wrap_put_var_realx (ncid, longxy_id    , longxy)
  call wrap_put_var_realx (ncid, latixy_id    , latixy)
  call wrap_put_var_realx (ncid, edgen_id     , edge(1))
  call wrap_put_var_realx (ncid, edgee_id     , edge(2))
  call wrap_put_var_realx (ncid, edges_id     , edge(3))
  call wrap_put_var_realx (ncid, edgew_id     , edge(4))
  call wrap_put_var_realx (ncid, landmask_id  , landmask)
  call wrap_put_var_realx (ncid, soil_color_id, soil_color)

! Close output file

  call wrap_close(ncid)

end program convert_soicol

!===============================================================================

      subroutine mkgaulat (a, w, k)

! ------------------------ code history ------------------------------
! source file       : mkgaulat.F
! purpose           : Gaussian latitudes
! date first created: July 1995 - lsm version 1
! by whom           : Gordon Bonan
! date last revised : Jan 2000 - lsm version 2
! by whom           : Mariana Vertenstein
! --------------------------------------------------------------------

      implicit none

      integer, parameter :: r8 = selected_real_kind(12)

! ------------------------ arguments ---------------------------------
      integer , intent(in) :: k       !number of latitudes pole to pole
      real(r8), intent(out):: a(k)    !sine of latitudes