!
      module ecmwf_module
!
! Provide daily ecmwf data for lower boundary of t,u,v,z (see lbc.F)
! Data file has 6-hourly data. To use this instead of the averaged
!   daily values, set ecmwf_daily below to false. HOWEVER, if this
!   is done, then gswm perturbations should NOT be added (see lbc.F).
!
      use params_module,only: nlon,nlonp4,glon1,nlat,dlat,dlon,nlevp1
      use nchist_module,only: handle_ncerr,nc_open
      use addfld_module,only: addfld
      implicit none
      integer :: ncid,ecmwf_nlev,ecmwf_ntime 
#include <netcdf.inc>
!
! Year, day of year, and hour from ecmwf file:
      integer,allocatable,dimension(:),save :: ecmwf_year,ecmwf_day,
     |  ecmwf_hour
!
! Output of this module is T,Z,U,V at 10 mb interpolated to current time.
      real,dimension(nlonp4,nlat) :: z_ecmwf, t_ecmwf, u_ecmwf, v_ecmwf
!
! use daily ecmwf am 3/09
      logical, parameter ::  ecmwf_daily = .true.
      contains
!-----------------------------------------------------------------------
      subroutine ecmwf_bndry(ncfile,istep,iyear,iday,isecs)
      implicit none
!
! Args:
      character(len=*),intent(in) :: ncfile
      integer,intent(in) :: 
     |  istep,   ! model time step (1->nstep)
     |  iyear,   ! current model year (4-digit)
     |  iday,    ! current model calendar day (1-367)
     |  isecs    ! current model time in seconds
!
! Local:
      integer :: i
      real :: hour
!
      hour = float(isecs)/3600.
!
! If first step, open and read coord vars of ecmwf data file:
      if (istep==1) then
        write(6,"('ecmwf_bndry: ncfile=',a)") ncfile
        call get_ecmwf_file(ncfile)
      endif
!
! Get daily ECMWF
      if(ecmwf_daily) then
         call get_ecmwf_daily(iyear,iday,hour)
      else	  
! Get ECMWF data at current model time:
         call get_ecmwf_time(iyear,iday,hour)
      endif	 
      end subroutine ecmwf_bndry
!-----------------------------------------------------------------------
      subroutine get_ecmwf_file(ncfile)
      implicit none
!
! Args:
      character(len=*),intent(in) :: ncfile
!
! Local:
      integer :: istat,ndims,idunlim,ngatts,nvars
      integer :: nlat_rd,nlon_rd
      integer :: id_lat,id_lon,id_time,id_lev
      integer :: idv_lev,idv_year,idv_day,idv_hour
      integer :: istart1(1),icount1(1)
      character(len=80) :: dskfile
      real,allocatable :: ecmwf_levs(:)
!
      dskfile = ' '
      call getfile(ncfile,dskfile)
      call nc_open(ncid,ncfile,'OLD','READ')
      istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim)
      write(6,"('get_ecmwf_file: ndims=',i3,' nvars=',i3,' ngatts=',i3,
     |  ' idunlim=',i3)") ndims,nvars,ngatts,idunlim
!
! Get and verify latitude dimension. 
! (read nlat_rd, compare with nlat, which is in params.h):
      istat = nf_inq_dimid(ncid,'lat',id_lat)
      istat = nf_inq_dimlen(ncid,id_lat,nlat_rd)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'get_ecmwf_file: error getting latitude dimension')
      if (nlat_rd /= nlat) then
        write(6,"(/,'>>> get_ecmwf_file: bad nlat_rd=',i3,
     |    ' -- should be nlat=',i3)") nlat_rd,nlat
        call shutdown('get_ecmwf_file')
      endif
!
! Get and verify longitude dimension. 
! (read nlon_rd, compare with nlon, which is in params.h):
      istat = nf_inq_dimid(ncid,'lon',id_lon)
      istat = nf_inq_dimlen(ncid,id_lon,nlon_rd)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'get_ecmwf_file: error getting longitude dimension')
      if (nlon_rd /= nlon) then
        write(6,"(/,'>>> get_ecmwf_file: bad nlon_rd=',i3,
     |    ' -- should be nlon=',i3)") nlon_rd,nlon
        call shutdown('get_ecmwf_file')
      endif
