program convert_nc

  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 :: numpft = 16       !number of plant types
  integer, parameter :: numgr = 6         !number of pft with grazing info

  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) :: dummy

  real(r8) :: landmask(nlon,nlat)         !fraction of land
  real(r8) :: pct_tree(nlon,nlat)         !fraction tree
  real(r8) :: pct_grass(nlon,nlat)        !fraction tree
  real(r8) :: pct_graing(nlon,nlat)       !fraction tree
  real(r8) :: pct_graz(nlon,nlat)         !fraction grazed land
  real(r8) :: pct_pft(nlon,nlat,0:numpft) !percent pft
  real(r8) :: pct_pft_gr(nlon,nlat,numgr) !percent shrub/grass pft grazed

  integer :: dimlon_id                    !netCDF dimension id
  integer :: dimlat_id                    !netCDF dimension id
  integer :: dimpft_id                    !netCDF dimension id
  integer :: dimgr_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 :: landmask_id                  !landmask id
  integer :: pct_tree_id                  !pct_tree id
  integer :: pct_grass_id                 !pct_grass id
  integer :: pct_graz_id                  !pct_graz id
  integer :: pct_pft_id                   !pct_pft id
  integer :: pct_pft_gr_id                !pct_pft_gr id

  integer :: i,j,m,l,k,iii                !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 :: dim3a_id(3)                  !netCDF dimension id for 3-d variables
  integer :: dim3b_id(3)                  !netCDF dimension id for 3-d variables
  integer :: status                       !status
  integer :: start_year                   !loop start
  integer :: end_year                     !loop end
  integer :: year_int                     !interval between years processed
  integer :: scen                         !select scenario to process
  integer :: CN                           !indicator for CN processing
    
  character(len=80) :: filei, fileo       !input,output filenames
  character(len=80) :: name, unit         !netCDF attributes
  character(len=80) :: name_in
, name_out
  character(len=35) :: in_1_hist
  character(len=35) :: in_2_hist
  character(len=35) :: out_1_hist
  character(len=35) :: out_2_hist
  character(len=35) :: in_1_2
  character(len=35) :: in_2_2
  character(len=35) :: out_1_2
  character(len=35) :: out_2_2
  character(len=35) :: in_1_3
  character(len=35) :: in_2_3
  character(len=35) :: out_1_3
  character(len=35) :: out_2_3
  character(len=35) :: midsect
  character(len=6) :: timestamp
!
!  check to see if processing CN or regular files
!
  write(6,*) ' Do you want the CN version?'
  write(6,*) ' -- this converts temperate/boreal DBT to tropical DBT'
  write(6,*) ' 1 for regular conversion'
  write(6,*) ' 2 for CN conversion'
  read(5,*)CN
!
!  ask for scenario to process and set data accordingly
!
  write(6,*) 'enter the scenario to process'
  write(6,*) ' Historical   = 1'
  write(6,*) ' Future - A2  = 2'
  write(6,*) ' Future - B1  = 3'
  write(6,*) ' Future - B2  = 4'
  write(6,*) ' Future - A1B = 5'
  read(5,*) scen
!
  write(6,*) 'enter the timestamp for the folders mmddyy'
  read(5,'(a6)') timestamp
!
  if(scen .eq. 1)then
   in_1_hist = 'output/pft-'
   in_2_hist = '-hist.5x.5'
   if (CN .eq. 1) then
    out_1_hist = 'historical_'
    out_2_hist = '.nc'
   else
    out_1_hist = 'historical_CN_'
    out_2_hist = '_CN.nc'
   endif
  elseif(scen .eq. 2)then
   in_1_2 = 'output/pft-'
   in_2_2 = '-A2.5x.5'
   if (CN .eq. 1) then
    out_1_2 = 'A2_'
    out_2_2 = '_A2.nc'
   else
    out_1_2 = 'A2_CN_'
    out_2_2 = '_CN_A2.nc'
   endif
  elseif(scen .eq. 3)then
   in_1_2 = 'output/pft-'
   in_2_2 = '-B1.5x.5'
   if (CN .eq. 1) then
    out_1_2 = 'B1_'
    out_2_2 = '_B1.nc'
   else
    out_1_2 = 'B1_CN_'
    out_2_2 = '_CN_B1.nc'
   endif
  elseif(scen .eq. 4)then
   in_1_2 = 'output/pft-'
   in_2_2 = '-B2.5x.5'
   if (CN .eq. 1) then
    out_1_2 = 'B2_'
    out_2_2 = '_B2.nc'
   else
    out_1_2 = 'B2_CN_'
    out_2_2 = '_CN_B2.nc'
   endif
  elseif(scen .eq. 5)then
   in_1_3 = 'output/pft-'
   in_2_3 = '-A1B.5x.5'
   if (CN .eq. 1) then
    out_1_3 = 'A1B_'
    out_2_3 = '_A1B.nc'
   else
    out_1_3 = 'A1B_CN_'
    out_2_3 = '_CN_A1B.nc'
   endif
  endif
  midsect = '/mksrf_pft_'
