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