!#################################################################
!
! user input 
!
!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!#################################################################

module user_input

  use GO, only : gol, goPr, goErr
  
  implicit none
  
  
  ! --- in/out --------------------------
  
  private
  
  public  ::  user_input_start
  

  ! --- const --------------------------------
    
  character(len=*), parameter  ::  mname = 'user_input'

  
contains


  ! =============================================================
  ! ===
  ! === start file
  ! ===
  ! =============================================================
  
  subroutine user_input_start( status )
  
    use GO         , only : gol, goPr, goErr
    use dims       , only : nregions
    use global_data, only : mass_dat
    use meteodata, only : m_dat

    ! --- in/out ------------------------------
    
    integer, intent(out)      ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/user_input_start'
    
    ! --- local -------------------------------
    
    integer       ::  n
    
    ! --- begin --------------------------------
    
    ! info ...

    call start_co2(status)
    IF_NOTOK_RETURN(status=1)    
    call start_c13(status)
    IF_NOTOK_RETURN(status=1)
    
    ! ok
    status=0
  
  end subroutine user_input_start

  subroutine start_co2(status)

  use GO,                  only : TrcFile, Init, Done, ReadRc
  use file_hdf   , only : THdfFile, TSds, Init, Done, ReadAttribute, ReadData, GetInfo
  use toolbox,     only : coarsen_emission
  use global_data, only : mass_dat
  use global_data,         only : rcfile
  use meteodata, only : m_dat
  use global_types,only : d3_data
  use ParTools,    only : ntracetloc, ntracet_ar, myid
  use Dims,        only : im,jm,lm,nregions, adv_scheme, nlon360,nlat180,at,bt
  use emission_common,only : EmisInputDir
  use chem_param, only : xmc, xmco2, xmair, names
  use GO,         only : goLoCase
  
  use meteodata       , only : lli, levi, sp_dat, m_dat
  use grid         , only : TllGridInfo, TLevelInfo,Init,Done,Fill3D
    
  implicit none

  !______________const______________________________
  character(len=*), parameter  ::  rname = mname//'/start_co2'
  !______________I/O______________________________
  integer, intent(out)      ::  status
  !______________local______________________________

  type(TrcFile)           ::  rcF
  real,dimension(:,:,:),pointer   ::m
  real,dimension(:,:,:,:),pointer   ::rm,rxm,rym,rzm
  real, dimension(:,:,:), pointer     ::  sp
  

  real,dimension(:,:,:,:),allocatable                  :: co2_file
  type(d3_data),dimension(nregions),target             :: co2_init

  integer,parameter                                    :: avg_field=1
  integer                                              :: dimheight, dimlat, dimlon, data_rank
  integer                                              :: region,l,imr,jmr,n, offsetn
  integer,dimension(4)                                 :: data_dims

  real                                                 :: co2_multiplier
  type(THdfFile)                   ::  hdf
  type(TSds)                       ::  sds_co2
  character(len=200) :: filename, varname


  type(TllGridInfo)                ::  lli_in
  type(TLevelInfo)                 ::  levi_in

  integer                          ::  xbeg_in, ybeg_in
  real                             ::  dx_in, dy_in
  integer                          ::  im_in, jm_in, lm_in
  real, allocatable                ::  at_in(:), bt_in(:)






    ! 34 (from standard) layer defintition:
    real, parameter :: at_tm5_std_34(0:34) = (/ &
    0.,              7.367743,        210.393890,      855.361755,     &
    2063.779785,     3850.913330,     6144.314941,     8802.356445,    &
    11632.758789,    14411.124023,    16899.468750,    17961.357422,   &
    18864.750000,    19584.330078,    20097.402344,    20384.480469,   &
    20429.863281,    20222.205078,    19755.109375,    19027.695313,   &
    18045.183594,    16819.474609,    15379.805664,    13775.325195,   &
    12077.446289,    10376.126953,    8765.053711,     7306.631348,    &
    6018.019531,     3960.291504,     1680.640259,     713.218079,     &
    298.495789,      95.636963,       0./)

    real, parameter  :: bt_tm5_std_34(0:34) = (/     &
    1.00000000,      0.99401945,      0.97966272,      0.95182151,     &
    0.90788388,      0.84737492,      0.77159661,      0.68326861,     &
    0.58616841,      0.48477158,      0.38389215,      0.33515489,     &
    0.28832296,      0.24393314,      0.20247590,      0.16438432,     &
    0.13002251,      0.09967469,      0.07353383,      0.05169041,     &
    0.03412116,      0.02067788,      0.01114291,      0.00508112,     &
    0.00181516,      0.00046139,      0.00007582,      0.00000000,     &
    0.00000000,      0.00000000,      0.00000000,      0.00000000,     &
    0.00000000,      0.00000000,      0.00000000/)



  write (gol,*) 'Initializing global mixing ratios from file' ; call goPr

  if(ntracetloc == 0) return
  do region=1,nregions
    imr=im(region)
    jmr=jm(region)
    allocate(co2_init(region)%d3(imr,jmr,lm(1)))
  enddo

  ! Get filename and variable name from rc-file

  filename = trim(EmisInputDir)//'/co2_init_2007.hdf'
  varname = 'co2_init_2007'

  call Init( rcF, rcfile, status )
  IF_NOTOK_RETURN(status=1)
  ! put out settings ?
  call ReadRc( rcF, 'init.co2.filename', filename, status, default=filename )
  IF_ERROR_RETURN(status=1)
  call ReadRc( rcF, 'init.co2.varname', varname, status, default=varname )
  IF_ERROR_RETURN(status=1)
  call ReadRc( rcF, 'init.co2.multiply', co2_multiplier, status, default=1.0 )
  IF_ERROR_RETURN(status=1)
  ! close rcfile:
  call Done( rcF, status  )
  IF_NOTOK_RETURN(status=1)


  write (gol,*) '.... Opening file with initial CO2 mixing ratios for reading: '//trim(filename) ; call goPr
 
  call Init( hdf, trim(filename), 'read', status )
  IF_NOTOK_RETURN(status=1)

  call Init( sds_co2, hdf, trim(varname), status )  
  IF_NOTOK_RETURN(status=1)


  call GetInfo( sds_co2, status,data_rank=data_rank, data_dims=data_dims )
  IF_NOTOK_RETURN(status=1)

  if(data_rank /= 3.and.data_rank /= 4) then
    write (gol,*) 'Rank of requested dataset ('//trim(varname)//') is not 3 or 4 (time,lon,lat,levels) ...error' ; call goPr
    IF_NOTOK_RETURN(status=1)
  endif

  allocate(co2_file(data_dims(1),data_dims(2),data_dims(3),1) )

  if(data_rank.eq.3) then 
     write (gol,*) '.... Retrieving variable '//trim(varname) ; call goPr
  endif

  if(data_rank.eq.4) then
     write (gol,*) '.... Retrieving first time step of variable '//trim(varname) ; call goPr
  endif

  call ReadData( sds_co2, co2_file, status)
  IF_NOTOK_RETURN(status=1)

  call Done( sds_co2, status )  
  IF_NOTOK_RETURN(status=1)

  write (gol,*) '.... Success' ; call goPr
  call Done( hdf, status )  
  IF_NOTOK_RETURN(status=1)
  write (gol,*) '.... File closed ' ; call goPr

  !WP! Multiply with specified factor

  co2_file = co2_file * co2_multiplier

  !WP Now grid the input file to the dimensions of the current model run

  !IV coarsen_emission does not regrid in the vertical, therefore must be replaced by Fill3D
  !call coarsen_emission('co2_init',data_dims(1),data_dims(2),data_dims(3),co2_file(:,:,:,1),co2_init,avg_field) ! not needed

  im_in=data_dims(1)
  jm_in=data_dims(2)
  lm_in=data_dims(3)

  !WP! We are assuming a 1x1 degree input file for now, this can be changed in the future

  if(data_dims(1) /= 360) then
    write (gol,'("input from initial concentration file implemented for 1 degree longitude only ...")'); call goErr
    TRACEBACK; status=1; return
  endif
  if(data_dims(2) /= 180) then
    write (gol,'("input from initial concentration file implemented for 1 degree latitude only ...")'); call goErr
    TRACEBACK; status=1; return
  endif

  !WP! We are assuming a 1x1 degree input file for now, this can be changed in the future

  xbeg_in = -180.0
  dx_in = 360.0/im_in
  ybeg_in = -90.0
  dy_in = 180.0/jm_in

  call Init( lli_in, xbeg_in+0.5*dx_in, dx_in, im_in, &
                   ybeg_in+0.5*dy_in, dy_in, jm_in, status )

  IF_NOTOK_RETURN(status=1)


  allocate( at_in(0:lm_in) )
  allocate( bt_in(0:lm_in) )

  select case ( lm_in )
      case ( lm(1) )  !WP! If vertical dimension equals current setup, use at, bt from current run
        at_in = at
        bt_in = bt
        levi_in = levi

      case ( 34 )
        at_in = at_tm5_std_34
        bt_in = bt_tm5_std_34
        call Init( levi_in, lm_in, at_in, bt_in, status )
        IF_NOTOK_RETURN(status=1)

       case default
          write (gol,'("input from initial concentration file implemented for 25 or 34 levels only ...")'); call goErr
          TRACEBACK; status=1; return

  end select
  deallocate( at_in )
  deallocate( bt_in )

  write (gol,*) '.... Global mean CO2 mixing ratio in input field     : ',sum(co2_file)/(im_in*jm_in*lm_in) ; call goPr

  do region=1,nregions
      m => m_dat(region)%data
      rm => mass_dat(region)%rm_t
      rxm => mass_dat(region)%rxm_t
      rym => mass_dat(region)%rym_t
      rzm => mass_dat(region)%rzm_t
      sp => sp_dat(region)%data

      call Fill3D( lli(region), levi, 'n', sp(1:im(region),1:jm(region),1), co2_init(region)%d3, &
                         lli_in, levi_in, co2_file(:,:,:,1), 'mass-aver', status )

      write (gol,*) '.... Mixing ratios coarsened to TM5 grid ' ; call goPr

      write (gol,*) '.... Mean CO2 mixing ratio in TM5 field region nr:',region,sum(co2_init(region)%d3)/(im(region)*jm(region)*lm(region)) ; call goPr


      
      do n=1,ntracetloc

         ! There is a bug with TM5 using the reduced grid.  Apparently
         ! negative concentrations are not permitted in this scenario.
         ! Thus for the forward simulation, we start each component
         ! with the background CO2 field, on which sinks can act.  To
         ! get the component concentration, the bg field must be
         ! subtracted.

         !WP! Get the global tracer number, we only initialize ico2_bg (if we
         !WP! run fwd) or we initialize all co2 tracers (in case of inv)

         offsetn = sum(ntracet_ar(0:myid-1))

         if (goLoCase(trim(names(offsetn+n)))=='co2'.or.goLoCase(trim(names(offsetn+n)))=='co2_bg') then 

            if (region==1) then
                write (gol,'(a,i4,a,a)') '.... Found tracer ',offsetn+n,' with name: ',trim(names(offsetn+n)) ; call goPr
                write (gol,'(a)')      '.... Initializing with CO2 from file'               ; call goPr
            endif
        
            rm(1:im(region),1:jm(region),1:lm(region),n) = &
              co2_init(region)%d3*1.e-6*m(1:im(region),1:jm(region),1:lm(region))*xmco2/xmair
        ! else

        !    rm(:,:,:,:) = 0.0

         endif
           

      enddo

