subroutine bcgen (mon1, iyr1, monn, iyrn, mon1rd,                &
                  iyr1rd, monnrd, iyrnrd, mon1clm, iyr1clm,      &
                  monnclm, iyrnclm, nlat, nlon, mon1out,         &
                  iyr1out, monnout, iyrnout, ncidin, outfilclim, &
                  outfilamip, nmon, oldttcalc, history)
   use prec
   use types

!------------------------------------------------------------------------
!    21 August 1997
!    Karl E. Taylor
!    PCMDI
!    taylor13@llnl.gov
!
!    14 February 2003
!    Jim Rosinski
!    NCAR
!    Modified to be CAM-specific, and use namelist input
!
!    Documentation provided by Karl Taylor
!
!  Purpose:

!    to create on some specified "target" grid an artificial mid-month 
!    sea ice fraction and SST data set (referred to here as the 
!    "boundary condition data" set) that, when linearly interpolated 
!    (in time), produces the observed monthly means (referred to here
!    as the "observed data").  

!  REFERENCE:  
!       
!    Taylor, K.E., D. Williamson, and F. Zwiers (2000): The sea surface
!          temperature and sea-ice concentration boundary conditions for
!          AMIP II simulations. PCMDI Report No. 60 and UCRL-MI-125597, 
!          Lawrence Livermore National Laboratory, Livermore, CA, 25 pp.
!          
!          pdf file available at 
!                http://www-pcmdi.llnl.gov/pcmdi/pubs/ab60.html

!  Input files:

!      monthly mean (observed) SSTs and sea ice fraction on the target grid

!      fortran namelist defining input and output filenames, and time periods
!      of input and output data.
    
!  Output files:

!         (1) The boundary condition data set for the years requested.
!               (The filenames contain the string "bc" and don't 
!                include the string "clim" (e.g. amipbc_sst_1978.nc).)
!         (2) A climatological boundary condition data set based on the
!               years specified for computing this climatology.
!               (The filenames contain the strings "bc" and "clim"
!               (e.g., amipbc_sic_360x180_1979_2001_clim.nc).)

!  Libraries used:
!      netcdf library: interface to netCDF files

!  Overview of algorithm:

!     The following assumes familiarity with the information contained
!     in the reference listed above.

!        Away from regions of near-freezing temperature, the SSTs are
!     generated by solving a set of N linear equations where N
!     is the number of months considered.  The mid-month (boundary
!     condition) temperatures for a given month depend on the observed
!     monthly-mean temperature of that month, and also the temperatures
!     of the preceeding month and following month.  Thus, 
!     the mid-month temperature for the first and last months 
!     formally require observed temperatures before and after 
!     the time period considered.  In practice the dependence 
!     on the temperatures outside the period considered is fairly weak. 
!     To close the problem mathematically, however, we must impose some
!     sort of condition on the temperatures before and following the
!     period of interest.  

!        If observations are available for several months before and 
!     after the period of interest, then these can be used along with 
!     imposition of a periodic boundary condition to fully constrain the
!     problem.  (By "periodic boundary condition," we mean  that the 
!     first month of the available observations is assumed to follow    
!     the last month).  The periodic boundary condition should only be 
!     imposed on an  observational period that comprises an integral 
!     number of years, so, for example, if the first month is March, the
!     last month must be February.

!        The periodic condition is a convenient way of closing the
!     problem mathematically, but is of course unrealistic.  In  
!     practice, however, it really affects only the mid-month 
!     temperatures of a few months: the temperatures for the 2 or 3 
!     months at the beginning and end of the observational period 
!     considered.  If the observational record extends well outside 
!     (both before and after) the period for which boundary condition 
!     data are needed, then this is not a problem.  If, however, we need
!     boundary condition data for every month for which observations are
!     available, then the months at the beginning and end will be  
!     sensitive to this assumption.

