  program clamp_rewrite

! Program to use CMOR routines to rewrite model results from the simulations
! for the Carbon-Land Model Intercomparison Project (C-LAMP).
!
! Forrest M. Hoffman <forrest@climatemodeling.org>
! Created: Fri Jun  8 11:50:01 EDT 2007


    use cmor_users_functions
    use netcdf
    use kind_mod
    use clm_mod, only: namelen, varnamelen, coord, var_data, &
      clm_read_coord, clm_copy_grid_data, clm_read_data, clm_convert_data, &
      clm_free_data
    implicit none

    integer, parameter :: maxf = 2500
    integer, parameter :: maxv = 200
    integer :: i, numf, numv

    integer, allocatable :: iaxis(:)          ! CMOR handle for axes (time, [hour,] [z-level,] latitude, longitude)
    integer              :: ivar              ! CMOR handle for variable
    integer              :: ierr              ! error flag
    integer              :: axis_num          ! axis counter
    integer              :: realization       ! model run realization
    character (len=128)  :: input_path        ! input path
    character (len=128)  :: fnames(maxf) = '' ! input file names
    integer              :: tshift(maxf) = 0  ! time shift for files
    character (len=128)  :: output_path       ! output path
    character (len=128)  :: input_table       ! CMOR table path
    character (len=32)   :: vnames(maxv) = '' ! input variable names
    character (len=256)  :: source = ''       ! source global attribute
    character (len=128)  :: experiment = ''   ! experiment_id global attribute
    logical              :: casa_flux_bug     ! scale fluxes due to bug in CASA'

    real(r4) :: imissing = 1.0e+36 ! missing-data flag on input
    real(r4) :: omissing = 1.0e+36 ! missing-data flag on output

    casa_flux_bug = .false.

    namelist /inparm/ input_table, output_path, experiment, source, &
      realization, input_path, fnames, tshift, vnames, casa_flux_bug

    open(10, file='namelist', status='old')
    ierr = 1
    do while(ierr /= 0)
      read(10, inparm, iostat=ierr)
      if (ierr < 0) then
        stop 'End of file on namelist read'
      end if
    end do
    close(10)

    ! Count number of input data files
    numf = 0
    do i = 1, maxf
       if (trim(fnames(i)) /= '') numf = numf + 1
    end do
    ! Count number of variables
    numv = 0
    do i = 1, maxv
       if (trim(vnames(i)) /= '') numv = numv + 1
    end do

    print *, 'Reading input files for coordinates and counts...'
    call clm_read_coord(numf, input_path, fnames, .true.)

    do i = 1, numv
      select case (vnames(i))
      ! gca
      case ('area')
         call clm_copy_grid_data(vnames(i))
      ! lbm
      case ('landmask')
         call clm_copy_grid_data(vnames(i))
      ! orog
      case ('topo')
         call clm_copy_grid_data(vnames(i))
      ! sftlf
      case ('landfrac')
         call clm_copy_grid_data(vnames(i))
      case default
         print *, 'Reading input files for variable ', vnames(i)
         call clm_read_data(numf, input_path, fnames, tshift, vnames(i))

         !print *,'var_data%varname = ',var_data%varname
         !print *,'var_data%ndims = ',var_data%ndims
         !print *,'var_data%time = ',var_data%time
         !print *,'var_data%long_name = ',trim(var_data%long_name)
         !print *,'var_data%missing_value = ',var_data%missing_value
         !print *,'var_data%units = ',trim(var_data%units)
      end select

      print *, 'Performing unit conversion on data...'
      call clm_convert_data(casa_flux_bug)

      print *, 'Initializing CMOR...'
      ierr = cmor_setup(inpath='./',netcdf_file_action='preserve', &
        set_verbosity=2, exit_control=2)

      print *, 'Identifying output data sets for CMOR...'
      ierr = cmor_dataset(                                              &
        outpath       = output_path,                                    &
        experiment_id = experiment,                                     &
        institution   = 'ORNL (Oak Ridge National Laboratory, Oak Ridge, Tennessee, USA)', &
        source        = trim(source),                                   &
        realization   = realization,                                    &
        calendar      = 'noleap',                                       &
        history       = 'Extracted from case '//trim(coord%case_id),    &
        comment       = 'Initial dataset: '//trim(coord%inidat)//       &
          '; surface dataset: '//trim(coord%surdat)//'; pft dataset: '//&
          trim(coord%pftdat)//'; rtm dataset: '//trim(coord%rtmdat),    &
        references    = 'http://www.climatemodeling.org/c-lamp',        &
        contact       = 'Forrest M. Hoffman <forrest@climatemodeling.org>')

      allocate(iaxis(var_data%ndims), stat=ierr)
      if (ierr /= 0) then
         print *, 'Cannot allocate iaxis'
         stop
      end if

      print *, 'Defining coordinates for CMOR output data...'
      axis_num = 1

      iaxis(axis_num) = cmor_axis(                                      &
        table       = input_table,                                      &
        table_entry = 'longitude',                                      &
        units       = 'degrees_east',                                   &
        length      = coord%xsize,                                      &
        coord_vals  = coord%x)
      axis_num = axis_num + 1

      iaxis(axis_num) = cmor_axis(                                      &
        table       = input_table,                                      &
        table_entry = 'latitude',                                       &
        units       = 'degrees_north',                                  &
        length      = coord%ysize,                                      &
        coord_vals  = coord%y)
      axis_num = axis_num + 1

      if (var_data%soil_layer_flag) then
         iaxis(axis_num) = cmor_axis(                                   &
           table       = input_table,                                   &
           table_entry = 'depth_soil',                                  &
           units       = 'm',                                           &
           length      = coord%zsoi_size,                               &
           coord_vals  = coord%zsoi)
         axis_num = axis_num + 1
      end if

      if (var_data%hour_flag) then
         iaxis(axis_num) = cmor_axis(                                   &
           table       = input_table,                                   &
           table_entry = 'hour',                                        &
           units       = 'hours',                                       &
           length      = coord%hr_size,                                 &
           coord_vals  = coord%hr)
         axis_num = axis_num + 1
      end if

      if (var_data%ndims > 2) then
        iaxis(axis_num) = cmor_axis(                                    &
          table       = input_table,                                    &
          table_entry = 'time',                                         &
          units       = 'days since 1798-01-01 00:00:00')
        axis_num = axis_num + 1
      end if

      print *, 'Defining CMOR output data variables...'
      if (var_data%positive == '') then
         ivar = cmor_variable(                                          &
           table         = input_table,                                 &
           table_entry   = var_data%out_varname,                        &
           original_name = var_data%varname,                            &
           units         = var_data%units,                              &
           missing_value = var_data%missing_value,                      &
           axis_ids      = iaxis)
      else
         ivar = cmor_variable(                                          &
           table         = input_table,                                 &
           table_entry   = var_data%out_varname,                        &
           original_name = var_data%varname,                            &
           units         = var_data%units,                              &
           positive      = var_data%positive,                           &
           missing_value = var_data%missing_value,                      &
           axis_ids      = iaxis)
      end if

      print *, 'Writing CMOR output...'
      select case (var_data%ndims)
         case (2)
            if (var_data%int_type) then
               ierr = cmor_write(                                       &
                 var_id    = ivar,                                      &
                 data      = var_data%int2d)
            else
               ierr = cmor_write(                                       &
                 var_id    = ivar,                                      &
                 data      = var_data%var2d)
            end if
         case (3)
            ierr = cmor_write(                                          &
              var_id    = ivar,                                         &
              data      = var_data%var3d,                               &
              time_vals = var_data%time,                                &
              time_bnds = var_data%time_bounds)
         case (4)
            ierr = cmor_write(                                          &
              var_id    = ivar,                                         &
              data      = var_data%var4d,                               &
              time_vals = var_data%time,                                &
              time_bnds = var_data%time_bounds)
         case default
            print *, 'Unable to handle data with ', var_data%ndims, 'dimensions'
            stop
      end select

      print *, 'Closing CMOR file(s)...'
      ierr = cmor_close()

      deallocate(iaxis)

      print *, 'Freeing data...'
      call clm_free_data()

    end do ! loop over variables

  end program clamp_rewrite