!
! Get levels dimension (should be 1, i.e., only 10 mb level):
      istat = nf_inq_dimid(ncid,'lev',id_lev)
      istat = nf_inq_dimlen(ncid,id_lev,ecmwf_nlev)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'get_ecmwf_file: error getting levels dimension')
!
! Verify that first level on file = 10 mb.
      istat = nf_inq_varid(ncid,"lev",idv_lev) ! day var id
      allocate(ecmwf_levs(ecmwf_nlev),stat=istat)
      istart1(1) = 1
      icount1(1) = ecmwf_nlev
      istat = nf_get_vara_double(ncid,idv_lev,istart1,icount1,
     |  ecmwf_levs)
      write(6,"('get_ecmwf_file: ecmwf_nlev=',i3,' ecmwf_levs(1)=',
     |  f8.2)") ecmwf_nlev,ecmwf_levs(1)         
      if (ecmwf_levs(1) /= 10.) then
        write(6,"('>>> ecmwf_levs(1)=',f8.2,' -- expected 10. mb')")
     |    ecmwf_levs(1)
        call shutdown('get_ecmwf_file')
      endif
      deallocate(ecmwf_levs)
!
! Get number of times (value of unlimited dimension) on the file:
      istat = nf_inq_unlimdim(ncid,id_time) ! id of unlimited record var
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'get_ecmwf_file: error getting unlimited dimension id')
      istat = nf_inq_dimlen(ncid,id_time,ecmwf_ntime) ! length of time var
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,
     |    'Error getting length of day record variable')
        call shutdown('get_ecmwf_file')
      else
!       write(6,"('get_ecmwf_file: ecmwf_ntime=',i5)") ecmwf_ntime
      endif
!
! Read vars year,day_of_year,hour:
!
      allocate(ecmwf_year(ecmwf_ntime),stat=istat)
      allocate(ecmwf_day (ecmwf_ntime),stat=istat)
      allocate(ecmwf_hour(ecmwf_ntime),stat=istat)

!     write(6,"('allocated ecmwf year,day,hour: ecmwf_ntime=',i4,
!    |  ' istat=',i4)") ecmwf_ntime,istat

      istart1(1) = 1
      icount1(1) = ecmwf_ntime
!
! Year:
      istat = nf_inq_varid(ncid,"year",idv_year)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'get_ecmwf_file: error getting id of var year')
      istat = nf_get_vara_int(ncid,idv_year,istart1,icount1,
     |  ecmwf_year)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'get_ecmwf_file: error reading var year')
!     write(6,"('get_ecmwf_file: year=',/,(10i6))") ecmwf_year
!
! Day (of year):
      istat = nf_inq_varid(ncid,"day",idv_day)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'get_ecmwf_file: error getting id of var day (day of year)')
      istat = nf_get_vara_int(ncid,idv_day,istart1,icount1,
     |  ecmwf_day)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'get_ecmwf_file: error reading var day (day of year)')
!     write(6,"('get_ecmwf_file: day of year=',/,(10i6))") ecmwf_day
!
! Hour:
      istat = nf_inq_varid(ncid,"hour",idv_hour)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'get_ecmwf_file: error getting id of var hour')
      istat = nf_get_vara_int(ncid,idv_hour,istart1,icount1,
     |  ecmwf_hour)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'get_ecmwf_file: error reading var hour')
!     write(6,"('get_ecmwf_file: hour=',/,(10i6))") ecmwf_hour

      end subroutine get_ecmwf_file
!-----------------------------------------------------------------------
      subroutine get_ecmwf_time(iyear,iday,hour)
      implicit none
!
! Args:
      integer,intent(in) :: iyear,iday
      real,intent(in) :: hour
!
! Local:
      integer :: i,itecmwf0,itecmwf1,ecmwfhr1
      real :: dummy=0.,hour0,hour1
!
      itecmwf0 = 0
      itecmwf1 = 0
!
! Check first time (if exact match, interp will not be necessary):
      if (ecmwf_year(1)==iyear.and.ecmwf_day(1)==iday.and.
     |  float(ecmwf_hour(1))==hour) then
        itecmwf0 = 1
      endif
!
! Search interior times:
      if (itecmwf0==0) then
        searchloop: do i=1,ecmwf_ntime-1