!        If you really want to simulate the full time period for which 
!     observed monthly means are available, the "end" effects
!     can be reduced by generating artifical "observed" temperatures
!     for a brief period preceeding and following the actual 
!     observational period.  In this code we do this as follows: The 
!     artificial "observed" temperature anomalies are assumed to decay  
!     to zero in the first few months prior to and following the  
!     observational period, so that the temperature approaches  
!     climatology.  The rate of decay is based on the the   
!     autocorrelation function for the monthly mean SST time-series    
!     (analyzed for the years 1979-1998).  The spatial mean 
!     (area-weighted over all ocean grid cells) of the autocorrelation 
!     function, at lags of 1, 2, 3, ..., and 8 months have been
!     computed and are specified within this code.  The correlations 
!     for lags greater than 8 months are small, and are assigned so 
!     that the correlation reaches zero (in a smooth way) at month 12. 
!     [Note: for sea ice fraction the correlations are based on 
!     observations for lags less than 6 months, and the correlation is 
!     assumed to be zero for lags greater than 7 months.]

!        In summary, to minimize dependence on the boundary end 
!     conditions (i.e., the values of "observations" specfied for the 
!     period before and after the interval where observations are 
!     actually available), it is best to use only the mid-month values 
!     that are generated for a subinterval away from the boundary.  In  
!     practice this means the mid-month values for the 2 or 3 months at 
!     the beginning and the end of the observational period should be  
!     excluded from use in forcing a GCM simulation.  

!        Consider the following example.  Suppose observations
!     are available for the months January, 1956 through August, 2002.
!     Suppose, further, we would like to prescribe the SST boundary
!     condition in an atmospheric GCM simulation of these same months
!     (1/1956 - 8/2002).  The recommended approach would be to generate
!     artificial observational data for the year preceeding 1/56 and the
!     year following 8/02.  We also require that we treat an integral 
!     number of years.  We could conservatively choose to extend the 
!     "observations" to 12/03 at the end, and to also tack on a year
!     at the beginning (i.e., 12 months starting at 1/55).  
!     Since the artifical "observational" data approaches climatology
!     outside the time period of interest, the user must also 
!     specify which years will contribute to the climatology. If there
!     is no significant temperature trend over the period considered,
!     the climatology could be justifiably based on any interval of
!     several (say, 10 or more) years.  

!        If, however, the data exhibit a noticable trend, then the
!     user should probably avoid trying to simulate the entire 
!     period because the current code is unable to handle this case
!     accurately; it assumes that the artificial data outside both ends
!     of the observational period "decay" toward (i.e., approach) the 
!     same climatology.  It would be better to assume that at the 
!     beginning of the period and at the end of the period, the values 
!     approached the different characteristic climatologies.

!        In any case, it is wise to avoid the last few months for which
!     observations are currently available.  This is because as 
!     observations become available for succeeding months, these
!     will replace the artificially generated "observations" and will
!     have some effect on boundary condition data for the few months 
!     near the end of the observational period.

!     The user can specify the beginning and ending months of the
!        following:
!            
!        1) the interval for which observational monthly mean data
!               are available.
!        2) the entire interval over which the analysis will be
!               performed, including buffer periods prior to and
!               following the observed data, which are recommended to 
!               total a year or two.
!        3) the interval over which the climatology is computed.
!        4) the interval defining which months will be included
!               as output (in the form of SST mid-month boundary
!               conditions).

!        Away from regions of near-freezing temperature, the mid-month
!     SSTs [represented here by S] satisfy:

!               A S = B

!      where S is a column vector of dimension N representing the 
!           mid-month values that constitute the SST boundary condition,
!           B is a column vector of dimension N, representing the 
!           observed monthly mean time-series, and A (except for 2 
!           elements) is a tri-diagnal matrix of dimension NxN.  The
!           two unconforming elements account for the assumption of 
!           periodicity. (The following shows the matrix structure.)
!

!          |b(1) c(1)             ...    a(1) | 
!          |a(2) b(2) c(2)                    |
!          | .   a(3) b(3) c(3)   ...         |
!      A = |                   .              |
!          |                   .              |
!          |                   .              |
!          |       ...    a(N-1) b(N-1) c(N-1)|
!          |c(N)   ...       .    a(N)   b(N) |
 
!      where 
!          a(i) = aa(i)/8
!          b(i) = 1 - aa(i)/8 - cc(i)/8
!          c(i) = cc(i)/8

