!
      subroutine readsource(ier)
!
! Read source history.
!
      use params_module,only: nlon,nlonp2,nlonp4,nlevp1,nlat
      use input_module,only: source,tempdir,source_start,output,start,
     |  f107,f107a,power,ctpoten,dynamo,potential_model
      use hist_module,only: nsource,nhist,ioutfile,ncid,iprint,h
      use fields_module,only: tlbc,ulbc,vlbc,tlbc_nm,ulbc_nm,vlbc_nm
      use addfld_module,only: addfld
      use mpi_module,only: lat0,lat1,lon0,lon1
#ifdef MPI
      use mpi_module,only: mp_periodic_f4d,mp_dynpot,
     |  mp_periodic_f2d,mp_bndlons_f2d,mp_bndlats_f2d
#endif
      use fields_module,only: itp,poten,un,f4d,nf4d_hist,i_tn,i_un,
     |  i_vn,dynpot
      implicit none
!
! Arg:
      integer,intent(out) :: ier
!
! Local:
      integer :: lu,nth,ncid_source,j,i,k,lat
      real :: fmin,fmax,fminj,fmaxj
#ifdef MPI
      real :: flbc(lon0:lon1,lat0:lat1,6)
      real :: flbc_2(lon0-2:lon1+2,lat0-2:lat1+2,6)
#endif
!
      ier = 0
!
! If source file was provided, open read-only, and close afterwards.
! If source file not provided, open first output file, and leave it
!   open for possible appending later by output_hist.
!
! Source was provided -- read source history from source file:
      if (nsource==1) then 
        call rdsource(source,tempdir,source_start,ncid,nth,
     |    .false.,iprint)
        nhist  = 0         ! no histories on output file
        ioutfile = 0       ! no output file currently in use
!
! Source file was not provided -- search 1st output file:
      else
        call rdsource(output(1),tempdir,start(:,1),ncid,nth,.true.,
     |    iprint)
        nhist  = nth       ! number of histories on current output file
        ioutfile = 1       ! current output file name is output(ioutfile)
      endif
!
      if (ncid==0) then
        ier = 1
        return
      endif
!
! Do mpi periodic points exchange for f4d(:) 
! Moved here from sub nc_rdhist because mpi necessary when f4d data
!  is allocated only for task-local subdomain block.
! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4
!
#ifdef MPI
      call mp_periodic_f4d(itp)