#ifdef slopes
      rxm(:,:,:,:) = 0.0
      rym(:,:,:,:) = 0.0
      rzm(:,:,:,:) = 0.0
#endif
      nullify(m)
      nullify(rm)
      nullify(rxm)
      nullify(rym)
      nullify(rzm)
      nullify(sp)      
  enddo

  deallocate(co2_file)


  do region=1,nregions
    deallocate(co2_init(region)%d3)
  enddo

  write (gol,'(a,a)')  rname ,' Done'               ; call goPr
  
  end subroutine start_co2















 subroutine start_c13(status)

  use GO,                  only : TrcFile, Init, Done, ReadRc
  use file_hdf   , only : THdfFile, TSds, Init, Done, ReadAttribute, ReadData, GetInfo
  use toolbox,     only : coarsen_emission
  use global_data, only : mass_dat
  use global_data,         only : rcfile
  use meteodata, only : m_dat
  use global_types,only : d3_data
  use ParTools,    only : ntracetloc, ntracet_ar, myid
  use Dims,        only : im,jm,lm,nregions, adv_scheme, nlon360,nlat180,at,bt
  use emission_common,only : C13InputDir,EmisInputDir
  use chem_param, only : xmc13, xm13co2,xmco2, xmair, names
  use GO,         only : goLoCase
  
  use meteodata       , only : lli, levi, sp_dat, m_dat
  use grid         , only : TllGridInfo, TLevelInfo,Init,Done,Fill3D
    
  implicit none

  !______________const______________________________
  character(len=*), parameter  ::  rname = mname//'/start_c13'
  !______________I/O______________________________
  integer, intent(out)      ::  status
  !______________local______________________________

  type(TrcFile)           ::  rcF
  real,dimension(:,:,:),pointer   ::m
  real,dimension(:,:,:,:),pointer   ::rm,rxm,rym,rzm
  real, dimension(:,:,:), pointer     ::  sp
  

  real,dimension(:,:,:,:),allocatable                  :: c13_file,co2_file
  type(d3_data),dimension(nregions),target             :: c13_init,co2_init

  integer,parameter                                    :: avg_field=1
  integer                                              :: dimheight, dimlat, dimlon, data_rank
  integer                                              :: region,l,imr,jmr,n, offsetn
  integer,dimension(4)                                 :: data_dims

  real                                                 :: c13_multiplier,co2_multiplier
  type(THdfFile)                   ::  hdf
  type(TSds)                       ::  sds_c13,sds_co2
  character(len=200) :: filename, varname




  type(TllGridInfo)                ::  lli_in
  type(TLevelInfo)                 ::  levi_in

  integer                          ::  xbeg_in, ybeg_in
  real                             ::  dx_in, dy_in
  integer                          ::  im_in, jm_in, lm_in
  real, allocatable                ::  at_in(:), bt_in(:)


    ! 34 (from standard) layer defintition:
    real, parameter :: at_tm5_std_34(0:34) = (/ &
    0.,              7.367743,        210.393890,      855.361755,     &
    2063.779785,     3850.913330,     6144.314941,     8802.356445,    &
    11632.758789,    14411.124023,    16899.468750,    17961.357422,   &
    18864.750000,    19584.330078,    20097.402344,    20384.480469,   &
    20429.863281,    20222.205078,    19755.109375,    19027.695313,   &
    18045.183594,    16819.474609,    15379.805664,    13775.325195,   &
    12077.446289,    10376.126953,    8765.053711,     7306.631348,    &
    6018.019531,     3960.291504,     1680.640259,     713.218079,     &
    298.495789,      95.636963,       0./)

    real, parameter  :: bt_tm5_std_34(0:34) = (/     &
    1.00000000,      0.99401945,      0.97966272,      0.95182151,     &
    0.90788388,      0.84737492,      0.77159661,      0.68326861,     &
    0.58616841,      0.48477158,      0.38389215,      0.33515489,     &
    0.28832296,      0.24393314,      0.20247590,      0.16438432,     &
    0.13002251,      0.09967469,      0.07353383,      0.05169041,     &
    0.03412116,      0.02067788,      0.01114291,      0.00508112,     &
    0.00181516,      0.00046139,      0.00007582,      0.00000000,     &
    0.00000000,      0.00000000,      0.00000000,      0.00000000,     &
    0.00000000,      0.00000000,      0.00000000/)

  
  write (gol,*) 'Initializing global mixing ratios from file' ; call goPr  

  if(ntracetloc == 0) return
  do region=1,nregions
    imr=im(region)
    jmr=jm(region)
    allocate(co2_init(region)%d3(imr,jmr,lm(1)))
  enddo

  ! Get filename and variable name from rc-file

  filename = trim(EmisInputDir)//'/co2_init_2007.hdf'
  varname = 'co2_init_2007'

  call Init( rcF, rcfile, status )
  IF_NOTOK_RETURN(status=1)
  ! put out settings ?
  call ReadRc( rcF, 'init.co2.filename', filename, status, default=filename )
  IF_ERROR_RETURN(status=1)
  call ReadRc( rcF, 'init.co2.varname', varname, status, default=varname )
  IF_ERROR_RETURN(status=1)
  call ReadRc( rcF, 'init.co2.multiply', co2_multiplier, status, default=1.0 )
  IF_ERROR_RETURN(status=1)
  ! close rcfile:
  call Done( rcF, status  )
  IF_NOTOK_RETURN(status=1)


  write (gol,*) '.... Opening file with initial CO2 mixing ratios for reading: '//trim(filename) ; call goPr
 
  call Init( hdf, trim(filename), 'read', status )
  IF_NOTOK_RETURN(status=1)

  call Init( sds_co2, hdf, trim(varname), status )  
  IF_NOTOK_RETURN(status=1)


  call GetInfo( sds_co2, status,data_rank=data_rank, data_dims=data_dims )
  IF_NOTOK_RETURN(status=1)

  if(data_rank /= 3.and.data_rank /= 4) then
    write (gol,*) 'Rank of requested dataset ('//trim(varname)//') is not 3 or 4 (time,lon,lat,levels) ...error' ; call goPr
    IF_NOTOK_RETURN(status=1)
  endif

  allocate(co2_file(data_dims(1),data_dims(2),data_dims(3),1) )

  if(data_rank.eq.3) then 
     write (gol,*) '.... Retrieving variable '//trim(varname) ; call goPr
  endif

  if(data_rank.eq.4) then
     write (gol,*) '.... Retrieving first time step of variable '//trim(varname) ; call goPr
  endif

  call ReadData( sds_co2, co2_file, status)
  IF_NOTOK_RETURN(status=1)

  call Done( sds_co2, status )  
  IF_NOTOK_RETURN(status=1)

  write (gol,*) '.... Success' ; call goPr
  call Done( hdf, status )  
  IF_NOTOK_RETURN(status=1)
  write (gol,*) '.... File closed ' ; call goPr

  !WP! Multiply with specified factor

  co2_file = co2_file * co2_multiplier

  !WP Now grid the input file to the dimensions of the current model run

  !IV coarsen_emission does not regrid in the vertical, therefore must be replaced by Fill3D
  !call coarsen_emission('co2_init',data_dims(1),data_dims(2),data_dims(3),co2_file(:,:,:,1),co2_init,avg_field) ! not needed

  im_in=data_dims(1)
  jm_in=data_dims(2)
  lm_in=data_dims(3)

  !WP! We are assuming a 1x1 degree input file for now, this can be changed in the future

  if(data_dims(1) /= 360) then
    write (gol,'("input from initial concentration file implemented for 1 degree longitude only ...")'); call goErr
    TRACEBACK; status=1; return
  endif
  if(data_dims(2) /= 180) then
    write (gol,'("input from initial concentration file implemented for 1 degree latitude only ...")'); call goErr
    TRACEBACK; status=1; return
  endif

  !WP! We are assuming a 1x1 degree input file for now, this can be changed in the future

  xbeg_in = -180.0
  dx_in = 360.0/im_in
  ybeg_in = -90.0
  dy_in = 180.0/jm_in

  call Init( lli_in, xbeg_in+0.5*dx_in, dx_in, im_in, &
                   ybeg_in+0.5*dy_in, dy_in, jm_in, status )

  IF_NOTOK_RETURN(status=1)


  allocate( at_in(0:lm_in) )
  allocate( bt_in(0:lm_in) )

  select case ( lm_in )
      case ( lm(1) )  !WP! If vertical dimension equals current setup, use at, bt from current run
        at_in = at
        bt_in = bt
        levi_in = levi

      case ( 34 )
        at_in = at_tm5_std_34
        bt_in = bt_tm5_std_34
        call Init( levi_in, lm_in, at_in, bt_in, status )
        IF_NOTOK_RETURN(status=1)

       case default
          write (gol,'("input from initial concentration file implemented for 25 or 34 levels only ...")'); call goErr
          TRACEBACK; status=1; return

  end select
  deallocate( at_in )
  deallocate( bt_in )

 ! print*,sum(co2_file),im_in,jm_in,lm_in
  write (gol,*) '.... Global mean CO2 mixing ratio in input field     : ',sum(co2_file)/(im_in*jm_in*lm_in) ; call goPr

  do region=1,nregions
      m => m_dat(region)%data
      rm => mass_dat(region)%rm_t
      rxm => mass_dat(region)%rxm_t
      rym => mass_dat(region)%rym_t
      rzm => mass_dat(region)%rzm_t
      sp => sp_dat(region)%data

      call Fill3D( lli(region), levi, 'n', sp(1:im(region),1:jm(region),1), co2_init(region)%d3, &
                         lli_in, levi_in, co2_file(:,:,:,1), 'mass-aver', status )

      write (gol,*) '.... Mixing ratios coarsened to TM5 grid ' ; call goPr

      write (gol,*) '.... Mean CO2 mixing ratio in TM5 field region nr:',region,sum(co2_init(region)%d3)/(im(region)*jm(region)*lm(region)) ; call goPr
      
      