!          aa(i) = 2*n(i)/(n(i)*n(i-1))
!          cc(i) = 2*n(i)/(n(i)*n(i+1))

!      where n(i) = number of days in month i.

!     Note that if all the months were of equal length, then aa=cc=1 
!          and a=1/8, b=3/4, and c=1/8.

!        For sea ice, the above procedure has to be modified because
!     the sea ice fraction is constrained to be between 0 and 1.  
!     Similarly, the water temperature cannot fall below its freezing
!     point, so this also places a physical constraint that is not 
!     always consistent with the above procedure.  In these cases 
!     the equations that must be solved have a similar structure as 
!     shown above, but the coefficients (a's, b's, and c's) depend on 
!     temperature or the sea ice fraction.  Thus, the equations in this
!     case are no longer linear, but they can be solved using an
!     iterative Newton-Raphson approach.

!     The Jacobian that is required under this approach is 
!     not generated analytically, but is approximated numerically.

!     There is another constraint necessary to ensure a unique
!     solution in the nonlinear case.  To see why the constraint
!     is needed, consider a grid cell where the ocean is ice-free
!     year around, except for 1 month, when the ice fraction is
!     10%.   The mid-month value for this month is not uniquely
!     determined because the cell could be covered by little
!     ice over the entire month, or by lots of ice for only a short
!     time in the middle of the month.  The algorithm relied on 
!     in this code tries to minimize the absolute difference between
!     mid-month values that exceed the maximum physically allowed
!     value (S_max=100% for sea ice) and S_max, while still yielding 
!     the correct monthly means: i.e.,
 
!                 where S > S_max, minimize (S - S_max) 
         
!     Similarly, for values that are less than the minimum physically
!     allowed value (S_min = 0% for sea ice):

!                 where S < S_min, minimize (S_min - S)

!     In the example above this results in the following mid-month 
!     values for sea ice (assuming for simplicity that all months
!     are of the same length):    S(i-1) = S(i+1) = -20%  and 
!     S(i) = 20%  where i is the month with mean sea ice fraction of 
!     10%.  These mid-month values when linearly interpolated give a 
!     sea ice fraction for month i that starts at 0%, linearly grows to
!     20% at the middle of the month, and then linearly decreases to 0%
!     at the end of the month.  For the month preceeding and the month
!     following month i, the linear interpolation leads to negative
!     values, but recall that the model's alogrithm will "clip" these
!     values, setting them to 0%.

!  Namelist details:
!
!    SPECIFY the first and last month and year for period in which 
!        observed monthly mean data will be read. (The first month
!        read must not preceed mon1, iyr1, and the last month must
!        not follow monn, iyrn).  For example:
           
!      mon1rd = 1            !AMIP
!      iyr1rd = 1956         !AMIP

!      monnrd = 8            !AMIP
!      iyrnrd = 2002         !AMIP

!    SPECIFY first and last month and year for entire period that will
!        be treated (i.e., the period of interest plus buffers at ends).
!        Note that the entire period treated should be an integral   
!        number of years.  For example:

!      mon1 = 1            !AMIP
!      iyr1 = 1955         !AMIP

!      monn = 12           !AMIP
!      iyrn = 2003         !AMIP

!     SPECIFY the first and last month and year that will be included in
!         climatological mean.  (This must be an interval within the 
!         observed period).  For example:

!      mon1clm = 1            !AMIP
!      iyr1clm = 1979         !AMIP

!      monnclm = 12           !AMIP
!      iyrnclm = 2000         !AMIP

!     SPECIFY the first and last month and year written to the output 
!         file.  (Try to begin at east a few months after the mon1rd,
!         iyr1rd, and end a few months before monnrd, iyrnrd, to avoid
!         sensitivity to the artificial data outside the observed 
!         period.)  For example:
 
!      mon1out = 1            !AMIP
!      iyr1out = 1956         !AMIP
!      
!      monnout = 6            !AMIP
!      iyrnout = 2002         !AMIP
!               
!    SPECIFY input dataset (as output by the program regrid) on the
!        target grid.  For example:
!
!      infil = 'regrid.T42.nc'
!
!    SPECIFY output climatological boundary condition dataset name.
!        For example:
!
!      outfilclim = 'outfilclim.T42.nc'
!               
!    SPECIFY output AMIP-style boundary condition dataset name
!
!      outfilclim = 'outfilclim.T42.nc'
!------------------------------------------------------------------------

   implicit none

   include 'netcdf.inc'
!
! Input arguments
!
   integer, intent(in) :: mon1        ! start month of period of interest plus buffer 
   integer, intent(in) :: iyr1        ! start year of period of interest plus buffer  
   integer, intent(in) :: monn        ! end month of period of interest plus buffer   
   integer, intent(in) :: iyrn        ! end year of period of interest plus buffer    
   integer, intent(in) :: mon1rd      ! start month to read from input data           
   integer, intent(in) :: iyr1rd      ! start year to read from input data            
   integer, intent(in) :: monnrd      ! end month to read from input data             
   integer, intent(in) :: iyrnrd      ! end year to read from input data              
   integer, intent(in) :: mon1clm     ! start month to use for climatology            
   integer, intent(in) :: iyr1clm     ! start year to use for climatology             
   integer, intent(in) :: monnclm     ! end month to use for climatology              
   integer, intent(in) :: iyrnclm     ! end year to use for climatology               
   integer, intent(in) :: nlat        ! number of latitudes (e.g. 64 for typical T42)
   integer, intent(in) :: nlon        ! number of longitudes (e.g. 128 for typical T42)
   integer, intent(in) :: mon1out     ! start month written to output file (default mon1rd)
   integer, intent(in) :: iyr1out     ! start year written to output file (default iyr1rd)
   integer, intent(in) :: monnout     ! end month written to output file (default monnrd)
   integer, intent(in) :: iyrnout     ! end year written to output file (default iyrnrd)
   integer, intent(in) :: ncidin      ! input file netcdf id
   integer, intent(in) :: nmon        ! total number of months of period of interest plus buffer

   logical, intent(in) :: oldttcalc   ! for bfb agreement with original code

   character*120, intent(in) :: outfilclim  ! output filename for climatology
   character*120, intent(in) :: outfilamip  ! output filename for year-by-year output
   character(len=*), intent(in) :: history  ! history attribute
!
! Local workspace
!
   integer :: i, j             ! lon, lat indices
   integer :: k                ! index value 1=ice, 2=sst
   integer :: m, mp, n, mm, l  ! month indices
   integer :: m1, m2           ! month indices
   integer :: ismc             ! number of monthly means smoothed (climatology)
   integer :: jsmc             ! number of grid cells affected by smoothing (climatology)
   integer :: ismf             ! number of year-by-year monthly means smoothed (year-by-year)
   integer :: jsmf             ! number of grid cells affected by smoothing (year-by-year)
   integer :: isea             ! number of ocean grid cells
   integer :: nmonout          ! number of output months
   integer :: prvyr, prvmo     ! previous year, month index from input file
   integer :: yr, mo           ! year, month index from input file

   logical :: dateok           ! flag indicates acceptable input date format
!
! netcdf info
!
   integer :: ncidclim = -1             ! output filehandle (climatology)
   integer :: ncidamip = -1             ! output filehandle (year-by-year)
   integer :: varid = -1                ! variable id
   integer :: lonid = -1                ! longitude id (input)
   integer :: latid = -1                ! latitude id (input)
   integer :: dateidclim = -1           ! date id (climatology)
   integer :: dateidamip = -1           ! date id (year-by-year)
   integer :: datesecidclim = -1        ! seconds of date id (climatology)
   integer :: datesecidamip = -1        ! seconds of date id (year-by-year)
   integer :: timeidclim = -1           ! time id (clim)
   integer :: timeidamip = -1           ! time id (year-by-year)
   integer :: dateidin = -1             ! date info on input file
   integer :: date(1) = -1              ! input date (yyyymmdd)
   integer :: start(3) = (/-1,-1,-1/)   ! for reading/writing non-contiguous data in netcdf file
   integer :: start1d(1) = -1           ! date start index on input file
   integer :: start1dtst(1) = -1        ! test date index on input file
   integer :: kount(3) = (/-1,-1,-1/)   ! for reading/writing non-contiguous data in netcdf file

   integer, parameter :: maxiter = 100  ! max iterations in solving time-diddling

   real(r8), parameter :: fac = 1.0     ! factor to multiply input data by (before any cropping).
   real(r8), parameter :: bbmin = 0.001 ! smallest diagonal element allowed in Jacobian
                                        ! i.e. will change 0.0005 to 0.001
   real(r8) :: obsclim(nlon,nlat,12)    ! climatological values
   real(r8) :: obsamip(nlon,nmon)       ! year-by-year values
   real(r8) :: arr(nlon,nmon)           ! sst or ice conc. (1 latitude)
   real(r8) :: centmon(nlon,nlat,12)    ! sst or ice conc. climatology (all latitudes)
   real(r8) :: a(nmon)                  ! lower diagonal matrix coeffs (year-by-year)
   real(r8) :: c(nmon)                  ! upper diagonal matrix coeffs (year-by-year)
   real(r8) :: ac(12)                   ! lower diagonal matrix coeffs (climatology)
   real(r8) :: cc(12)                   ! upper diagonal matrix coeffs (climatology)
   real(r8) :: xlon(nlon)               ! longitudes on input file (between 0 and 360)
   real(r8) :: xlat(nlat)               ! latitudes on input file (between -90 and 90)

   integer :: monlen(12,2)              ! length of months

   data monlen/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, &
               31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/

   character(len=80) :: center = 'NCAR' ! research center creating data

   type(icesstparms) :: icesst(2)       ! struct containing variable-specific info
!
! Set ice parms. dt value: Not much smoothing => 0->100 becomes 2->98 between months
!
   icesst(1) = icesstparms ('ICEFRAC', -1, &                                ! inname, invarid

                            'ice_cov', &                                    ! climname
                            'ice_cov_prediddle', &                          ! climname_pre
                            'BCS Pseudo Sea-ice concentration', &           ! climlongname
                            'Sea-ice concentration before time diddling', & ! climlongname_pre

                            'ice_cov', &                                    ! amipname
                            'ice_cov_prediddle', &                          ! amipname_pre
                            'BCS Pseudo Sea-ice concentration', &           ! amiplongname
                            'Sea-ice concentration before time diddling', & ! amiplongname_pre

                            'fraction', &                                   ! units

                            -1, -1, -1, -1, &                               ! varids
                       
                            1.e-4,  0.,  1., 1.e-4, 0.02, &      ! conv, tmin, tmax, varmin, dt
                  (/0.53, 0.23, 0.13, 0.08, 0.05, 0.02, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/)) ! correl
!
! With Karl's bugfix to tt calc, convergence criterion for ice can be tightened.
! Leave as is for now since halving this number can lead to big excursions in ice 
! concentration.
!
   if (oldttcalc) then
      icesst(1)%dt = (icesst(1)%tmax - icesst(1)%tmin) / 50.
   else
      icesst(1)%dt = (icesst(1)%tmax - icesst(1)%tmin) / 100.
      icesst(1)%dt = (icesst(1)%tmax - icesst(1)%tmin) / 50.
   end if
!
! Set sst parms.  Cropping SST at 1000. effectively means do not crop
!
   icesst(2) = icesstparms ('SST', -1, &                                    ! inname, invarid  
                                                                                               
                            'SST_cpl', &                                    ! climname         
                            'SST_cpl_prediddle', &                          ! climname_pre     
                            'BCS Pseudo SST', &                             ! climlongname     
                            'SST before time diddling', &                   ! climlongname_pre 
                                                                                               
                            'SST_cpl', &                                    ! amipname         
                            'SST_cpl_prediddle', &                          ! amipname_pre     
                            'BCS Pseudo SST', &                             ! amiplongname     
                            'SST before time diddling', &                   ! amiplongname_pre 
                                                                                               
                            'deg_C', &                                      ! units            
                                                                                               
                            -1, -1, -1, -1, &                               ! varids           
                                                                                               
                            1.e-3, -1.8, 1000., 0., 0.001, &      ! conv,tmin, tmax, varmin, dt
                   (/0.68, 0.46, 0.33, 0.26, 0.20, 0.17, 0.14, 0.11, 0.08, 0.05, 0.02, 0.00/)) ! correl
             
! *****************************************************************

   monlen(2,2) = 28       ! JR for no leap year

!
! Copy dimension vars from input to output
!
   call wrap_nf_inq_varid (ncidin, 'lon', lonid)
   call wrap_nf_inq_varid (ncidin, 'lat', latid)
   call wrap_nf_get_var_double (ncidin, lonid, xlon)
   call wrap_nf_get_var_double (ncidin, latid, xlat)
!
! Open output files and define metadata
!
   call setup_outfile (ncidclim, outfilclim, 'CLIM', dateidclim, datesecidclim, &
                       timeidclim, icesst, nlon, nlat, xlon,                    &
                       xlat, 0, 1, history)
   call setup_outfile (ncidamip, outfilamip, 'AMIP', dateidamip, datesecidamip, &
                       timeidamip, icesst, nlon, nlat, xlon,                    &
                       xlat, iyr1out, mon1out, history)
!
! Determine starting index of data to be read.  Assume this and subsequent data are monthly
! values, centered in the middle of the month.
!
   start1d(1) = 1
   call wrap_nf_inq_varid (ncidin, 'date', dateidin)
   do while (.true.)
      call wrap_nf_get_vara_int (ncidin, dateidin, start1d, 1, date)
      if (date(1)/10000 == iyr1rd .and. mod(date(1)/100,100) == mon1rd) then
         write(6,*)'Input data at iyr1rd, mon1rd=', iyr1rd, mon1rd, ' are record # ', start1d(1)
         exit
      end if
      start1d(1) = start1d(1) + 1
   end do
!
! Ensure that enough data are present, and that they increment properly in time
!
   start1dtst(1) = start1d(1)
   call wrap_nf_get_vara_int (ncidin, dateidin, start1dtst, 1, date)
   prvyr = date(1)/10000
   prvmo = mod(date(1)/100,100)

   do while (.true.)
      start1dtst(1) = start1dtst(1) + 1
      call wrap_nf_get_vara_int (ncidin, dateidin, start1dtst, 1, date)
      yr = date(1)/10000
      mo = mod(date(1)/100,100)
!
! Exit loop when date read in extend beyond requested data
!
      if (yr > iyrnrd) exit
      if (yr == iyrnrd .and. mo >= monnrd) then
         write(6,*)'Input yr,mo=', yr, mo,' will be the last date read from input file'
         exit
      end if
      dateok = (yr == prvyr .and. (mo == prvmo+1 .and. mo >= 1 .and. mo <= 12)) .or. &
               yr == prvyr+1 .and. mo == 1
      if (.not. dateok) then
         write(6,*)'Bad date info read: yr,mo=',yr,mo
         stop 999
      end if
      write(6,*)'Input yr,mo=', yr, mo, ' OK'
      prvyr = yr
      prvmo = mo
   end do
!
! Set kount array for netcdf calls
!
   kount(1) = nlon
   kount(2) = 1
   kount(3) = 1
!
! Loop 1 to 2 to do ice first then sst
!
   do 22 k=1,2
      write(6,*) 'Output from monthly boundary condition generator '
      write(6,*) 'clim var name = ', icesst(k)%climname
      write(6,*) 'climlong name = ', icesst(k)%climlongname
      write(6,*) 'amip var name = ', icesst(k)%amipname
      write(6,*) 'amiplong name = ', icesst(k)%amiplongname
      write(6,*) 'units = ', icesst(k)%units
      write(6,*) 'center = ', center
      write(6,*) 'nlon = ', nlon
      write(6,*) 'nlat = ', nlat
      write(6,*) 'nmon = ', nmon
      write(6,*) 'mon1rd = ', mon1rd
      write(6,*) 'iyr1rd = ', iyr1rd
      write(6,*) 'monnrd = ', monnrd
      write(6,*) 'iyrnrd = ', iyrnrd
      write(6,*) 'mon1clm = ', mon1clm
      write(6,*) 'iyr1clm = ', iyr1clm
      write(6,*) 'monnclm = ', monnclm
      write(6,*) 'iyrnclm = ', iyrnclm
      write(6,*) 'mon1   = ', mon1
      write(6,*) 'iyr1   = ', iyr1
      write(6,*) 'monn   = ', monn
      write(6,*) 'iyrn   = ', iyrn
      write(6,*) 'mon1out = ', mon1out
      write(6,*) 'iyr1out = ', iyr1out
      write(6,*) 'monnout = ', monnout
      write(6,*) 'iyrnout = ', iyrnout
      write(6,'(a,1pe14.7)') 'conv = ', icesst(k)%conv
      write(6,*) 'maxiter = ', maxiter
      write(6,'(a,1pe14.7)') 'bbmin = ', bbmin
      write(6,'(a,1pe14.7)') 'varmin = ', icesst(k)%varmin
      write(6,'(a,1p,12(e10.3))') 'correl = ', icesst(k)%correl
      write(6,*) 'monlen = ', monlen
      write(6,'(a,1pe14.7)') 'tmin = ', icesst(k)%tmin
      write(6,'(a,1pe14.7)') 'tmax = ', icesst(k)%tmax
      write(6,'(a,1pe14.7)') 'dt = ', icesst(k)%dt
      write(6,'(a,1pe14.7)') 'fac =  ', fac
      write(6,*) 'oldttcalc=', oldttcalc

! Calculate jacobian elements for climatological years

      m = 12
      mp = 1
      do n=1,12
         mm = m
         m = mp
         mp = mp + 1
         if (mp == (12+1)) mp = 1
         ac(n) = 2.0_r8*monlen(m,1)/(monlen(m,1) + monlen(mm,1))
         cc(n) = 2.0_r8*monlen(m,1)/(monlen(m,1) + monlen(mp,1))
      end do

! Calculate jacobian elements for all years

      i = iyr1
      j = 1
      if (((mod(i,4) == 0) .and. (mod(i,100) /= 0)) .or. (mod(i,400) == 0)) j = 2

      m = mod((mon1+12-2), 12) + 1
      mp = mod((mon1-1), 12) + 1

      do n=1,nmon
         mm = m 
         m  = mp
         mp = mp + 1
         if (mp == (12+1)) then
            mp = 1
            i = i + 1
            j = 1
            if (((mod(i,4) == 0) .and. (mod(i,100) /= 0)) .or. (mod(i,400) == 0)) j = 2
         end if
         a(n) = 2.0_r8*monlen(m,j)/(monlen(m,j) + monlen(mm,j))
         c(n) = 2.0_r8*monlen(m,j)/(monlen(m,j) + monlen(mp,j))
      end do
!
! Get input variable id
!
      call wrap_nf_inq_varid (ncidin, icesst(k)%inname, icesst(k)%invarid)
!
! Monthly mean computation section: save memory by computing and writing one 
! latitude slice at a time
!
      ismc = 0
      jsmc = 0
      ismf = 0
      jsmf = 0
      isea = 0

      do j=1,nlat
         m1 = (iyr1rd-iyr1)*12 + mon1rd - mon1 + 1
         m2 = (iyrnrd-iyr1)*12 + monnrd - mon1 + 1
!
! Read in all time values of the j'th latitude slice
!
         start(1) = 1
         start(2) = j
         start(3) = start1d(1)
!
! Note data read in to array index starting at m1 not 1 to allow for padding at start of
! period.
!
         do l=m1,m2
            call wrap_nf_get_vara_double (ncidin, icesst(k)%invarid, start, kount, arr(1,l))
            start(3) = start(3) + 1
         end do
!
! Save off obs for writing to output file
!
         obsamip(:nlon,1:m1-1) = 0.
         obsamip(:nlon,m1:m2) = arr(:nlon,m1:m2)
         obsamip(:nlon,m2+1:) = 0.
         
         write(6,*) 'latitude ', j, ' read successfully.'

! Compute climatological means

         call calcclim (j, nlon, nlat, nmon, ismc,                    &
                        jsmc, mon1clm, iyr1clm, iyr1, mon1,           & 
                        iyr1rd, mon1rd, iyrnrd, monnrd, monnclm,      &
                        iyrnclm, maxiter, icesst(k)%tmin,             &
                           icesst(k)%tmax, icesst(k)%dt,              &
                        icesst(k)%conv, bbmin, isea, obsclim, arr,    &
                        centmon, ac, cc)

! Solve for full mid-month values

         call calcfull (j, nlon, nlat, nmon, ismf,                    &
                        jsmf, iyr1rd, mon1rd, iyr1, mon1,             &
                        iyrn, monn, iyrnrd, monnrd, maxiter,          &
                        icesst(k)%tmin, icesst(k)%tmax, icesst(k)%dt, &
                           icesst(k)%conv, bbmin,                     &
                        obsclim, arr, centmon, a, c,                  &
                        icesst(k)%correl, oldttcalc)
!
! Add climatology back on to output data
!         
         call finish (nlon, nlat, mon1, nmon, centmon, &
                      arr, j)
!
! Write 1 latitude of pre- and post-time-diddled variable to output file.
! 
         start(1) = 1
         start(2) = j
         nmonout = (iyrnout - iyr1out)*12 + (monnout - mon1out + 1)
!
! l is the output time index for the netcdf file
! m indexes into output arrays with a buffer offset
!
         m = (iyr1out - iyr1)*12 + (mon1out - mon1)
         do l=1,nmonout
            m = m + 1
            start(3) = l
            varid = icesst(k)%varidamip
            call wrap_nf_put_vara_double (ncidamip, varid, start, kount, arr(1,m))
            varid = icesst(k)%varidamip_pre
            call wrap_nf_put_vara_double (ncidamip, varid, start, kount, obsamip(1,m))
         end do

         do l=1,12
            start(3) = l
            varid = icesst(k)%varidclim
            call wrap_nf_put_vara_double (ncidclim, varid, start, kount, centmon(1,j,l))
            varid = icesst(k)%varidclim_pre
            call wrap_nf_put_vara_double (ncidclim, varid, start, kount, obsclim(1,j,l))
         end do
      end do     ! j=1,nlat

      write(6,*) ' '
      write(6,*) ismc, ' climatological monthly means were smoothed'
      write(6,*) jsmc, ' grid cells were affected'
      write(6,*) ' '
      write(6,*) ismf, ' monthly means were smoothed'
      write(6,*) jsmf, ' grid cells were affected'
      write(6,*) ' '

22 continue     ! 1=ice 2=sst
!
! Write model-readable date and datesec info to output files
!
   call output_dateinfo (ncidclim, dateidclim, datesecidclim, timeidclim, 0, &
                         1, 0, 12, monlen(1,1))
   call output_dateinfo (ncidamip, dateidamip, datesecidamip, timeidamip, iyr1out, &
                         mon1out, iyrnout, monnout, monlen(1,1))

   call wrap_nf_close (ncidclim)
   call wrap_nf_close (ncidamip)
   call wrap_nf_close (ncidin)

   return
end subroutine bcgen

subroutine finish (nlon, nlat, mon1, nmon, centmon, &
                   arr, j)
!
! Add climatology back on to year-by-year data
!
   use prec
   implicit none

   integer, intent(in) :: nlon, nlat   ! number of lons, lats
   integer, intent(in) :: mon1         ! starting month index of arr
   integer, intent(in) :: nmon         ! total number of months
   integer, intent(in) :: j            ! latitude index

   real(r8), intent(in) :: centmon(nlon,nlat,12)  ! climatology
   real(r8), intent(inout) :: arr(nlon,nmon)      ! year-by-year data

   integer :: i     ! longitude index
   integer :: m     ! month index cycling 1-12
   integer :: n     ! month index for entire period

! Calculate full mid-month time-series

   do n=1,nmon
      m = mod((n+mon1-2), 12) + 1
      do i=1,nlon
         arr(i,n) = arr(i,n) + centmon(i,j,m)
      end do
   end do

   return
end subroutine finish