!
! Periodic points for t,u,v lbc
! These are dimensioned (lond0:lond1,latd0:latd1) (see fields.F):
! real :: flbc(lon0:lon1,lat0:lat1,6)
!  
! am 8/09 set halo points
!    
      flbc_2(:,:,1) = tlbc(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc_2(:,:,2) = ulbc(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc_2(:,:,3) = vlbc(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc_2(:,:,4) = tlbc_nm(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc_2(:,:,5) = ulbc_nm(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc_2(:,:,6) = vlbc_nm(lon0-2:lon1+2,lat0-2:lat1+2)

      call mp_bndlons_f2d(flbc_2,lon0,lon1,lat0,lat1,6)

      tlbc(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,1)
      ulbc(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,2)
      vlbc(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,3)
      tlbc_nm(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,4)
      ulbc_nm(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,5)
      vlbc_nm(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,6)
      
      flbc_2(:,:,1) = tlbc(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc_2(:,:,2) = ulbc(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc_2(:,:,3) = vlbc(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc_2(:,:,4) = tlbc_nm(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc_2(:,:,5) = ulbc_nm(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc_2(:,:,6) = vlbc_nm(lon0-2:lon1+2,lat0-2:lat1+2)

      call mp_bndlats_f2d(flbc_2,lon0,lon1,lat0,lat1,6)

      tlbc(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,1)
      ulbc(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,2)
      vlbc(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,3)
      tlbc_nm(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,4)
      ulbc_nm(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,5)
      vlbc_nm(lon0-2:lon1+2,lat0-2:lat1+2) = flbc_2(:,:,6)
      
      flbc(:,:,1) = tlbc(lon0:lon1,lat0:lat1)
      flbc(:,:,2) = ulbc(lon0:lon1,lat0:lat1)
      flbc(:,:,3) = vlbc(lon0:lon1,lat0:lat1)
      flbc(:,:,4) = tlbc_nm(lon0:lon1,lat0:lat1)
      flbc(:,:,5) = ulbc_nm(lon0:lon1,lat0:lat1)
      flbc(:,:,6) = vlbc_nm(lon0:lon1,lat0:lat1)

      call mp_periodic_f2d(flbc,lon0,lon1,lat0,lat1,6)

      tlbc(lon0:lon1,lat0:lat1) = flbc(:,:,1)
      ulbc(lon0:lon1,lat0:lat1) = flbc(:,:,2)
      vlbc(lon0:lon1,lat0:lat1) = flbc(:,:,3)
      tlbc_nm(lon0:lon1,lat0:lat1) = flbc(:,:,4)
      ulbc_nm(lon0:lon1,lat0:lat1) = flbc(:,:,5)
      vlbc_nm(lon0:lon1,lat0:lat1) = flbc(:,:,6)
!
! Define dynpot(nlonp1,0:nlatp1,nlevp1) from electric potential
! subdomains, which were read from source history. Do this only
! if dynamo input flag is set:
!
      if (dynamo > 0.or.
     |   (dynamo <= 0 .and.trim(potential_model) /= 'NONE')) then
        call mp_dynpot(itp)
      else
        write(6,"('readsource: dynamo=',i2,' --> zero out ',
     |    'electric potential and dynpot...')") dynamo
        poten(:,:,:,itp) = 0. ! zero out electric potential
        dynpot = 0. ! whole array
      endif
#else
      do lat=lat0,lat1
	do i=1,2
	  ulbc(i,lat)	     = ulbc(nlon+i,lat)
	  ulbc(nlonp2+i,lat) = ulbc(i+2,lat)
	  vlbc(i,lat)	     = vlbc(nlon+i,lat)
	  vlbc(nlonp2+i,lat) = vlbc(i+2,lat)
	  ulbc_nm(i,lat)	= ulbc_nm(nlon+i,lat)
	  ulbc_nm(nlonp2+i,lat) = ulbc_nm(i+2,lat)
	  vlbc_nm(i,lat)	= vlbc_nm(nlon+i,lat)
	  vlbc_nm(nlonp2+i,lat) = vlbc_nm(i+2,lat)
	enddo
      enddo
      call set_periodic_f4d(itp)      
!
! Set_dynpot for non-mpi runs:
      call set_dynpot(itp)
#endif
!
! Change some minor species (initial run only):
      if (nsource==1) call fix_minor
!
!     do j=lat0,lat1
!       do i=1,nf4d_hist
!         call fminmax(f4d(i)%data(:,:,j,itp),(lon1-lon0+1)*nlevp1,
!    |      fmin,fmax)
!         write(6,"('readsource: j=',i3,' i=',i3,' f4d(i) min,max=',
!    |      2e12.4)") j,i,fmin,fmax
!       enddo
!     enddo
!
! These addfld calls will not work, unless the following conditional
! is in addfld after the time2write call:
!     if (.not.wrsech.and.istep /= 0) return
!
!     do j=lat0,lat1
!       call addfld('rdsrc_tn',' ',' ',
!    |    f4d(i_tn)%data(:,lon0:lon1,j,itp),'lon',lon0,lon1,
!    |    'lev',1,nlevp1,j)
!       call addfld('rdsrc_un',' ',' ',
!    |    f4d(i_un)%data(:,lon0:lon1,j,itp),'lon',lon0,lon1,
!    |    'lev',1,nlevp1,j)
!       call addfld('rdsrc_vn',' ',' ',
!    |    f4d(i_vn)%data(:,lon0:lon1,j,itp),'lon',lon0,lon1,
!    |    'lev',1,nlevp1,j)
!     enddo
      end subroutine readsource
!-------------------------------------------------------------------
      subroutine rdsource(filepath,tempdir,mtime,ncid,nthist,
     |  reopen_append,iprint)
!
! Acquire and read source history:
!
! On input:
!   filepath = mss path to history file
!   tempdir  = path to a temporary directory
!   mtime(3) = model time of requested source history
!   reopen_append: if true, reopen the file for later writing after
!                  reading the history.
!   iprint: if > 0, report to stdout
!
! On output:
!   ncid   = file id of history file
!   nthist = source history is nth history on the file 
!   global history structure h is defined (see nc_rdhist)
!
      use nchist_module,only: nc_open,nc_close,nc_rdhist
      implicit none
! 
! Args:
      character(len=*),intent(in) :: filepath,tempdir
      integer,intent(in) :: mtime(3),iprint
      integer,intent(out) :: nthist,ncid
      logical,intent(in) :: reopen_append
!
! Local:
      integer ::
     |  mday,mhour,mmin, ! model day,hour,minute from header
     |  j,               ! latitude loop index
     |  ier              ! error flag
      real :: dum,rj
      character(len=120) :: diskfile
!
! Acquire source file:
!     call mkdiskflnm(filepath,diskfile)
!     call getms(filepath,diskfile,tempdir)
      diskfile = ' '
      call getfile(filepath,diskfile)
      write(6,"('Acquired source history file ',a,
     |  /,'  (disk file is ',a,')')") trim(filepath),trim(diskfile)
!
! Open existing netcdf file for read-only:
      call nc_open(ncid,diskfile,'OLD','READ')
      if (ncid==0) then
        write(6,"(/,'>>> rdsource: error opening ',a,' as a ',
     |    'netcdf file.')") trim(diskfile)
!       return
        call shutdown('open netcdf source history file')
      endif
!
! Search for and read the source history:
      write(6,"('rdsource calling nc_rdhist..')")
      call nc_rdhist(ncid,diskfile,mtime,nthist,ier)
      if (ier > 0) then
        write(6,"(/,'>>> ERROR return from nc_rdhist reading')")
        write(6,"('  source file ',a,' from ncid=',i8,' mtime=',
     |    3i4)") trim(diskfile),ncid,mtime
        stop 'nc_rdhist'
      else
        write(6,"('Completed reading SOURCE history.')")
      endif
!
! Optionally initialize AR and NAT to global means:
! 8/04: init_ar_nat not ready yet 
!       (see ~roble/timegcm/tgcm24dinos1/modsrc.yrrunth) 
!     call init_ar_nat
!
      call nc_close(ncid)
!
! Reopen file for writing if necessary:
      if (reopen_append) call nc_open(ncid,diskfile,'OLD','WRITE')
      end subroutine rdsource
!-----------------------------------------------------------------------
      subroutine set_periodic_f4d(itx)
!
! Set periodic points for all f4d fields (serial or non-mpi only):
!
      use params_module,only: nlonp4
      use fields_module,only: f4d,nf4d_hist
      implicit none
      integer,intent(in) :: itx
      integer :: n

! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4
      do n=1,nf4d_hist
        f4d(n)%data(:,1,:,itx) = f4d(n)%data(:,nlonp4-3,:,itx)
        f4d(n)%data(:,2,:,itx) = f4d(n)%data(:,nlonp4-2,:,itx)
!
        f4d(n)%data(:,nlonp4-1,:,itx) = f4d(n)%data(:,3,:,itx)
        f4d(n)%data(:,nlonp4  ,:,itx) = f4d(n)%data(:,4,:,itx)
      enddo
      end subroutine set_periodic_f4d
!-----------------------------------------------------------------------
      subroutine set_dynpot(itx)
!
! This is called for non-MPI runs only:
! Define dynpot(nlonp1,0:nlatp1,nlevp1)  (3d electric potential geographic)
!   from 4d field poten(levd0:levd1,lond0:lond1,latd0:latd1,2) that has
!   been read from source history.
!
      use params_module,only: nlat,nlonp1,nlevp1
      use fields_module,only: poten,dynpot
      implicit none
      integer,intent(in) :: itx
!
! Local:
      integer :: k,i,j
!
      do j=1,nlat
!
! 8/13/07 btf: changed nlonp4 to nlonp1 for dynpot lon dimension
        do i=1,nlonp1
          do k=1,nlevp1
            dynpot(i,j,k) = poten(k,i,j,itx)
          enddo ! k=1,nlevp1
        enddo ! i=1,nlonp1
      enddo ! j=1,nlat
      end subroutine set_dynpot
!-----------------------------------------------------------------------
      subroutine fix_minor
!
! Optionally change minor species read from the source history. 
! Multiply the 3d field by the given factor, if the factor is /= 1.
!
      use fields_module,only: f4d,nf4d_hist
      implicit none
      real,parameter ::
     |  fac_co2 = 1.,
     |  fac_co  = 1.,
     |  fac_h2o = 1.,
     |  fac_h2  = 1.,
     |  fac_hox = 1.,
     |  fac_ch4 = 1.,
     |  fac_n4s = 1.,
     |  fac_noz = 1.
      integer :: i
!
      do i=1,nf4d_hist
        select case (trim(f4d(i)%short_name))
          case('CO2')
            if (fac_co2 /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_co2
              write(6,"('fix_minor: multiplied field ',a,' by fac_co2=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_co2,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('CO2_NM')
            if (fac_co2 /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_co2
              write(6,"('fix_minor: multiplied field ',a,' by fac_co2=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_co2,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('CO')
            if (fac_co /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_co
              write(6,"('fix_minor: multiplied field ',a,' by fac_co=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_co,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('CO_NM')
            if (fac_co /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_co
              write(6,"('fix_minor: multiplied field ',a,' by fac_co=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_co,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('H2O')
            if (fac_h2o /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_h2o
              write(6,"('fix_minor: multiplied field ',a,' by fac_h2o=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_h2o,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('H2O_NM')
            if (fac_h2o /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_h2o
              write(6,"('fix_minor: multiplied field ',a,' by fac_h2o=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_h2o,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('H2')
            if (fac_h2 /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_h2
              write(6,"('fix_minor: multiplied field ',a,' by fac_h2=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_h2,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('H2_NM')
            if (fac_h2 /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_h2
              write(6,"('fix_minor: multiplied field ',a,' by fac_h2=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_h2,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('HOX')
            if (fac_hox /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_hox
              write(6,"('fix_minor: multiplied field ',a,' by fac_hox=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_hox,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('HOX_NM')
            if (fac_hox /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_hox
              write(6,"('fix_minor: multiplied field ',a,' by fac_hox=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_hox,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('CH4')
            if (fac_ch4 /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_ch4
              write(6,"('fix_minor: multiplied field ',a,' by fac_ch4=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_ch4,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('CH4_NM')
            if (fac_ch4 /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_ch4
              write(6,"('fix_minor: multiplied field ',a,' by fac_ch4=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_ch4,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('N4S')
            if (fac_n4s /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_n4s
              write(6,"('fix_minor: multiplied field ',a,' by fac_n4s=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_n4s,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('N4S_NM')
            if (fac_n4s /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_n4s
              write(6,"('fix_minor: multiplied field ',a,' by fac_n4s=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_n4s,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('NOZ')
            if (fac_noz /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_noz
              write(6,"('fix_minor: multiplied field ',a,' by fac_noz=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_noz,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case('NOZ_NM')
            if (fac_noz /= 1.) then
              f4d(i)%data = f4d(i)%data*fac_noz
              write(6,"('fix_minor: multiplied field ',a,' by fac_noz=',
     |          e12.4,' 3d min,max=',2e12.4)") trim(f4d(i)%short_name),
     |          fac_noz,minval(f4d(i)%data),maxval(f4d(i)%data)
            endif
          case default
        end select
      enddo ! i=1,nf4d_hist
      end subroutine fix_minor