!-----------------------------------------------------------------

! Determine input and output file names
!
!  enter the years over which to create the files
!
  write(6,*) 'enter the start year'
  read(5,*)start_year
  write(6,*) 'enter the end year'
  read(5,*)end_year
  write(6,*) 'enter interval between years'
  read(5,*)year_int
 
 do iii = start_year,end_year,year_int
  write(6,*)'On year ', iii
  if(scen .eq. 1)then
   write(name_in,'(a11,i4,a10)')in_1_hist,iii,in_2_hist
   if (CN .eq. 1) then
    write(name_out,'(a11,a6,a11,i4,a3)')out_1_hist,timestamp,midsect,iii,out_2_hist
   else
    write(name_out,'(a14,a6,a11,i4,a6)')out_1_hist,timestamp,midsect,iii,out_2_hist
   endif
  elseif(scen .gt. 1 .and. scen .lt. 4)then
   write(name_in,'(a11,i4,a8)')in_1_2,iii,in_2_2
   if (CN .eq. 1) then
    write(name_out,'(a3,a6,a11,i4,a6)')out_1_2,timestamp,midsect,iii,out_2_2
   else
    write(name_out,'(a6,a6,a11,i4,a9)')out_1_2,timestamp,midsect,iii,out_2_2
   endif
  else
   write(name_in,'(a11,i4,a9)')in_1_3,iii,in_2_3
   if (CN .eq. 1) then
    write(name_out,'(a4,a6,a11,i4,a7)')out_1_3,timestamp,midsect,iii,out_2_3
   else
    write(name_out,'(a7,a6,a11,i4,a10)')out_1_3,timestamp,midsect,iii,out_2_3
   endif
  endif
 
  filei = name_in
  fileo = name_out
!
! -----------------------------------------------------------------
! Determine grid for input data 
!
! Data are 0.5 x 0.5 degree, stored in latitude bands,
! from south to north. In a given latitude band, data begin
! at the dateline and proceed eastward. So first data
! point (x(1,1)) is a box centered at 89.75S, 179.75W
!
!   89.5S  ---------------------
!          |         |         |
!          |    x    |    x    |
!          |  (1,1)  |  (2,1)  |
!          |         |         |
!   90.0S  ---------------------
!        180.W     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
! -----------------------------------------------------------------

! Define global attributes

  call wrap_create (fileo, nf_clobber, ncid)
  call wrap_put_att_text (ncid, nf_global, 'data_type', 'pft_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, 'pft', numpft+1, dimpft_id)
  call wrap_def_dim (ncid, 'graz', numgr  , dimgr_id)

! Define grid 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 pft and grazing variables

  name = 'land mask'
  unit = 'unitless'
  dim2_id(1) = dimlon_id
  dim2_id(2) = dimlat_id
  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)

  name = 'percent pft'
  unit = 'unitless'
  dim3a_id(1) = dimlon_id
  dim3a_id(2) = dimlat_id
  dim3a_id(3) = dimpft_id
  call wrap_def_var (ncid ,'PCT_PFT' ,nf_float, 3, dim3a_id, pct_pft_id)
  call wrap_put_att_text (ncid, pct_pft_id, 'long_name', name)
  call wrap_put_att_text (ncid, pct_pft_id, 'units'    , unit)

  name = 'tree percent'
  unit = 'unitless'
  dim2_id(1) = dimlon_id
  dim2_id(2) = dimlat_id
  call wrap_def_var (ncid ,'PCT_TREE' ,nf_float, 2, dim2_id, pct_tree_id)
  call wrap_put_att_text (ncid, pct_tree_id, 'long_name', name)
  call wrap_put_att_text (ncid, pct_tree_id, 'units'    , unit)

  name = 'grass percent'
  unit = 'unitless'
  dim2_id(1) = dimlon_id
  dim2_id(2) = dimlat_id
  call wrap_def_var (ncid ,'PCT_GRASS' ,nf_float, 2, dim2_id, pct_grass_id)
  call wrap_put_att_text (ncid, pct_grass_id, 'long_name', name)
  call wrap_put_att_text (ncid, pct_grass_id, 'units'    , unit)

  name = 'gazing percent'
  unit = 'unitless'
  dim2_id(1) = dimlon_id
  dim2_id(2) = dimlat_id
  call wrap_def_var (ncid ,'TOT_GRAZING_PCT' ,nf_float, 2, dim2_id, pct_graz_id)
  call wrap_put_att_text (ncid, pct_graz_id, 'long_name', name)
  call wrap_put_att_text (ncid, pct_graz_id, 'units'    , unit)

  name = 'percent pft'
  unit = 'unitless'
  dim3b_id(1) = dimlon_id
  dim3b_id(2) = dimlat_id
  dim3b_id(3) = dimgr_id
  call wrap_def_var (ncid ,'PCT_PFT_GRAZED' ,nf_float, 3, dim3b_id, pct_pft_gr_id)
  call wrap_put_att_text (ncid, pct_pft_gr_id, 'long_name', name)
  call wrap_put_att_text (ncid, pct_pft_gr_id, 'units'    , unit)