!
! Check for exact match at i:
          if (ecmwf_year(i)==iyear.and.ecmwf_day(i)==iday.and.
     |      float(ecmwf_hour(i))==hour) then
            itecmwf0 = i
            exit searchloop
          endif
!
! Check for exact match at i+1:
          if (ecmwf_year(i+1)==iyear.and.ecmwf_day(i+1)==iday.and.
     |      float(ecmwf_hour(i+1))==hour) then
            itecmwf0 = i+1 
            exit searchloop
          endif
!
! Check for bracket between i and i+1 (interp will be necessary):
! Note if ecmwf_hour(i) is 18, then ecmwf_hour(i+1) is 0, but need
! to test between 18 and 24 (not 18 and 0), see ecmwfhr1:
          if (itecmwf0==0) then
            ecmwfhr1 = ecmwf_hour(i+1)
            if (ecmwf_hour(i)==18) ecmwfhr1 = 24
            if (ecmwf_year(i)==iyear.and.ecmwf_day(i)==iday.and.
     |        hour > float(ecmwf_hour(i)).and.
     |        hour < float(ecmwfhr1)) then
              itecmwf0 = i
              itecmwf1 = i+1
              exit searchloop
            endif
          endif
        enddo searchloop
      endif
!
! Check for exact match at last time:
      if (itecmwf0==0) then
        if (ecmwf_year(ecmwf_ntime)==iyear.and.
     |      ecmwf_day(ecmwf_ntime)==iday.and.
     |      hour==ecmwf_hour(ecmwf_ntime)) then
          itecmwf0 = ecmwf_ntime 
        endif
      endif
!
! Time not found:
      if (itecmwf0==0) then
        write(6,"('>>> get_ecmwf_time: could not find iyear=',i4,
     |    ' iday=',i3,' hour=',f8.3)") iyear,iday,hour
        call shutdown('get_ecmwf_time')
      endif
!
!     write(6,"(/,'get_ecmwf_time: found iyear,iday,hour=',i5,i3,f8.3,
!    |  ' between itecmwf0,1=',2i4)") iyear,iday,hour,itecmwf0,itecmwf1
!     write(6,"('ecmwf_year,day,hour at itecmwf0=',3i4)")
!    |  ecmwf_year(itecmwf0),ecmwf_day(itecmwf0),ecmwf_hour(itecmwf0)
!     if (itecmwf1 > 0) 
!    |  write(6,"('ecmwf_year,day,hour at itecmwf1=',3i4)")
!    |    ecmwf_year(itecmwf1),ecmwf_day(itecmwf1),ecmwf_hour(itecmwf1)
!
! Read ecmwf data and do interpolation if necessary:
! (4th arg, ecmwf_hour, is ignored if itecmwf1==0)
      if (itecmwf1==0) then ! exact time match at itecmwf0
        hour0 = float(ecmwf_hour(itecmwf0))
        hour1 = dummy
        call get_ecmwf_data(itecmwf0,itecmwf1,hour0,hour1,hour)
      else                 ! will interpolate between it0,it1
        ecmwfhr1 = ecmwf_hour(itecmwf1)
        if (ecmwfhr1==0) ecmwfhr1 = 24
        hour0 = float(ecmwf_hour(itecmwf0))
        hour1 = float(ecmwfhr1)
        call get_ecmwf_data(itecmwf0,itecmwf1,hour0,hour1,hour)
      endif

      end subroutine get_ecmwf_time
!-----------------------------------------------------------------------
      subroutine get_ecmwf_daily(iyear,iday,hour)
      implicit none
!
! Args:
      integer,intent(in) :: iyear,iday
      real,intent(in) :: hour
!
! Local:
      integer :: i,itecmwf0,itecmwf1,ecmwfhr1
      real :: hour0,hour1,hour_int
!
      itecmwf0 = 0
      itecmwf1 = 0
!
! Search interior times for model day:
      if (itecmwf0==0) then
        searchloop: do i=1,ecmwf_ntime
          if (ecmwf_year(i)==iyear.and.ecmwf_day(i)==iday.and.
     |      float(ecmwf_hour(i))==0.) then
            itecmwf0 = i
            exit searchloop
          endif
        enddo searchloop
      endif
