convert_vegtype.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

program make_surftype

  implicit none
#include <netcdf.inc>

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

  integer, parameter :: r8 = selected_real_kind(12)

  integer, parameter :: nlon = 720  !input grid : longitude points
  integer, parameter :: nlat = 360  !input grid : latitude  points
  integer, parameter :: nlsm = 28   !number of LSM surface types

  real(r8) :: longxy(nlon,nlat)     !longitude dimension array (2d)        
  real(r8) :: latixy(nlon,nlat)     !longitude dimension array (2d)
  real(r8) :: lon(nlon)            
  real(r8) :: lat(nlat)
  real(r8) :: edge(4)
  real(r8) :: dx,dy

  integer :: surtyp(nlon,nlat)      !input surface type

  integer :: dimlon_id                    !netCDF dimension id
  integer :: dimlat_id                    !netCDF dimension id
  integer :: nvegmax                !maximum value for input surface data

  integer :: lon_id                 !longitude array id
  integer :: lat_id                 !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 :: surftyp_id             !surface type id  

  integer :: i,j,k,m,n              !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                 !netCDF status

  integer :: miss=99999             !missing data indicator
  integer :: surftyp_i(nlon,nlat)   !input surface type (Olson or LSM)
  integer :: surftyp_o(nlon,nlat)   !output surface type (LSM)
  integer :: in2lsm(100)            !LSM surface type for each input type
  integer :: i_w                    !grid cell to west
  integer :: i_e                    !grid cell to east
  integer :: j_n                    !grid cell to north
  integer :: j_s                    !grid cell to south
  
  character(len=80) :: filei, fileo
  character(len=80) :: fvegtyp
  character(len=80) :: name,unit


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

! Determine input and output file names

  filei = 'olson.dat'
  fileo = 'olson.nc'

! -----------------------------------------------------------------
! Read in input data. 
! Data can be in one of two forms:
! OLSON vegetation types or LSM vegetation types. 
!  o The LSM vegetation types are from 0 to [nlsm]. 
!  o The OLSON vegetation types are from 0 to 73 and 
!    are mapped to LSM types.
!
! Data are 1/2 x 1/2 degree, stored in latitude bands, 
! from north to south. In a given latitude band, data begin 
! at Dateline (180W) and proceed eastward. So first data  
! point (x(1,1)) is a box centered at 89.75N, 179.75W. 
!
!   90.0N  ---------------------  
!          |         |         |
!          |    x    |    x    |
!          |  (1,1)  |  (2,1)  |
!          |         |         |
!   89.5N  --------------------- 
!        180.0W    179.5W    179.0W
! -----------------------------------------------------------------

! 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(3)+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','surface_type_data')

! Define dimensions

  call wrap_def_dim (ncid, 'LON' , nlon, dimlon_id)
  call wrap_def_dim (ncid, 'LAT' , nlat, dimlat_id)

! Define grid variables

  dim1_id(1) = dimlon_id

  dim2_id(1) = dimlon_id
  dim2_id(2) = dimlat_id

  name = 'lon'
  unit = 'degrees east'
  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'
  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 = '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 = '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 = '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 = '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)

  name = 'surface_type'
  unit = 'unitless'
  call wrap_def_var (ncid ,'SURFACE_TYPE' ,nf_int, 2, dim2_id, surftyp_id)
  call wrap_put_att_text (ncid, surftyp_id, 'long_name', name)
  call wrap_put_att_text (ncid, surftyp_id, 'units'    , unit)

  status = nf_enddef(ncid)

! -----------------------------------------------------------------
! read in formatted surface data
! -----------------------------------------------------------------

  open (unit=ndata,file=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
  nvegmax = 0
  do j = 1, nlat
     do i = 1, nlon
        read (ndata,*) surftyp_i(i,j)
        nvegmax = max(nvegmax,surftyp_i(i,j))
     end do
  end do
  close(ndata)

! -----------------------------------------------------------------
! Olson data : Convert input surface types to LSM surface types
! -----------------------------------------------------------------

  if (nvegmax > nlsm) then

     do j = 1, nlat
        do i = 1, nlon

           k = surftyp_i(i,j)
           
! There are several values (2, 6, 8) in the data that are not defined. 
! Use neighboring cells in this order: west, east, north, south

           i_w = max(   1,i-1)
           i_e = min(nlon,i+1)
           j_n = min(   1,j-1)
           j_s = max(nlat,j+1)
           
           if (k==2 .or. k==6 .or. k==8) k = surftyp_i(i_w,j  )
           if (k==2 .or. k==6 .or. k==8) k = surftyp_i(i_e,j  )
           if (k==2 .or. k==6 .or. k==8) k = surftyp_i(i  ,j_n)
           if (k==2 .or. k==6 .or. k==8) k = surftyp_i(i  ,j_s)
           
! Split antarctica (17) into polar desert (69) and ice (70)
 
           if (k == 17) then
              if (j <= 313) then
                 k = 69
              else
                 k = 70
              end if
           end if

! 61 (eastern south taiga) will be classified as needleleaf deciduous tree. 
! Change 61 to 20 (main taiga = needleleaf evergreen tree) based on longitude

           if (k==61 .and. i<=576) k = 20 

! 61 (eastern south taiga) will be classified needleleaf deciduous tree. 
! Create additional needleleaf deciduous tree from 21 (main taiga) and
! 60 (southern taiga) based on longitude

           if (k==21 .and. i>=555) k = 61  
           if (k==60 .and. i>=582) k = 61   

! Change 26 (warm mixed) to broad-leaved humid forest based on latitude