! End of definitions

  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,*) dummy,dummy,landmask(i,j),(pct_pft(i,j,m),m=0,numpft), &
                       pct_tree(i,j),pct_grass(i,j),pct_graz(i,j),            &
                       (pct_pft_gr(i,j,l),l=1,numgr)
        if (landmask(i,j) == 100.) then
!
!   landmask change is part of original CCSM conversion
!
         landmask(i,j) = 1.
!
!  Assign the CN change
!
         
         if (CN .ne. 1) then
          if(lat(j) .gt. -23.5 .and. lat(j) .lt. 23.5) then
           pct_pft(i,j,6) = pct_pft(i,j,6) + pct_pft(i,j,7) + pct_pft(i,j,8)
           pct_pft(i,j,7) = 0._r8
           pct_pft(i,j,8) = 0._r8
          endif
         endif
!
        endif
        
     end do
  end do
  close(ndata)

! Write variables for original data

  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, pct_pft_id , pct_pft)
  call wrap_put_var_realx (ncid, pct_tree_id , pct_tree)
  call wrap_put_var_realx (ncid, pct_grass_id , pct_grass)
  call wrap_put_var_realx (ncid, pct_graz_id , pct_graz)
  call wrap_put_var_realx (ncid, pct_pft_gr_id , pct_pft_gr)

  call wrap_close(ncid)

 enddo
end program convert_nc

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

subroutine wrap_create (path, cmode, ncid)
  implicit none
  include 'netcdf.inc'
  integer, parameter :: r8 = selected_real_kind(12)
  character(len=*) path
  integer cmode, ncid, ret
  ret = nf_create (path, cmode, ncid)
  if (ret.ne.NF_NOERR) call handle_error (ret)
end subroutine wrap_create

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

subroutine wrap_def_dim (nfid, dimname, len, dimid)
  implicit none
  include 'netcdf.inc'
  integer, parameter :: r8 = selected_real_kind(12)
  integer :: nfid, len, dimid
  character(len=*) :: dimname
  integer ret
  ret = nf_def_dim (nfid, dimname, len, dimid)
  if (ret.ne.NF_NOERR) call handle_error (ret)
end subroutine wrap_def_dim

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

subroutine wrap_def_var (nfid, name, xtype, nvdims, vdims, varid)
  implicit none
  include 'netcdf.inc'
  integer, parameter :: r8 = selected_real_kind(12)
  integer :: nfid, xtype, nvdims, varid
  integer :: vdims(nvdims)
  character(len=*) :: name
  integer ret
  ret = nf_def_var (nfid, name, xtype, nvdims, vdims, varid)
  if (ret.ne.NF_NOERR) call handle_error (ret)
end subroutine wrap_def_var

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

subroutine wrap_put_att_text (nfid, varid, attname, atttext)
  implicit none
  include 'netcdf.inc'
  integer, parameter :: r8 = selected_real_kind(12)
  integer :: nfid, varid
  character(len=*) :: attname, atttext
  integer :: ret, siz
  siz = len_trim(atttext)
  ret = nf_put_att_text (nfid, varid, attname, siz, atttext)
  if (ret.ne.NF_NOERR) call handle_error (ret)
end subroutine wrap_put_att_text

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

subroutine wrap_put_var_realx (nfid, varid, arr)
  implicit none
  include 'netcdf.inc'
  integer, parameter :: r8 = selected_real_kind(12)
  integer :: nfid, varid
  real(r8) :: arr(*)
  integer :: ret
  ret = nf_put_var_double (nfid, varid, arr)
  if (ret.ne.NF_NOERR) call handle_error (ret)
end subroutine wrap_put_var_realx

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

subroutine wrap_put_var_int (nfid, varid, arr)
  implicit none
  include 'netcdf.inc'
  integer, parameter :: r8 = selected_real_kind(12)
  integer :: nfid, varid
  integer :: arr(*)
  integer :: ret
  ret = nf_put_var_int (nfid, varid, arr)
  if (ret.ne.NF_NOERR) call handle_error (ret)
end subroutine wrap_put_var_int
  
!===============================================================================

subroutine wrap_close (ncid)
  implicit none
  include 'netcdf.inc'
  integer, parameter :: r8 = selected_real_kind(12)
  integer :: ncid
  integer :: ret
  ret = nf_close (ncid)
  if (ret.ne.NF_NOERR) then
     write(6,*)'WRAP_CLOSE: nf_close failed for id ',ncid
     call handle_error (ret)
  end if
end subroutine wrap_close

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

subroutine handle_error(ret)
  implicit none
  include 'netcdf.inc'
  integer :: ret
  if (ret .ne. nf_noerr) then
     write(6,*) 'NCDERR: ERROR: ',nf_strerror(ret)
     call abort
  endif
end subroutine handle_error

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