!
! check if exact time (12UT) or need day before or day after this model day
! for interpolation (only for first and last day of ECMWF don't interpolate
! ecmwf_ntime-3 since we get the itecmwf0 for 0 UT and the file last time is 18 UT
       if(hour==12.) then
           itecmwf1 = 0  
       elseif(hour > 12. .and. itecmwf0 < ecmwf_ntime-3) then 
           itecmwf1 =  itecmwf0	+ 4
	   hour0    =  0.
	   hour1    = 24.
	   hour_int = hour-12.
	   if(hour_int < 0.or.hour_int >12.) then
	       write(6,*) 'error in subroutine get_ecmwf_daily: ',
     |            'hour_int= ',hour_int
	       stop
	   endif
       elseif(hour <= 12..and. itecmwf0 > 4) then
           itecmwf1 =  itecmwf0	 
           itecmwf0 =  itecmwf0	- 4  
	   hour0    =  0.
	   hour1    = 24.
	   hour_int = hour+12.
	   if(hour_int > 24.or.hour_int<12.) then
	       write(6,*) 'error in subroutine get_ecmwf_daily: ',
     |            'hour_int= ',hour_int
	       stop
	   endif
       endif	       
!
! Time not found:
      if (itecmwf0==0) then
        write(6,"('>>> get_ecmwf_time: could not find iyear=',i4,
     |    ' iday=',i3)") iyear,iday
        call shutdown('get_ecmwf_time')
      endif
!
!     write(6,"(/,'get_ecmwf_time: found iyear,iday,hour=',i5,i3,f8.3,
!    |  ' between itecmwf0,1=',2i4)") iyear,iday,hour,itecmwf0,itecmwf1
!     write(6,"('ecmwf_year,day,hour at itecmwf0=',3i4)")
!    |  ecmwf_year(itecmwf0),ecmwf_day(itecmwf0),ecmwf_hour(itecmwf0)
!     if (itecmwf1 > 0) 
!    |  write(6,"('ecmwf_year,day,hour at itecmwf1=',3i4)")
!    |    ecmwf_year(itecmwf1),ecmwf_day(itecmwf1),ecmwf_hour(itecmwf1)
!
! Read ecmwf data and do daily average
      call get_ecmwf_data_daily(itecmwf0,itecmwf1,
     |    hour0,hour1,hour_int)

      end subroutine get_ecmwf_daily
!-----------------------------------------------------------------------
      subroutine get_ecmwf_data(it0,it1,hour0,hour1,hour)
!
! Read ecmwf data at time it0. If it1 > 0, then read ecmwf data
! at both it0 and it1, and interpolate to hour.
!
      implicit none
!
! Args:
      integer,intent(in) :: it0,it1
      real,intent(in) :: hour0,hour1,hour
!
! Local: 
      integer :: i,j,k
      real,dimension(nlon,nlat) :: tecmwf0,tecmwf1, zecmwf0,zecmwf1,
     |  uecmwf0,uecmwf1, vecmwf0,vecmwf1
      real,dimension(nlonp4,nlevp1) :: tecmwfik,zecmwfik,uecmwfik,
     |  vecmwfik
      real,parameter :: grav=9.8 ! m/s^2
!
! External (util.F):
      real,external :: time_interp
!
! Read data at it0:
      call rd_ecmwf_data(it0,tecmwf0,zecmwf0,uecmwf0,vecmwf0)

!     write(6,"('get_ecmwf_data: it0,1=',2i4,' hour=',f7.3,
!    |  ' tecmwf0 min,max=',2e12.4,' zecmwf0 min,max=',2e12.4)")
!    |  it0,it1,hour,minval(tecmwf0),maxval(tecmwf0),
!    |               minval(zecmwf0),maxval(zecmwf0)
!     write(6,"('get_ecmwf_data: it0,1=',2i4,' hour=',f7.3,
!    |  ' uecmwf0 min,max=',2e12.4,' vecmwf0 min,max=',2e12.4)") 
!    |  it0,it1,hour,minval(uecmwf0),maxval(uecmwf0),
!    |               minval(vecmwf0),maxval(vecmwf0)
!
! If it1 == 0, then no need to interpolate: 
! (exact time match was found at it0):
      if (it1==0) then
        t_ecmwf(3:nlon+2,:) = tecmwf0(:,:)
        z_ecmwf(3:nlon+2,:) = zecmwf0(:,:)
        u_ecmwf(3:nlon+2,:) = uecmwf0(:,:)
        v_ecmwf(3:nlon+2,:) = vecmwf0(:,:)
!
! If it1 > 0, then interpolate to needed hour:
      elseif (it1 > 0) then
        call rd_ecmwf_data(it1,tecmwf1,zecmwf1,uecmwf1,vecmwf1)
!
! Double check interpolates:
        if (hour0 > hour1) then
          write(6,"('>>> get_ecmwf_data: hour0 must be < hour1: hour0=',
     |      f8.3,' hour1=',f8.3)")
          call shutdown('get_ecmwf_data')
        elseif (hour < hour0.or.hour > hour1) then
          write(6,"('>>> get_ecmwf_data: hour must be between hour0,1:',
     |      ' hour0,1=',2f8.3,' hour=',f8.3)") hour0,hour1,hour
          call shutdown('get_ecmwf_data')
        endif
!
! Interpolate to hour:
! Util.F: real function time_interp(f0,f1,time0,time1,time)
!
        do j=1,nlat
          do i=3,nlon+2
            t_ecmwf(i,j) = time_interp(tecmwf0(i-2,j),tecmwf1(i-2,j),
     |        hour0,hour1,hour)
            z_ecmwf(i,j) = time_interp(zecmwf0(i-2,j),zecmwf1(i-2,j),
     |        hour0,hour1,hour)
            u_ecmwf(i,j) = time_interp(uecmwf0(i-2,j),uecmwf1(i-2,j),
     |        hour0,hour1,hour)
            v_ecmwf(i,j) = time_interp(vecmwf0(i-2,j),vecmwf1(i-2,j),
     |        hour0,hour1,hour)
          enddo
        enddo
      endif ! it1 > 0
!
! Periodic points:
! lons 1,2 <- nlon+1,nlon+2, and nlonp+3,nlon+4 <- 3,4
!
      t_ecmwf(1:2,:) = t_ecmwf(nlon+1:nlon+2,:)
      t_ecmwf(nlon+3:nlon+4,:) = t_ecmwf(3:4,:)
      z_ecmwf(1:2,:) = z_ecmwf(nlon+1:nlon+2,:)
      z_ecmwf(nlon+3:nlon+4,:) = z_ecmwf(3:4,:)
      u_ecmwf(1:2,:) = u_ecmwf(nlon+1:nlon+2,:)
      u_ecmwf(nlon+3:nlon+4,:) = u_ecmwf(3:4,:)
      v_ecmwf(1:2,:) = v_ecmwf(nlon+1:nlon+2,:)
      v_ecmwf(nlon+3:nlon+4,:) = v_ecmwf(3:4,:)
!
! real,dimension(nlonp4,nlevp1) :: tecmwfik,zecmwfik
! Save on sec hist redundantly in vertical for tgcmproc_f90:
!
      do j=1,nlat
        do i=3,nlon+2
          tecmwfik(i,:) = t_ecmwf(i,j)
          zecmwfik(i,:) = z_ecmwf(i,j)
          uecmwfik(i,:) = u_ecmwf(i,j)
          vecmwfik(i,:) = v_ecmwf(i,j)
        enddo
        do k=1,nlevp1
          tecmwfik(1:2,k) = tecmwfik(nlon+1:nlon+2,j)
          tecmwfik(nlon+3:nlon+4,k) = tecmwfik(3:4,j)
          zecmwfik(1:2,k) = zecmwfik(nlon+1:nlon+2,j)
          zecmwfik(nlon+3:nlon+4,k) = zecmwfik(3:4,j)
          uecmwfik(1:2,k) = uecmwfik(nlon+1:nlon+2,j)
          uecmwfik(nlon+3:nlon+4,k) = uecmwfik(3:4,j)
          vecmwfik(1:2,k) = vecmwfik(nlon+1:nlon+2,j)
          vecmwfik(nlon+3:nlon+4,k) = vecmwfik(3:4,j)
        enddo
!
        call addfld('ECMWF_T','ECMWF: TN at 10 mb','K',
     |    tecmwfik,'lon',1,nlonp4,'lev',1,nlevp1,j)
        call addfld('ECMWF_Z','ECMWF: Z at 10 mb','cm',
     |    zecmwfik,'lon',1,nlonp4,'lev',1,nlevp1,j)
        call addfld('ECMWF_U','ECMWF: U at 10 mb','cm/s',
     |    uecmwfik,'lon',1,nlonp4,'lev',1,nlevp1,j)
        call addfld('ECMWF_V','ECMWF: V at 10 mb','cm/s',
     |    vecmwfik,'lon',1,nlonp4,'lev',1,nlevp1,j)
      enddo
!
      end subroutine get_ecmwf_data
!-----------------------------------------------------------------------
      subroutine get_ecmwf_data_daily(it0,it1,hour0,hour1,hour)
!
! Read ecmwf data at time it0, it0+1, it0+2, it0+3; then do a daily average
!
      implicit none
!
! Args:
      integer,intent(in) :: it0,it1
      real,intent(in) :: hour0,hour1,hour
!
! Local: 
      integer :: i,j,k
      real,dimension(nlon,nlat) :: tecmwf0,tecmwf1, tecmwf2,tecmwf3,
     |  zecmwf0,zecmwf1,zecmwf2,zecmwf3,
     |  uecmwf0,uecmwf1,uecmwf2,uecmwf3,
     |  vecmwf0,vecmwf1,vecmwf2,vecmwf3,
     |  tecmwf_d0,zecmwf_d0,uecmwf_d0,vecmwf_d0,
     |  tecmwf_d1,zecmwf_d1,uecmwf_d1,vecmwf_d1
      real,dimension(nlonp4,nlevp1) :: tecmwfik,zecmwfik,uecmwfik,
     |  vecmwfik
!
! External (util.F):
      real,external :: time_interp
!
!
! Read data at it0:
      call rd_ecmwf_data(it0,tecmwf0,zecmwf0,uecmwf0,vecmwf0)
      call rd_ecmwf_data(it0+1,tecmwf1,zecmwf1,uecmwf1,vecmwf1)
      call rd_ecmwf_data(it0+2,tecmwf2,zecmwf2,uecmwf2,vecmwf2)
      call rd_ecmwf_data(it0+3,tecmwf3,zecmwf3,uecmwf3,vecmwf3)
!      
!If it1 == 0, then no need to interpolate: 
! (exact time match was found at it0):      
      if (it1==0) then
      do j=1,nlat
        do i=3,nlon+2
            t_ecmwf(i,j) = 0.25*(tecmwf0(i-2,j)+tecmwf1(i-2,j)+
     |        tecmwf2(i-2,j)+tecmwf3(i-2,j))
            z_ecmwf(i,j) = 0.25*(zecmwf0(i-2,j)+zecmwf1(i-2,j)+
     |        zecmwf2(i-2,j)+zecmwf3(i-2,j))
            u_ecmwf(i,j) = 0.25*(uecmwf0(i-2,j)+uecmwf1(i-2,j)+
     |        uecmwf2(i-2,j)+uecmwf3(i-2,j))
            v_ecmwf(i,j) = 0.25*(vecmwf0(i-2,j)+vecmwf1(i-2,j)+
     |        vecmwf2(i-2,j)+vecmwf3(i-2,j))
        enddo
      enddo
!
! If it1 > 0, then interpolate to needed hour:
      elseif (it1 > 0) then
        ! calculate avarge of it0
        do j=1,nlat
          do i=1,nlon
              tecmwf_d0(i,j) = 0.25*(tecmwf0(i,j)+tecmwf1(i,j)+
     |  	tecmwf2(i,j)+tecmwf3(i,j))
              zecmwf_d0(i,j) = 0.25*(zecmwf0(i,j)+zecmwf1(i,j)+
     |  	zecmwf2(i,j)+zecmwf3(i,j))
              uecmwf_d0(i,j) = 0.25*(uecmwf0(i,j)+uecmwf1(i,j)+
     |  	uecmwf2(i,j)+uecmwf3(i,j))
              vecmwf_d0(i,j) = 0.25*(vecmwf0(i,j)+vecmwf1(i,j)+
     |  	vecmwf2(i,j)+vecmwf3(i,j))
          enddo
        enddo
        
        call rd_ecmwf_data(it1  ,tecmwf0,zecmwf0,uecmwf0,vecmwf0)
        call rd_ecmwf_data(it1+1,tecmwf1,zecmwf1,uecmwf1,vecmwf1)
        call rd_ecmwf_data(it1+2,tecmwf2,zecmwf2,uecmwf2,vecmwf2)
        call rd_ecmwf_data(it1+3,tecmwf3,zecmwf3,uecmwf3,vecmwf3)
	
        ! calculate avarge of it1
        do j=1,nlat
          do i=1,nlon
              tecmwf_d1(i,j) = 0.25*(tecmwf0(i,j)+tecmwf1(i,j)+
     |  	tecmwf2(i,j)+tecmwf3(i,j))
              zecmwf_d1(i,j) = 0.25*(zecmwf0(i,j)+zecmwf1(i,j)+
     |  	zecmwf2(i,j)+zecmwf3(i,j))
              uecmwf_d1(i,j) = 0.25*(uecmwf0(i,j)+uecmwf1(i,j)+
     |  	uecmwf2(i,j)+uecmwf3(i,j))
              vecmwf_d1(i,j) = 0.25*(vecmwf0(i,j)+vecmwf1(i,j)+
     |  	vecmwf2(i,j)+vecmwf3(i,j))
          enddo
        enddo
!
! Interpolate to hour:
! Util.F: real function time_interp(f0,f1,time0,time1,time)
        do j=1,nlat
          do i=3,nlon+2
            t_ecmwf(i,j) = time_interp(tecmwf_d0(i-2,j),
     |        tecmwf_d1(i-2,j),hour0,hour1,hour)
            z_ecmwf(i,j) = time_interp(zecmwf_d0(i-2,j),
     |        zecmwf_d1(i-2,j),hour0,hour1,hour)
            u_ecmwf(i,j) = time_interp(uecmwf_d0(i-2,j),
     |        uecmwf_d1(i-2,j),hour0,hour1,hour)
            v_ecmwf(i,j) = time_interp(vecmwf_d0(i-2,j),
     |        vecmwf_d1(i-2,j),hour0,hour1,hour)
          enddo
        enddo
	
      endif
!
! Periodic points:
! lons 1,2 <- nlon+1,nlon+2, and nlonp+3,nlon+4 <- 3,4
!
      t_ecmwf(1:2,:) = t_ecmwf(nlon+1:nlon+2,:)
      t_ecmwf(nlon+3:nlon+4,:) = t_ecmwf(3:4,:)
      z_ecmwf(1:2,:) = z_ecmwf(nlon+1:nlon+2,:)
      z_ecmwf(nlon+3:nlon+4,:) = z_ecmwf(3:4,:)
      u_ecmwf(1:2,:) = u_ecmwf(nlon+1:nlon+2,:)
      u_ecmwf(nlon+3:nlon+4,:) = u_ecmwf(3:4,:)
      v_ecmwf(1:2,:) = v_ecmwf(nlon+1:nlon+2,:)
      v_ecmwf(nlon+3:nlon+4,:) = v_ecmwf(3:4,:)
!
! real,dimension(nlonp4,nlevp1) :: tecmwfik,zecmwfik
! Save on sec hist redundantly in vertical for tgcmproc_f90:
!
      do j=1,nlat
        do i=3,nlon+2
          tecmwfik(i,:) = t_ecmwf(i,j)
          zecmwfik(i,:) = z_ecmwf(i,j)
          uecmwfik(i,:) = u_ecmwf(i,j)
          vecmwfik(i,:) = v_ecmwf(i,j)
        enddo
        do k=1,nlevp1
          tecmwfik(1:2,k) = tecmwfik(nlon+1:nlon+2,j)
          tecmwfik(nlon+3:nlon+4,k) = tecmwfik(3:4,j)
          zecmwfik(1:2,k) = zecmwfik(nlon+1:nlon+2,j)
          zecmwfik(nlon+3:nlon+4,k) = zecmwfik(3:4,j)
          uecmwfik(1:2,k) = uecmwfik(nlon+1:nlon+2,j)
          uecmwfik(nlon+3:nlon+4,k) = uecmwfik(3:4,j)
          vecmwfik(1:2,k) = vecmwfik(nlon+1:nlon+2,j)
          vecmwfik(nlon+3:nlon+4,k) = vecmwfik(3:4,j)
        enddo
!
        call addfld('ECMWF_T','ECMWF: TN at 10 mb','K',
     |    tecmwfik,'lon',1,nlonp4,'lev',1,nlevp1,j)
        call addfld('ECMWF_Z','ECMWF: Z at 10 mb','cm',
     |    zecmwfik,'lon',1,nlonp4,'lev',1,nlevp1,j)
        call addfld('ECMWF_U','ECMWF: U at 10 mb','cm/s',
     |    uecmwfik,'lon',1,nlonp4,'lev',1,nlevp1,j)
        call addfld('ECMWF_V','ECMWF: V at 10 mb','cm/s',
     |    vecmwfik,'lon',1,nlonp4,'lev',1,nlevp1,j)
      enddo
!
      end subroutine get_ecmwf_data_daily
!-----------------------------------------------------------------------
      subroutine rd_ecmwf_data(it,tecmwf,zecmwf,uecmwf,vecmwf)
!
! Read ecmwf T,Z,U,V at time index it.
!
      implicit none
!
! Args:
      integer,intent(in) :: it
      real,dimension(nlon,nlat),intent(out) :: 
     |  tecmwf,zecmwf,uecmwf,vecmwf
!
! Local:
      integer :: istat,start_4d(4),count_4d(4),
     |  idv_t,idv_z,idv_u,idv_v
      character(len=120) :: char120
!
      start_4d(1:3) = 1
      start_4d(4) = it
      count_4d(1) = nlon
      count_4d(2) = nlat
      count_4d(3) = ecmwf_nlev
      count_4d(4) = 1
!
! Read T:
      istat = nf_inq_varid(ncid,"T_GDS0_ISBL",idv_t)
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,
     |    'Error getting var id of ecmwf T_GDS0_ISBL')
        call shutdown('rd_ecmwf_data')
      endif
      istat = nf_get_vara_double(ncid,idv_t,start_4d,count_4d,tecmwf)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_get_vara_double',
     |    ' for field tecmwf: it=',i4)") it
        call handle_ncerr(istat,char120)
      else
