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