#ifdef slopes
      rxm(:,:,:,:) = 0.0
      rym(:,:,:,:) = 0.0
      rzm(:,:,:,:) = 0.0
#endif
      nullify(m)
      nullify(rm)
      nullify(rxm)
      nullify(rym)
      nullify(rzm)
      nullify(sp)   
  enddo    

  write (gol,*) 'Initializing global 13C mixing ratios from file' ; call goPr

  if(ntracetloc == 0) return
  do region=1,nregions
    imr=im(region)
    jmr=jm(region)
    allocate(c13_init(region)%d3(imr,jmr,lm(1)))
  enddo

  ! Get filename and variable name from rc-file

  filename = trim(C13InputDir)//'/d13c_init_2007_offset2.hdf'
  varname = 'd13c_init_2007'

  call Init( rcF, rcfile, status )
  IF_NOTOK_RETURN(status=1)
  ! put out settings ?
  call ReadRc( rcF, 'init.c13.filename', filename, status, default=filename )
  IF_ERROR_RETURN(status=1)
  call ReadRc( rcF, 'init.c13.varname', varname, status, default=varname )
  IF_ERROR_RETURN(status=1)
  call ReadRc( rcF, 'init.c13.multiply', c13_multiplier, status, default=1.0 )
  IF_ERROR_RETURN(status=1)
  ! close rcfile:
  call Done( rcF, status  )
  IF_NOTOK_RETURN(status=1)


  write (gol,*) '.... Opening file with initial d13C ratios for reading: '//trim(filename) ; call goPr

  call Init( hdf, trim(filename), 'read', status )
  IF_NOTOK_RETURN(status=1)

  call Init( sds_c13, hdf, trim(varname), status )  
  IF_NOTOK_RETURN(status=1)


  call GetInfo( sds_c13, status,data_rank=data_rank, data_dims=data_dims )
  IF_NOTOK_RETURN(status=1)

  if(data_rank /= 3.and.data_rank /= 4) then
    write (gol,*) 'Rank of requested dataset ('//trim(varname)//') is not 3 or 4 (time,lon,lat,levels) ...error' ; call goPr
    IF_NOTOK_RETURN(status=1)
  endif

  allocate(c13_file(data_dims(1),data_dims(2),data_dims(3),1) )

  if(data_rank.eq.3) then 
     write (gol,*) '.... Retrieving variable '//trim(varname) ; call goPr
  endif

  if(data_rank.eq.4) then
     write (gol,*) '.... Retrieving first time step of variable '//trim(varname) ; call goPr
  endif

  call ReadData( sds_c13, c13_file, status)
  IF_NOTOK_RETURN(status=1)

  call Done( sds_c13, status )  
  IF_NOTOK_RETURN(status=1)

  write (gol,*) '.... Success' ; call goPr
  call Done( hdf, status )  
  IF_NOTOK_RETURN(status=1)
  write (gol,*) '.... File closed ' ; call goPr

  !WP! Multiply with specified factor

  c13_file = c13_file * c13_multiplier

  !WP Now grid the input file to the dimensions of the current model run

  !IV coarsen_emission does not regrid in the vertical, therefore must be replaced by Fill3D
  !call coarsen_emission('co2_init',data_dims(1),data_dims(2),data_dims(3),co2_file(:,:,:,1),co2_init,avg_field) ! not needed

  im_in=data_dims(1)
  jm_in=data_dims(2)
  lm_in=data_dims(3)

  !WP! We are assuming a 1x1 degree input file for now, this can be changed in the future

  if(data_dims(1) /= 360) then
    write (gol,'("input from initial concentration file implemented for 1 degree longitude only ...")'); call goErr
    TRACEBACK; status=1; return
  endif
  if(data_dims(2) /= 180) then
    write (gol,'("input from initial concentration file implemented for 1 degree latitude only ...")'); call goErr
    TRACEBACK; status=1; return
  endif

  !WP! We are assuming a 1x1 degree input file for now, this can be changed in the future

  xbeg_in = -180.0
  dx_in = 360.0/im_in
  ybeg_in = -90.0
  dy_in = 180.0/jm_in

  call Init( lli_in, xbeg_in+0.5*dx_in, dx_in, im_in, &
                   ybeg_in+0.5*dy_in, dy_in, jm_in, status )

  IF_NOTOK_RETURN(status=1)


  allocate( at_in(0:lm_in) )
  allocate( bt_in(0:lm_in) )

  select case ( lm_in )
      case ( lm(1) )  !WP! If vertical dimension equals current setup, use at, bt from current run
        at_in = at
        bt_in = bt
        levi_in = levi

      case ( 34 )
        at_in = at_tm5_std_34
        bt_in = bt_tm5_std_34
        call Init( levi_in, lm_in, at_in, bt_in, status )
        IF_NOTOK_RETURN(status=1)

       case default
          write (gol,'("input from initial concentration file implemented for 25 or 34 levels only ...")'); call goErr
          TRACEBACK; status=1; return

  end select
  deallocate( at_in )
  deallocate( bt_in )

  !print*,sum(c13_file),im_in,jm_in,lm_in
  
  write (gol,*) '.... Global mean 13C mixing ratio in input field     : ',sum(c13_file)/(im_in*jm_in*lm_in) ; call goPr
  
  
  do region=1,nregions
      m => m_dat(region)%data
      rm => mass_dat(region)%rm_t
      rxm => mass_dat(region)%rxm_t
      rym => mass_dat(region)%rym_t
      rzm => mass_dat(region)%rzm_t
      sp => sp_dat(region)%data


      call Fill3D( lli(region), levi, 'n', sp(1:im(region),1:jm(region),1), c13_init(region)%d3, &
                         lli_in, levi_in, c13_file(:,:,:,1), 'mass-aver', status )

      write (gol,*) '.... Mean 13C mixing ratio in TM5 field region nr:',region,sum(c13_init(region)%d3)/(im(region)*jm(region)*lm(region)) ; call goPr




      
      do n=1,ntracetloc

         ! There is a bug with TM5 using the reduced grid.  Apparently
         ! negative concentrations are not permitted in this scenario.
         ! Thus for the forward simulation, we start each component
         ! with the background C13 field, on which sinks can act.  To
         ! get the component concentration, the bg field must be
         ! subtracted.

         !WP! Get the global tracer number, we only initialize ic13_bg (if we
         !WP! run fwd) or we initialize all c13 tracers (in case of inv)

         offsetn = sum(ntracet_ar(0:myid-1))
         !print*,myid,n,offsetn+n,trim(names(offsetn+n))
         if (goLoCase(trim(names(offsetn+n)))=='co2c13'.or.goLoCase(trim(names(offsetn+n)))=='co2c13_bg') then 

            if (region==1) then
                !print*,'.... Found tracer ',offsetn+n,' with name: ',trim(names(offsetn+n))
                write (gol,'(a,i4,a,a)') '.... Found tracer ',offsetn+n,' with name: ',trim(names(offsetn+n)) ; call goPr
                write (gol,'(a)')      '.... Initializing with C13 from file'               ; call goPr
            endif
            
            rm(1:im(region),1:jm(region),1:lm(region),n) = &
              ((c13_init(region)%d3/1000.)+1.)*0.011112*co2_init(region)%d3*1.e-6*m(1:im(region),1:jm(region),1:lm(region))*xm13co2/xmair
       !  else

        !    rm(:,:,:,:) = 0.0

         endif
           

      enddo

#ifdef slopes
      rxm(:,:,:,:) = 0.0
      rym(:,:,:,:) = 0.0
      rzm(:,:,:,:) = 0.0
#endif
      nullify(m)
      nullify(rm)
      nullify(rxm)
      nullify(rym)
      nullify(rzm)
      nullify(sp)      
  enddo
  deallocate(co2_file)
  deallocate(c13_file)
  

  if (myid==0) write (gol,'(a)')   '.... All other tracers initialized to 0.0! '               ; call goPr

  do region=1,nregions
    deallocate(c13_init(region)%d3)
  enddo

  write (gol,'(a,a)')  rname ,' Done'               ; call goPr
  
  end subroutine start_c13


end module user_input