!       write(6,"('rd_ecmwf_data: read tecmwf: it=',i4,' min,max=',
!    |    2e12.4)") it,minval(tecmwf),maxval(tecmwf)
      endif
!
! Read Z:
      istat = nf_inq_varid(ncid,"Z_GDS0_ISBL",idv_z)
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,
     |    'Error getting var id of ecmwf Z_GDS0_ISBL')
        call shutdown('rd_ecmwf_data')
      endif
      istat = nf_get_vara_double(ncid,idv_z,start_4d,count_4d,zecmwf)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_get_vara_double',
     +    ' for field zecmwf')")
        call handle_ncerr(istat,char120)
      else
!       write(6,"('rd_ecmwf_data: read zecmwf: it=',i4,' min,max=',
!    |    2e12.4)")  it,minval(zecmwf),maxval(zecmwf)
      endif
!
! Read U:
      istat = nf_inq_varid(ncid,"U_GDS0_ISBL",idv_u)
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,
     |    'Error getting var id of ecmwf U_GDS0_ISBL')
        call shutdown('rd_ecmwf_data')
      endif
      istat = nf_get_vara_double(ncid,idv_u,start_4d,count_4d,uecmwf)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_get_vara_double',
     +    ' for field uecmwf')")
        call handle_ncerr(istat,char120)
      else
!       write(6,"('rd_ecmwf_data: read uecmwf: it=',i4,' min,max=',
!    |    2e12.4)")  it,minval(uecmwf),maxval(uecmwf)
      endif
!
! Read V:
      istat = nf_inq_varid(ncid,"V_GDS0_ISBL",idv_v)
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,
     |    'Error getting var id of ecmwf V_GDS0_ISBL')
        call shutdown('rd_ecmwf_data')
      endif
      istat = nf_get_vara_double(ncid,idv_v,start_4d,count_4d,vecmwf)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_get_vara_double',
     +    ' for field vecmwf')")
        call handle_ncerr(istat,char120)
      else
!       write(6,"('rd_ecmwf_data: read uecmwf: it=',i4,' min,max=',
!    |    2e12.4)")  it,minval(uecmwf),maxval(vecmwf)
      endif
      end subroutine rd_ecmwf_data
!-----------------------------------------------------------------------
      end module ecmwf_module