program makezones

!Take the site data and aggregate by time over the spinup period (1901-1920)
!  and "chop" it up into 24 longitundinal zones
use netcdf
implicit none

integer ng
integer v, n, i, z, y,m,j
character(len=4) yst
character(len=4) mst
character(len=4) zst
character(len=50) metvars(7)
character(len=200) fname
character(len=100) met_prefix
double precision data_in(7,1,1,744)
double precision temp
integer*2 data_zone(7,1,744,1)
integer*2 temp_zone(744,1)
double precision longxy(1,1), latixy(1,1)
integer count_zone(24), ncid_out(24)
integer startyear, endyear
integer starti(3), counti(3), dimid(2), starti_out
integer ierr, ncid, varid, varids_out(24,50), ndaysm(12)
integer mask(720,360)
double precision dtime(1460)
real add_offsets(7), scale_factors(7), data_ranges(7,2)

data ndaysm /31,28,31,30,31,30,31,31,30,31,30,31/

ng=1
startyear = 1999
endyear   = 2006

data_ranges(:,1) =  1.0e9
data_ranges(:,2) = -1.0e9    
print*, 'getting scale factors'

metvars(3)='TBOT'
metvars(6)='PSRF'
metvars(2)='FSDS'
metvars(5)='FLDS'
metvars(4)='RH'
metvars(1)='PRECTmms'
metvars(7)='WIND'

do y=startyear,endyear
  do m=1,12
    write(yst,'(I4)') y
    write(mst,'(I4)') 1000+m
    fname = yst // '-' // mst(3:4) // '.nc'
    ierr = nf90_open(trim(fname), NF90_NOWRITE, ncid)
    do v=1,7
       ierr = nf90_inq_varid(ncid, trim(metvars(v)), varid)
       starti(1:3)  = 1
       counti(1)  =   1
       counti(2)  =   1
       counti(3)  =   ndaysm(m)*24    !hourly
       ierr = nf90_get_var(ncid, varid, data_in(v,1:1,1:1,1:counti(3)), &
                             starti, counti)
       if (maxval(data_in(v,1,1,1:counti(3))) .gt. data_ranges(v,2)) then 
         data_ranges(v,2) = maxval(data_in(v,1,1,1:counti(3)))
       end if
       if (minval(data_in(v,1,1,1:counti(3))) .lt. data_ranges(v,1)) then
         data_ranges(v,1) = minval(data_in(v,1,1,1:counti(3)))
       end if
   end do
   ierr = nf90_close(ncid)
  end do
end do

do v=1,7
  temp = data_ranges(v,2)-data_ranges(v,1)
  data_ranges(v,2) = data_ranges(v,2)+0.02*temp
  data_ranges(v,1) = data_ranges(v,1)-0.02*temp
  add_offsets(v) = (data_ranges(v,2)+data_ranges(v,1))/2.
  scale_factors(v) = (data_ranges(v,2)-data_ranges(v,1))/2**16
  write(*,'(A12,1x,3(g13.6,1x))') trim(metvars(v)), data_ranges(v,:), scale_factors(v)
end do

starti_out=1
z=1
do y=startyear,endyear
  do m=1,12
    write(mst,'(I4)') 1000+m
    print*, y, m
    count_zone(:)=0
    write(yst,'(I4)') y
    fname = yst // '-' // mst(3:4) // '.nc'
    ierr = nf90_open(trim(fname), NF90_NOWRITE, ncid)
    ierr = nf90_inq_varid(ncid, 'LONGXY', varid)
    ierr = nf90_get_var(ncid, varid, longxy)
    ierr = nf90_inq_varid(ncid, 'LATIXY', varid)
    ierr = nf90_get_var(ncid, varid, latixy)
    do v=1,7
       ierr = nf90_inq_varid(ncid, trim(metvars(v)), varid)
       starti(1:3)  = 1
       counti(1)  =   1
       counti(2)  =   1
       counti(3)  =   ndaysm(m)*24    !hourly
       ierr = nf90_get_var(ncid, varid, data_in(v,1:1,1:1,1:counti(3)), &
                             starti, counti)
       data_zone(v,1,1:ndaysm(m)*24,1) = &
            nint((data_in(v,1,1,1:ndaysm(m)*24)-add_offsets(v))/ &
            scale_factors(v)) 
    end do
    ierr = nf90_close(ncid)
    
    if (y .eq. startyear .and. m .eq. 1) then 
       fname = './all_hourly_test.nc'
       ierr = nf90_create(trim(fname), NF90_CLOBBER, ncid_out(z))
       ierr = nf90_def_dim(ncid_out(z), 'n', 1, dimid(2))
       ierr = nf90_def_dim(ncid_out(z), 'DTIME', (endyear-startyear+1)*8760, dimid(1))
       ierr = nf90_def_var(ncid_out(z), 'DTIME', NF90_DOUBLE, dimid(1), &
            varids_out(z,1))
       ierr = nf90_put_att(ncid_out(z), varids_out(z,1), 'long_name', &
            'Day of Year')
       ierr = nf90_put_att(ncid_out(z), varids_out(z,1), 'units', &
            'Days since ' // yst // '-01-01 00:00')
       ierr = nf90_def_var(ncid_out(z), 'LONGXY', NF90_FLOAT, dimid(2), &
            varids_out(z,2))
       ierr = nf90_def_var(ncid_out(z), 'LATIXY', NF90_FLOAT, dimid(2), &
            varids_out(z,3))
       do v=1,7
          ierr = nf90_def_var(ncid_out(z), trim(metvars(v)), NF90_SHORT, &
               dimid(1:2), varids_out(z,4+v))
          ierr = nf90_put_att(ncid_out(z), varids_out(z,4+v), 'add_offset', &
               add_offsets(v))
          ierr = nf90_put_att(ncid_out(z), varids_out(z,4+v), 'scale_factor', &
               scale_factors(v))
       end do
    end if
    ierr = nf90_enddef(ncid_out(z))
    ierr = nf90_put_var(ncid_out(z), varids_out(z,2), longxy)
    ierr = nf90_put_var(ncid_out(z), varids_out(z,3), latixy)

    do i=1,ndaysm(m)*24
       dtime(i) = (starti_out+i-1)/24.0-(0.5/24)
    end do
    starti(1) = starti_out
    counti(1) = ndaysm(m)*24
    ierr = nf90_put_var(ncid_out(z), varids_out(z,1), &
         dtime(1:counti(1)), starti(1:1), counti(1:1))
    starti(2)=1
    starti(1)=starti_out
    counti(2)=1
    counti(1)=ndaysm(m)*24.
    do v=1,7
       temp_zone = data_zone(v,z,:,:)
       ierr = nf90_put_var(ncid_out(z), varids_out(z,4+v), &
            temp_zone(1:counti(1), 1:1), starti(1:2), counti(1:2))
    end do
    starti_out = starti_out+ndaysm(m)*24 
 end do
end do
ierr = nf90_close(ncid_out(z))



end program makezones
