#include "options.inc"

c
c--------------------------------------------------------
c      templates to read surface wind stress (tau) and
c      surface "friction velocity" (ustar) data
c      from external files.
c
c     Notes:
c      "friction velocity" is needed for turbulence closure 
c      models and parameterisations in spflame.
c      wind stress should be in  dyn/cm**2
c      friction velocity should be (cm/s)**3
c
c      Four subroutines are defined:
c      "read_dims_tau" and "read_tau" both called
c      in prep_sbc (in prep_main.F) and
c      "read_dims_ustar" and read_ustar called
c      in prep_sbc as well.
c      Header of subroutines should not be touched.
c      CPP directives are set in prep_template.F
c
c                            c.eden  nov 2002
c--------------------------------------------------------
c
#ifdef read_tau_ecmwf
c
c--------------------------------------------------------
c    use a 3 year monthly mean climatology from the european
c    ECMWF model (same data used in the DYNAMO project)
c--------------------------------------------------------
c
      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt,io
      character*32 stamp
      real (kind=4) :: dpm
      integer (kind=4) nx_,ny_
      print*,' reading monthly wind stress data from ECMWF'
      open(10,file='data/ecmwf/ecmwf_BB.tau',form='unformatted')
      read (10); read (10) stamp, dpm, nx_, ny_
      nx=nx_;ny=ny_; nt=12
      close(10)
      end subroutine read_dims_tau


      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt)

      integer io,n,k,kk,i,j
      real cleanmsk(nx,ny)
      character*32 stamp,s
      real (kind=4) dpm(nt), x4(nx),y4(ny),u(nx,ny),v(nx,ny),
     &               mask4(nx,ny)
      integer (kind=4) k4
      toffset=0.0
      open(10,file='data/ecmwf/ecmwf_BB.tau',form='unformatted')
      do n=1,nt
        read (10); read (10) s, dpm(n), k4, k4, k4, x4,y4,u,v
        time(n)=dpm(n)
        tau(:,:,n,1)=u
        tau(:,:,n,2)=v
        x=x4;y=y4;
      enddo
      close(10)
c   read land mask
      open(10,file='data/ecmwf/ecmwf_BB.mask',form='unformatted')
      read (10) ; read (10); read (10) mask4
      mask=mask4
      close(10)
      cleanmsk = mask
      do j=2, ny-1
        do i=2, nx-1
          if (cleanmsk(i,j)/=0. .and. 
     &        cleanmsk(i+1,j)==0. .and. cleanmsk(i-1,j)==0. .and. 
     &        cleanmsk(i,j+1)==0. .and. cleanmsk(i,j-1)==0.)
     &        mask(i,j) = 0.
        end do
      end do
c   convert N/m**2 to dyn/cm**2 
      tau= 10.*tau
c   apply mask
      do n=1,nt
        where (mask==0.) tau(:,:,n,1) = spval
        where (mask==0.) tau(:,:,n,2) = spval
      enddo
      end subroutine read_tau
c
c
#elif defined read_tau_ncep_monthly
c
c--------------------------------------------------------
c    use the 3 year monthly mean ECMWF climatology plus
c    wind stress anomalies from NCEP/NCAR data.
c    Same forcing was used in Eden+Willebrand, 01, J.Clim.
c--------------------------------------------------------
c
      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt,io
      character*32 stamp
      real (kind=4) :: dpm
      integer (kind=4) nx_,ny_
#include "netcdf.inc"
       integer :: iret, ncid,dims(4),vid

      print*,' reading monthly wind stress data from ECMWF'
      print*,' plus wind stress anomalies from NCEP/NCAR'
      open(10,file='data/ecmwf/ecmwf_BB.tau',form='unformatted')
      read (10); read (10) stamp, dpm, nx_, ny_
      nx=nx_;ny=ny_; 
      close(10)
      iret=nf_open('data/ncep/taux_anom.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'uflx',vid)
      iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt)
      call ncclos (ncid, iret)
      end subroutine read_dims_tau


      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      use time_type_module
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt),mtime(12)
      integer io,n,k,kk,i,j
      real cleanmsk(nx,ny)
      character*32 stamp,s
      real :: mtau(nx,ny,12,2) 
      real (kind=4) dpm(nt), x4(nx),y4(ny),u(nx,ny),v(nx,ny),
     &               mask4(nx,ny)
      integer (kind=4) k4
#include "netcdf.inc"
      integer :: ncid,txid,tyid,dims(4),lonid,latid,iret,ncid2
      integer :: nid
      character*80 name
      real (kind=4) ::  spv
      real (kind=4), allocatable :: x2(:),y2(:),tau2(:,:,:),tx(:)
      real (kind=4), allocatable :: regx(:,:),regy(:,:),nmask(:,:)
      real, allocatable   :: work_t(:)
      integer :: init_t=1, lenw_t,nx2,ny2,nn
      type(time_type) :: t1

      open(10,file='data/ecmwf/ecmwf_BB.tau',form='unformatted')
      do n=1,12
        read (10); read (10) s, dpm(n), k4, k4, k4, x4,y4,u,v
        mtime(n)=dpm(n)
        mtau(:,:,n,1)=u
        mtau(:,:,n,2)=v
        x=x4;y=y4;
      enddo
      close(10)
c   read land mask
      open(10,file='data/ecmwf/ecmwf_BB.mask',form='unformatted')
      read (10) ; read (10); read (10) mask4
      mask=mask4
      close(10)
      cleanmsk = mask
      do j=2, ny-1
        do i=2, nx-1
          if (cleanmsk(i,j)/=0. .and. 
     &        cleanmsk(i+1,j)==0. .and. cleanmsk(i-1,j)==0. .and. 
     &        cleanmsk(i,j+1)==0. .and. cleanmsk(i,j-1)==0.)
     &        mask(i,j) = 0.
        end do
      end do
c   convert N/m**2 to dyn/cm**2 and apply mask
      do n=1,12
        mtau(:,:,n,1)= 10.*mtau(:,:,n,1)
        mtau(:,:,n,2)= 10.*mtau(:,:,n,2)
        where (mask==0.) mtau(:,:,n,1) = spval
        where (mask==0.) mtau(:,:,n,2) = spval
      enddo
c
c   read land/sea mask for NCEP data
c
      iret=nf_open('data/ncep/lsmask.19294.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'lsmask',txid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx2)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny2)
      allocate( nmask(nx2,ny2) )
      iret= nf_get_vara_real (ncid,txid ,(/1,1,1/),(/nx2,ny2,1/),nmask)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      call ncclos (ncid, iret)
c
c   read NCEP/NCAR wind stress anomalies
c
      iret=nf_open('data/ncep/taux_anom.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_open('data/ncep/tauy_anom.cdf',NF_nowrite,ncid2)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'uflx',txid)
      iret=nf_inq_varid(ncid2,'vflx',tyid)
      iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx2)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny2)
      nx2=nx2+2 ! append some rows
      ny2=ny2+4

      allocate(x2(nx2),y2(ny2), tau2(nx2,ny2,2), tx(nx2),
     &         regx(nx,ny), regy(nx,ny) )
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
      iret=NF_INQ_VARID (NCID, name, lonid)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
      iret=NF_INQ_VARID (NCID, name, latid)

      iret= nf_get_vara_real (ncid,lonid ,1,nx2-2,x2(2:nx2-1))
      x2(nx2)=2*x2(nx2-1)-x2(nx2-2)
      x2(1)=x2(2)-(x2(3)-x2(2))
      iret= nf_get_vara_real (ncid,latid ,1,ny2-4,y2(3:ny2-2)) 
      y2(ny2-1)=2*y2(ny2-2)-y2(ny2-3)
      y2(2)=y2(3)-(y2(4)-y2(3))
      y2(ny2)=2*y2(ny2-1)-y2(ny2-2)
      y2(1)=y2(2)-(y2(3)-y2(2))
      y2=y2(ny2:1:-1)
      iret = nf_get_att_real(ncid, txid, 'missing_value', spv)

c   read each month and interpolate on ECMWF grid
      print*,' interpolating data '
      lenw_t = 10*max(nx,ny)
      allocate( work_t(lenw_t) )
      t1 =  set_date(1948, 1, 1, 0, 0, 0)
      toffset=t1%days
      do n=1,nt
        tau2=spv
        iret= nf_get_vara_real (ncid,txid ,
     &         (/1,1,n/),(/nx2-2,ny2-4,1/),tau2(2:nx2-1,3:ny2-2,1))
        iret= nf_get_vara_real (ncid2,tyid ,
     &         (/1,1,n/),(/nx2-2,ny2-4,1/),tau2(2:nx2-1,3:ny2-2,2))

        where( nmask == -1) tau2(2:nx2-1,3:ny2-2,1) = spv
        where( nmask == -1) tau2(2:nx2-1,3:ny2-2,2) = spv
c     acount for appended rows and switch latitudinal ordering
        tau2(nx2,:,:)=tau2(2,:,:)
        tau2(1,:,:)  =tau2(nx2-1,:,:)
        tau2(:,:,:) = tau2(:,ny2:1:-1,:)
c    convert special values
        where( tau2 == spv) tau2=spval
!      convert N/m**2 to dyn/cm**2  and change sign
        where( tau2 /= spval) tau2=tau2*(-10.)
c     interpolate
        call tranlon (tau2(1,1,1), nx2, nx2-2, ny2, x2, x,nx, tx)
        call tranlon (tau2(1,1,2), nx2, nx2-2, ny2, x2, x,nx, tx)
        call ctf (tau2(1,1,1), nx2, ny2, tx, y2, 
     &                  regx(1,1),nx, ny, 1, nx, 1, ny,
     &                  x, y, init_t, work_t, lenw_t,spval)
        init_t=0 ! we can use the same interpolation weights
        call ctf (tau2(1,1,2), nx2, ny2, tx, y2, 
     &                  regy(1,1),nx, ny, 1, nx, 1, ny,
     &                  x, y, init_t, work_t, lenw_t,spval)
c     add mean
        nn = mod(n-1,12)+1
        where( regx == spval) regx=0.
        where( regy == spval) regy=0.
        tau(:,:,n,1)=mtau(:,:,nn,1)+regx(:,:)
        tau(:,:,n,2)=mtau(:,:,nn,2)+regy(:,:)
        where (mask==0.) tau(:,:,n,1) = spval
        where (mask==0.) tau(:,:,n,2) = spval
        time(n)=mtime(nn)
c     care about time
      enddo
      call ncclos (ncid, iret)
      call ncclos (ncid2, iret)
      print*,' done '
      deallocate(x2,y2,tx,tau2,regx,regy, work_t,nmask)
      end subroutine read_tau
c
#elif defined read_tau_nao_proxy
c
c--------------------------------------------------------
c    use the 3 year monthly mean ECMWF climatology plus
c    NAO related time series of wind stress anomalies
c    (which is given by the product of the monthly NAO-index
c    and the monthly regression patterns of taux/tauy with
c    the NAO-index). 
c    Same forcing was used in Eden+Jung 01, J.Clim.
c
c    NOTE: Start and end year of the NAO time series to be used
c    have to be set in the common block tau_nao_reg.
c    The NAO-index is given for 1865 to 2001 a subset of this
c    period can be set with the common block tau_nao_reg
c    See also prep_template_sflx.F 
c--------------------------------------------------------
c
      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt,io
      character*32 stamp
      real (kind=4) :: dpm
      integer (kind=4) nx_,ny_
#include "netcdf.inc"
      integer :: ncid,iret,nid,dims
       integer :: start_y,end_y
       common / tau_nao_reg/ start_y,end_y

      print*,' reading monthly wind stress data from ECMWF'
      print*,' plus NAO related wind stress anomalies from NCEP/NCAR'
      open(10,file='data/ecmwf/ecmwf_BB.tau',form='unformatted')
      read (10); read (10) stamp, dpm, nx_, ny_
      nx=nx_;ny=ny_; nt=12
      close(10)
      iret=nf_open('data/NAO_reg/nao_monthly.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'nao',nid)
      iret=NF_INQ_VARDIMID (NCID,nid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs, nt)
      call ncclos (ncid, iret)
c     start/end_y gives starting and end year of NAO time series
        start_y = max(1865,1865); end_y = min(1948,nt/12+1864) ! 
c        start_y = max(1947,1865); end_y = min(5948,nt/12+1864) ! 
      nt = (end_y-start_y+1)*12
      print*,' from Jan. of year ',start_y
      print*,' to Dec. of year   ',end_y
      print*,' nt=',nt
      end subroutine read_dims_tau

      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      use time_type_module
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt),mtime(12)
      integer io,n,k,kk,i,j
      real cleanmsk(nx,ny)
      character*32 stamp,s
      real :: mtau(nx,ny,12,2) 
      real (kind=4) dpm(nt), x4(nx),y4(ny),u(nx,ny),v(nx,ny),
     &               mask4(nx,ny)
      integer (kind=4) k4
#include "netcdf.inc"
      integer :: ncid,txid,tyid,dims(4),lonid,latid,iret
      integer :: nid
      character*80 name
      real (kind=4) ::  spv,add_offset,scale_factor
      real (kind=4), allocatable :: x2(:),y2(:),tau2(:,:,:),tx(:)
      real (kind=4), allocatable :: regx(:,:,:),regy(:,:,:)
      real, allocatable   :: work_t(:)
      integer :: init_t=1, lenw_t,nx2,ny2,nn
      type(time_type) :: t1
      integer :: start_y,end_y
      common / tau_nao_reg/ start_y,end_y

      open(10,file='data/ecmwf/ecmwf_BB.tau',form='unformatted')
      do n=1,12
        read (10); read (10) s, dpm(n), k4, k4, k4, x4,y4,u,v
        mtime(n)=dpm(n)
        mtau(:,:,n,1)=u
        mtau(:,:,n,2)=v
        x=x4;y=y4;
      enddo
      close(10)
c   read land mask
      open(10,file='data/ecmwf/ecmwf_BB.mask',form='unformatted')
      read (10) ; read (10); read (10) mask4
      mask=mask4
      close(10)
      cleanmsk = mask
      do j=2, ny-1
        do i=2, nx-1
          if (cleanmsk(i,j)/=0. .and. 
     &        cleanmsk(i+1,j)==0. .and. cleanmsk(i-1,j)==0. .and. 
     &        cleanmsk(i,j+1)==0. .and. cleanmsk(i,j-1)==0.)
     &        mask(i,j) = 0.
        end do
      end do
c   convert N/m**2 to dyn/cm**2 and apply mask
      do n=1,12
        mtau(:,:,n,1)= 10.*mtau(:,:,n,1)
        mtau(:,:,n,2)= 10.*mtau(:,:,n,2)
        where (mask==0.) mtau(:,:,n,1) = spval
        where (mask==0.) mtau(:,:,n,2) = spval
      enddo
c
c   read NAO monthly regression patterns
c
      iret=nf_open('data/NAO_reg/reg_tau.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'taux_reg',txid)
      iret=nf_inq_varid(ncid,'tauy_reg',tyid)
      iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx2)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny2)
      nx2=nx2+2 ! append some rows
      ny2=ny2+4
      print*,' nx2 = ',nx2, 'ny2=',ny2

      allocate(x2(nx2),y2(ny2), tau2(nx2,ny2,2), tx(nx2),
     &         regx(nx,ny,12), regy(nx,ny,12) )
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
      iret=NF_INQ_VARID (NCID, name, lonid)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
      iret=NF_INQ_VARID (NCID, name, latid)

      iret= nf_get_vara_real (ncid,lonid ,1,nx2-2,x2(2:nx2-1))
      x2(nx2)=2*x2(nx2-1)-x2(nx2-2)
      x2(1)=x2(2)-(x2(3)-x2(2))
      iret= nf_get_vara_real (ncid,latid ,1,ny2-4,y2(3:ny2-2)) 
      y2(ny2-1)=2*y2(ny2-2)-y2(ny2-3)
      y2(2)=y2(3)-(y2(4)-y2(3))
      y2(ny2)=2*y2(ny2-1)-y2(ny2-2)
      y2(1)=y2(2)-(y2(3)-y2(2))
      y2=y2(ny2:1:-1)
      iret = nf_get_att_real(ncid, txid, 'missing_value', spv)
      add_offset=0.0
      iret= nf_get_att_real (ncid,txid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,txid ,'scale_factor',scale_factor)

c   read each month and interpolate on ECMWF grid
      print*,' interpolating data '
      lenw_t = 10*max(nx,ny)
      allocate( work_t(lenw_t) )
      do n=1,12
        tau2=spv
        iret= nf_get_vara_real (ncid,txid ,
     &         (/1,1,n/),(/nx2-2,ny2-4,1/),tau2(2:nx2-1,3:ny2-2,1))
        iret= nf_get_vara_real (ncid,tyid ,
     &         (/1,1,n/),(/nx2-2,ny2-4,1/),tau2(2:nx2-1,3:ny2-2,2))
        tau2(nx2,:,:)=tau2(2,:,:)
        tau2(1,:,:)  =tau2(nx2-1,:,:)
        tau2(:,:,:) = tau2(:,ny2:1:-1,:)
        where( tau2 == spv) tau2=spval
!      convert N/m**2 to dyn/cm**2  and change sign
        where( tau2 /= spval) tau2=(tau2*scale_factor+add_offset)*(-10.)
        call tranlon (tau2(1,1,1), nx2, nx2-2, ny2, x2, x,nx, tx)
        call tranlon (tau2(1,1,2), nx2, nx2-2, ny2, x2, x,nx, tx)
        call ctf (tau2(1,1,1), nx2, ny2, tx, y2, 
     &                  regx(1,1,n),nx, ny, 1, nx, 1, ny,
     &                  x, y, init_t, work_t, lenw_t,spval)
        init_t=0 ! we can use the same interpolation weights
        call ctf (tau2(1,1,2), nx2, ny2, tx, y2, 
     &                  regy(1,1,n),nx, ny, 1, nx, 1, ny,
     &                  x, y, init_t, work_t, lenw_t,spval)
      enddo
      call ncclos (ncid, iret)
      print*,' done '

c   read NAO index, multiply with regression pattern and add mean
      iret=nf_open('data/NAO_reg/nao_monthly.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'nao',nid)
      t1 =  set_date(start_y, 1, 1, 0, 0, 0)
      toffset=t1%days
      do n=1,nt
        nn = mod(n-1,12)+1
        time(n)=mtime(nn)
        iret= nf_get_vara_real (ncid,nid ,n+(start_y-1865)*12,1,spv)
c        print*,' n=',n,nn,' nao(n)=',spv
        tau(:,:,n,1)=mtau(:,:,nn,1)+spv*regx(:,:,nn)
        tau(:,:,n,2)=mtau(:,:,nn,2)+spv*regy(:,:,nn)
        where (mask==0.) tau(:,:,n,1) = spval
        where (mask==0.) tau(:,:,n,2) = spval
      enddo
      call ncclos (ncid, iret)
      deallocate(x2,y2,tx,tau2,regx,regy, work_t)
      end subroutine read_tau
c
c
#elif defined read_tau_nao_ideal
c
c--------------------------------------------------------
c    use the 3 year monthly mean ECMWF climatology plus
c    NAO related wind stress anomalies (NAO index of +3)
c    (which is given by 
c    the monthly regression patterns of taux/tauy with
c    the NAO-index). 
c    Same forcing was used in Eden+Willebrand 01, J.Clim.
c    See also prep_template_sflx.F 
c--------------------------------------------------------
c
      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt,io
      character*32 stamp
      real (kind=4) :: dpm
      integer (kind=4) nx_,ny_

      print*,' reading monthly wind stress data from ECMWF'
      print*,' plus contant NAO related wind stress anomalies '
      open(10,file='data/ecmwf/ecmwf_BB.tau',form='unformatted')
      read (10); read (10) stamp, dpm, nx_, ny_
      nx=nx_;ny=ny_; nt=12
      close(10)
      print*,' nt=',nt
      end subroutine read_dims_tau

      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      use time_type_module
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt),mtime(12)
      integer io,n,k,kk,i,j
      real cleanmsk(nx,ny)
      character*32 stamp,s
      real :: mtau(nx,ny,12,2) 
      real (kind=4) dpm(nt), x4(nx),y4(ny),u(nx,ny),v(nx,ny),
     &               mask4(nx,ny)
      integer (kind=4) k4
#include "netcdf.inc"
      integer :: ncid,txid,tyid,dims(4),lonid,latid,iret
      integer :: nid
      character*80 name
      real (kind=4) ::  spv,add_offset,scale_factor
      real (kind=4), allocatable :: x2(:),y2(:),tau2(:,:,:),tx(:)
      real (kind=4), allocatable :: regx(:,:,:),regy(:,:,:)
      real, allocatable   :: work_t(:)
      integer :: init_t=1, lenw_t,nx2,ny2,nn
      type(time_type) :: t1

      open(10,file='data/ecmwf/ecmwf_BB.tau',form='unformatted')
      do n=1,12
        read (10); read (10) s, dpm(n), k4, k4, k4, x4,y4,u,v
        mtime(n)=dpm(n)
        mtau(:,:,n,1)=u
        mtau(:,:,n,2)=v
        x=x4;y=y4;
      enddo
      close(10)
c   read land mask
      open(10,file='data/ecmwf/ecmwf_BB.mask',form='unformatted')
      read (10) ; read (10); read (10) mask4
      mask=mask4
      close(10)
      cleanmsk = mask
      do j=2, ny-1
        do i=2, nx-1
          if (cleanmsk(i,j)/=0. .and. 
     &        cleanmsk(i+1,j)==0. .and. cleanmsk(i-1,j)==0. .and. 
     &        cleanmsk(i,j+1)==0. .and. cleanmsk(i,j-1)==0.)
     &        mask(i,j) = 0.
        end do
      end do
c   convert N/m**2 to dyn/cm**2 and apply mask
      do n=1,12
        mtau(:,:,n,1)= 10.*mtau(:,:,n,1)
        mtau(:,:,n,2)= 10.*mtau(:,:,n,2)
        where (mask==0.) mtau(:,:,n,1) = spval
        where (mask==0.) mtau(:,:,n,2) = spval
      enddo
c
c   read NAO monthly regression patterns
c
      iret=nf_open('data/NAO_reg/reg_tau.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'taux_reg',txid)
      iret=nf_inq_varid(ncid,'tauy_reg',tyid)
      iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx2)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny2)
      nx2=nx2+2 ! append some rows
      ny2=ny2+4
      print*,' nx2 = ',nx2, 'ny2=',ny2

      allocate(x2(nx2),y2(ny2), tau2(nx2,ny2,2), tx(nx2),
     &         regx(nx,ny,12), regy(nx,ny,12) )
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
      iret=NF_INQ_VARID (NCID, name, lonid)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
      iret=NF_INQ_VARID (NCID, name, latid)

      iret= nf_get_vara_real (ncid,lonid ,1,nx2-2,x2(2:nx2-1))
      x2(nx2)=2*x2(nx2-1)-x2(nx2-2)
      x2(1)=x2(2)-(x2(3)-x2(2))
      iret= nf_get_vara_real (ncid,latid ,1,ny2-4,y2(3:ny2-2)) 
      y2(ny2-1)=2*y2(ny2-2)-y2(ny2-3)
      y2(2)=y2(3)-(y2(4)-y2(3))
      y2(ny2)=2*y2(ny2-1)-y2(ny2-2)
      y2(1)=y2(2)-(y2(3)-y2(2))
      y2=y2(ny2:1:-1)
      iret = nf_get_att_real(ncid, txid, 'missing_value', spv)
      add_offset=0.0
      iret= nf_get_att_real (ncid,txid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,txid ,'scale_factor',scale_factor)

c   read each month and interpolate on ECMWF grid
      print*,' interpolating data '
      lenw_t = 10*max(nx,ny)
      allocate( work_t(lenw_t) )
      do n=1,12
        tau2=spv
        iret= nf_get_vara_real (ncid,txid ,
     &         (/1,1,n/),(/nx2-2,ny2-4,1/),tau2(2:nx2-1,3:ny2-2,1))
        iret= nf_get_vara_real (ncid,tyid ,
     &         (/1,1,n/),(/nx2-2,ny2-4,1/),tau2(2:nx2-1,3:ny2-2,2))
        tau2(nx2,:,:)=tau2(2,:,:)
        tau2(1,:,:)  =tau2(nx2-1,:,:)
        tau2(:,:,:) = tau2(:,ny2:1:-1,:)
        where( tau2 == spv) tau2=spval
!      convert N/m**2 to dyn/cm**2  and change sign
        where( tau2 /= spval) tau2=(tau2*scale_factor+add_offset)*(-10.)
        call tranlon (tau2(1,1,1), nx2, nx2-2, ny2, x2, x,nx, tx)
        call tranlon (tau2(1,1,2), nx2, nx2-2, ny2, x2, x,nx, tx)
        call ctf (tau2(1,1,1), nx2, ny2, tx, y2, 
     &                  regx(1,1,n),nx, ny, 1, nx, 1, ny,
     &                  x, y, init_t, work_t, lenw_t,spval)
        init_t=0 ! we can use the same interpolation weights
        call ctf (tau2(1,1,2), nx2, ny2, tx, y2, 
     &                  regy(1,1,n),nx, ny, 1, nx, 1, ny,
     &                  x, y, init_t, work_t, lenw_t,spval)
      enddo
      call ncclos (ncid, iret)
      print*,' done '

c   multiply 3 with regression pattern and add mean
      t1 =  set_date(1999, 1, 1, 0, 0, 0)
      toffset=t1%days
      do n=1,nt
        nn = mod(n-1,12)+1
        print*,' n=',n,nn
        time(n)=mtime(nn)
        tau(:,:,n,1)=mtau(:,:,nn,1)+3.0*regx(:,:,nn)
        tau(:,:,n,2)=mtau(:,:,nn,2)+3.0*regy(:,:,nn)
        where (mask==0.) tau(:,:,n,1) = spval
        where (mask==0.) tau(:,:,n,2) = spval
      enddo
      deallocate(x2,y2,tx,tau2,regx,regy, work_t)
      end subroutine read_tau
c
c
#elif defined read_tau_ncep_daily
c
c--------------------------------------------------------
c     use daily wind stress data from the
c     NCEP/NCAR reanalysis project
c     (two parameter have to be set: start_y and end_y
c      denoting the start and end year taken from the
c      reanalysis data, e.g. start_y=1990, end_y=1998)
c--------------------------------------------------------

      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4),nn,k
      character (len=80) :: file
      integer :: start_y,end_y        
      common / tau_ncep/ start_y,end_y
      start_y = 1998
      end_y   = 1999
      print*,' reading daily wind stress data from NCEP '
      print*,' for years ',start_y,' to ',end_y
      nt=0
      do k=start_y,end_y
       write(file,'("data/ncep_daily/uflx.sfc.gauss.",i4,".nc")') k
       iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=nf_inq_varid(ncid,'uflx',varid)
       iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
       iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
       iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
       iret=NF_INQ_DIMLEN (NCID, DIMs(3), nn)
       nt=nt+nn
       call ncclos (ncid, iret)
      enddo
      end subroutine read_dims_tau


      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      use time_type_module
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt)
#include "netcdf.inc"
      integer :: iret,ncid,varid,corner(4),edges(4),k,dims(4),nn
      integer :: n1,i,j,n,xid,yid
      character (len=80) :: file
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:,:,:)
      real (kind=4) :: add_offset,scale_factor
#else
      real , allocatable :: v4(:,:,:)
      real  :: add_offset,scale_factor
#endif
      type(time_type) :: t1,t2,t3
      integer :: start_y,end_y
      common / tau_ncep/ start_y,end_y
c   read land mask
      iret=nf_open('data/ncep_daily/lsmask.19294.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'lsmask',varid)
      corner=(/1,1,1,1/); edges = (/nx,ny,1,1/)
      allocate( v4(nx,ny,1) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      do j=1,ny; do i=1,nx
         mask(i,j)=v4(i,ny-j+1,1)
      enddo; enddo
      deallocate(v4)
      ! read grid info
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),file) 
      iret=nf_inq_varid(ncid,file,xid)
      allocate(v4(nx,1,1) )
      iret= nf_get_vara_real (ncid,xid ,1,nx,v4); x=v4(1:nx,1,1)
      deallocate(v4)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),file) 
      iret=nf_inq_varid(ncid,file,yid)
      allocate(v4(ny,1,1) )
      iret= nf_get_vara_real (ncid,yid ,1,ny,v4); 
      do j=1,ny
         y(j)=v4(ny-j+1,1,1)
      enddo
      deallocate(v4)
      call ncclos (ncid, iret)
      ! time vector
      t1 =  set_date(start_y, 1, 1, 0, 0, 0)
      toffset=t1%days
      n1=1
      do k=start_y,end_y
       write(file,'("data/ncep_daily/uflx.sfc.gauss.",i4,".nc")') k
       iret=nf_open(file,NF_nowrite,ncid)
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=nf_inq_varid(ncid,'uflx',varid)
       iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
       iret=NF_INQ_DIMLEN (NCID, DIMs(3), nn)
       call ncclos (ncid, iret)
       do n=n1,n1+nn-1
         time(n) =  365.0/(nn*1.0)
       enddo
       n1=n1+nn
      enddo
      n1=1
      do k=start_y,end_y
       print*,' reading year ',k
       call sub_flush(6)
       write(file,'("data/ncep_daily/uflx.sfc.gauss.",i4,".nc")') k
       iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=nf_inq_varid(ncid,'uflx',varid)
       iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
       iret=NF_INQ_DIMLEN (NCID, DIMs(3), nn)
       corner=(/1,1,1,1/); edges = (/nx,ny,nn,1/)
       allocate( v4(nx,ny,nn) )
       iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
       add_offset=0.0
       iret= nf_get_att_real (ncid,varid ,'add_offset',add_offset)
       scale_factor=1.0
       iret= nf_get_att_real (ncid,varid ,'scale_factor',scale_factor)
c    convert N/m**2 to dyn/cm**2 and change sign
       do n=n1,n1+nn-1; do j=1,ny; do i=1,nx
          if (mask(i,ny-j+1) == 0.0 ) then
           tau(i,ny-j+1,n,1) = 
     &       (v4(i,j,n-n1+1)*scale_factor+add_offset)*(-10.0)
          else
           tau(i,ny-j+1,n,1) = spval
          endif
       enddo; enddo; enddo
       call ncclos (ncid, iret)
       write(file,'("data/ncep_daily/vflx.sfc.gauss.",i4,".nc")') k
       iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=nf_inq_varid(ncid,'vflx',varid)
       iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
       add_offset=0.0
       iret= nf_get_att_real (ncid,varid ,'add_offset',add_offset)
       scale_factor=1.0
       iret= nf_get_att_real (ncid,varid ,'scale_factor',scale_factor)
c    convert N/m**2 to dyn/cm**2 and change sign
       do n=n1,n1+nn-1; do j=1,ny; do i=1,nx
          if (mask(i,ny-j+1) == 0.0 ) then
           tau(i,ny-j+1,n,2) = 
     &        (v4(i,j,n-n1+1)*scale_factor+add_offset)*(-10.0)
          else
           tau(i,ny-j+1,n,2) = spval
          endif
       enddo; enddo; enddo
       call ncclos (ncid, iret)
       deallocate(v4)
       n1=n1+nn
      enddo
      end subroutine read_tau


#elif defined read_tau_ncep_daily_large_file
c
c--------------------------------------------------------
c     use daily wind stress data from the
c     NCEP/NCAR reanalysis project
c     start_y : year of forcing
c     fraction_y : year is divided into n fractions
c     nfraction_y : actual fraction within year
c--------------------------------------------------------

      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4),nn,k
      character (len=80) :: file
      integer :: start_y, fraction_y, nfraction_y,full_y
      common / tau_ncep/ start_y,fraction_y,nfraction_y,full_y
      namelist /ncep_time/ start_y,fraction_y,nfraction_y

      call getunit(ncid,'namelist.ncep_time','fsr')
      read(ncid, nml = ncep_time) 
      close(ncid)

      print*,' reading daily wind stress data from NCEP '
      print*,' year ',start_y,' in ',fraction_y,' fractions'
      print*,' current fraction ',nfraction_y
      write(file,'("data/ncep_daily/uflx.sfc.gauss.",i4,".nc")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'uflx',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nn)
      full_y = nn
      print*,' days in year:',full_y
      if (mod(full_y,fraction_y)/=0) then
        print*,'fraction must be integer divider of 365'
        print*,'but ', mod(full_y,fraction_y)
        stop
      endif
      nt=full_y/fraction_y ! must be an integer
      call ncclos (ncid, iret)
      end subroutine read_dims_tau


      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      use time_type_module
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt)
#include "netcdf.inc"
      integer :: iret,ncid,varid,corner(4),edges(4),k,dims(4),nn
      integer :: n1,i,j,n,xid,yid
      character (len=80) :: file
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:,:,:)
      real (kind=4) :: add_offset,scale_factor
#else
      real , allocatable :: v4(:,:,:)
      real  :: add_offset,scale_factor
#endif
      type(time_type) :: t1,t2,t3
      integer :: start_y, fraction_y, nfraction_y,full_y
      common / tau_ncep/ start_y,fraction_y,nfraction_y,full_y
c   read land mask
      iret=nf_open('data/ncep_daily/lsmask.19294.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'lsmask',varid)
      corner=(/1,1,1,1/); edges = (/nx,ny,1,1/)
      allocate( v4(nx,ny,1) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      do j=1,ny; do i=1,nx
         mask(i,j)=v4(i,ny-j+1,1)
      enddo; enddo
      deallocate(v4)
      ! read grid info
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),file) 
      iret=nf_inq_varid(ncid,file,xid)
      allocate(v4(nx,1,1) )
      iret= nf_get_vara_real (ncid,xid ,1,nx,v4); x=v4(1:nx,1,1)
      deallocate(v4)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),file) 
      iret=nf_inq_varid(ncid,file,yid)
      allocate(v4(ny,1,1) )
      iret= nf_get_vara_real (ncid,yid ,1,ny,v4); 
      do j=1,ny
         y(j)=v4(ny-j+1,1,1)
      enddo
      deallocate(v4)
      call ncclos (ncid, iret)
      ! time vector
      t1 =  set_date(start_y, 1, 1, 0, 0, 0)
      toffset=t1%days + (nfraction_y-1)*(full_y/fraction_y)
      write(file,'("data/ncep_daily/uflx.sfc.gauss.",i4,".nc")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'uflx',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      call ncclos (ncid, iret)
      do n=1,nt
       time(n) =  365.0/(full_y*1.0)
      enddo


      print*,' reading year ',start_y
      call sub_flush(6)
      write(file,'("data/ncep_daily/uflx.sfc.gauss.",i4,".nc")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'uflx',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nn)
      n1=1+(nfraction_y-1)*(full_y/fraction_y)
      nn=nfraction_y*(full_y/fraction_y)
     
      corner=(/1,1,n1,1/); edges = (/nx,ny,nt,1/)
      allocate( v4(nx,ny,nt) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      add_offset=0.0
      iret= nf_get_att_real (ncid,varid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,varid ,'scale_factor',scale_factor)
c    convert N/m**2 to dyn/cm**2 and change sign
      do n=1,nt; do j=1,ny; do i=1,nx
          if (mask(i,ny-j+1) == 0.0 ) then
           tau(i,ny-j+1,n,1) = 
     &       (v4(i,j,n)*scale_factor+add_offset)*(-10.0)
          else
           tau(i,ny-j+1,n,1) = spval
          endif
      enddo; enddo; enddo
      call ncclos (ncid, iret)

      write(file,'("data/ncep_daily/vflx.sfc.gauss.",i4,".nc")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'vflx',varid)
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      add_offset=0.0
      iret= nf_get_att_real (ncid,varid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,varid ,'scale_factor',scale_factor)
c    convert N/m**2 to dyn/cm**2 and change sign
      do n=1,nt; do j=1,ny; do i=1,nx
          if (mask(i,ny-j+1) == 0.0 ) then
           tau(i,ny-j+1,n,2) = 
     &        (v4(i,j,n)*scale_factor+add_offset)*(-10.0)
          else
           tau(i,ny-j+1,n,2) = spval
          endif
      enddo; enddo; enddo
      call ncclos (ncid, iret)
      deallocate(v4)

      end subroutine read_tau

#elif defined read_tau_T511_daily_large_file
c
c--------------------------------------------------------
c     use daily wind stress data from the
c     ECMWF model with T511 resolution
c     three parameter have to be set: start_y, fraction_y and nfraction_y
c     start_y : year of forcing
c     fraction_y : year is divided into n fractions
c     nfraction_y : actual fraction within year
c--------------------------------------------------------

      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4),nn,k
      character (len=80) :: file
      integer :: start_y, fraction_y, nfraction_y,full_y
      common / tau_ncep/ start_y,fraction_y,nfraction_y,full_y
      namelist /T511_time/ start_y,fraction_y,nfraction_y

      call getunit(ncid,'namelist.T511_time','fsr')
      read(ncid, nml = T511_time) 
      close(ncid)
      print*,' reading daily wind stress data from ECMWF/T511'
      print*,' year ',start_y,' in ',fraction_y,' fractions'
      print*,' current fraction ',nfraction_y
      write(file,'("data/ecmwf_T511/EWSS_",i4,".cdf")') start_y
      print*,'reading from file ',file
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'EWSS',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nn)
      full_y = nn
      print*,' days in year:',full_y
      if (mod(full_y,fraction_y)/=0) then
        print*,'fraction must be integer divider of 365'
        print*,'but ', mod(full_y,fraction_y)
        stop
      endif
      nt=full_y/fraction_y ! must be an integer
      call ncclos (ncid, iret)
      end subroutine read_dims_tau


      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      use time_type_module
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt)
#include "netcdf.inc"
      integer :: iret,ncid,varid,corner(4),edges(4),k,dims(4),nn
      integer :: n1,i,j,n,xid,yid
      character (len=80) :: file
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:,:,:)
      real (kind=4) :: add_offset,scale_factor
#else
      real , allocatable :: v4(:,:,:)
      real  :: add_offset,scale_factor
#endif
      real :: mask2(nx,ny),fxa
      type(time_type) :: t1,t2,t3
      integer :: start_y, fraction_y, nfraction_y,full_y
      common / tau_ncep/ start_y,fraction_y,nfraction_y,full_y
c   read land mask
      iret=nf_open('data/ecmwf_T511/land-sea-mask.nc.new',
     &              NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'tcw',varid)
      corner=(/1,1,1,1/); edges = (/nx,ny,1,1/)
      allocate( v4(nx,ny,1) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      add_offset=0.0
      iret= nf_get_att_real (ncid,varid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,varid ,'scale_factor',scale_factor)
      do j=1,ny; do i=1,nx
        mask(i,j)=0.0
        fxa = v4(i,ny-j+1,1)*scale_factor + add_offset
        if (fxa.gt. 0.1)  mask(i,j)=1.
        mask2(i,j)=0.0
        if (fxa.gt. 0.5)  mask2(i,j)=1.
      enddo; enddo
      deallocate(v4)
      call ncclos (ncid, iret)

c     we have to manipulate the land mask a bit
      mask(705:712,288:300)=mask2(705:712,288:300) ! this is Scotland
      mask(635:645,306:312)=mask2(635:645,306:312) ! East greenland
      mask(585:600,310:318)=mask2(585:600,310:318)

      ! read grid info from other file (because of corrupt netcdf-lib in Kiel)
      ! this lib cannot read double variables into floats
      iret=nf_open('data/ecmwf_T511/land-sea-mask.nc',NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'lsm',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),file) 
      iret=nf_inq_varid(ncid,file,xid)
      allocate(v4(nx,1,1) )
      print*,' reading axis ',file
      iret= nf_get_vara_real(ncid,xid,1,nx,v4); 
      x=v4(1:nx,1,1)
      deallocate(v4)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),file) 
      iret=nf_inq_varid(ncid,file,yid)
      allocate(v4(ny,1,1) )
      print*,' reading axis ',file
      iret= nf_get_vara_real (ncid,yid ,1,ny,v4); 
      do j=1,ny
         y(j)=v4(ny-j+1,1,1)
      enddo
      deallocate(v4)
      call ncclos (ncid, iret)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      ! time vector
      t1 =  set_date(start_y, 1, 1, 0, 0, 0)
      toffset=t1%days + (nfraction_y-1)*(full_y/fraction_y)
      write(file,'("data/ecmwf_T511/EWSS_",i4,".cdf")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'EWSS',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      call ncclos (ncid, iret)
      do n=1,nt
       time(n) =  365.0/(full_y*1.0)
      enddo

      print*,' reading year ',start_y
      call sub_flush(6)
      write(file,'("data/ecmwf_T511/EWSS_",i4,".cdf")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'EWSS',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nn)
      n1=1+(nfraction_y-1)*(full_y/fraction_y)
      nn=nfraction_y*(full_y/fraction_y)
     
      corner=(/1,1,n1,1/); edges = (/nx,ny,nt,1/)
      allocate( v4(nx,ny,nt) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
c    convert N/m**2 to dyn/cm**2 
      do n=1,nt; do j=1,ny; do i=1,nx
          if (mask(i,j) /= 1.0 ) then
           tau(i,j,n,1) =v4(i,j,n)*10.0
          else
           tau(i,j,n,1) = spval
          endif
      enddo; enddo; enddo
      call ncclos (ncid, iret)

      write(file,'("data/ecmwf_T511/NSSS_",i4,".cdf")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'NSSS',varid)
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
c    convert N/m**2 to dyn/cm**2 
      do n=1,nt; do j=1,ny; do i=1,nx
          if (mask(i,j) /= 1.0 ) then
           tau(i,j,n,2) =v4(i,j,n)*10.0
          else
           tau(i,j,n,2) = spval
          endif
      enddo; enddo; enddo
      call ncclos (ncid, iret)
      deallocate(v4)

      end subroutine read_tau

#elif defined read_tau_T511_daily_large_file_smooth
c
c--------------------------------------------------------
c     use smoothed daily wind stress data from the
c     ECMWF model with T511 resolution
c     three parameter have to be set: start_y, fraction_y and nfraction_y
c     start_y : year of forcing
c     fraction_y : year is divided into n fractions
c     nfraction_y : actual fraction within year
c--------------------------------------------------------

      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4),nn,k
      character (len=80) :: file
      integer :: start_y, fraction_y, nfraction_y,full_y
      common / tau_ncep/ start_y,fraction_y,nfraction_y,full_y
      namelist /T511_time/ start_y,fraction_y,nfraction_y

      call getunit(ncid,'namelist.T511_time','fsr')
      read(ncid, nml = T511_time) 
      close(ncid)
      print*,' reading daily wind stress data from ECMWF/T511'
      print*,' year ',start_y,' in ',fraction_y,' fractions'
      print*,' current fraction ',nfraction_y
      write(file,'("data/ecmwf_T511/smooth/EWSS_",i4,".cdf")') start_y
      print*,'reading from file ',file
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'EWSS',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nn)
      full_y = nn
      print*,' days in year:',full_y
      if (mod(full_y,fraction_y)/=0) then
        print*,'fraction must be integer divider of 365'
        print*,'but ', mod(full_y,fraction_y)
        stop
      endif
      nt=full_y/fraction_y ! must be an integer
      call ncclos (ncid, iret)
      end subroutine read_dims_tau


      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      use time_type_module
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt)
#include "netcdf.inc"
      integer :: iret,ncid,varid,corner(4),edges(4),k,dims(4),nn
      integer :: n1,i,j,n,xid,yid
      character (len=80) :: file
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:,:,:)
      real (kind=4) :: add_offset,scale_factor,sp
#else
      real , allocatable :: v4(:,:,:)
      real  :: add_offset,scale_factor,sp
#endif
      integer :: indp
      real :: mask2(nx,ny),fxa
      type(time_type) :: t1,t2,t3
      integer :: start_y, fraction_y, nfraction_y,full_y
      common / tau_ncep/ start_y,fraction_y,nfraction_y,full_y
c   read land mask
      iret=nf_open('data/ecmwf_T511/land-sea-mask.nc.new',
     &              NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'tcw',varid)
      corner=(/1,1,1,1/); edges = (/nx,ny,1,1/)
      allocate( v4(nx,ny,1) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      add_offset=0.0
      iret= nf_get_att_real (ncid,varid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,varid ,'scale_factor',scale_factor)
      do j=1,ny; do i=1,nx
        mask(i,j)=0.0
        fxa = v4(i,ny-j+1,1)*scale_factor + add_offset
        if (fxa.gt. 0.1)  mask(i,j)=1.
        mask2(i,j)=0.0
        if (fxa.gt. 0.5)  mask2(i,j)=1.
      enddo; enddo
      deallocate(v4)
      call ncclos (ncid, iret)

c     we have to manipulate the land mask a bit
      mask(705:712,288:300)=mask2(705:712,288:300) ! this is Scotland
      mask(635:645,306:312)=mask2(635:645,306:312) ! East greenland
      mask(585:600,310:318)=mask2(585:600,310:318)

      ! read grid info from other file (because of corrupt netcdf-lib in Kiel)
      ! this lib cannot read double variables into floats
      iret=nf_open('data/ecmwf_T511/land-sea-mask.nc',NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'lsm',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),file) 
      iret=nf_inq_varid(ncid,file,xid)
      allocate(v4(nx,1,1) )
      print*,' reading axis ',file
      iret= nf_get_vara_real(ncid,xid,1,nx,v4); 
      x=v4(1:nx,1,1)
      deallocate(v4)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),file) 
      iret=nf_inq_varid(ncid,file,yid)
      allocate(v4(ny,1,1) )
      print*,' reading axis ',file
      iret= nf_get_vara_real (ncid,yid ,1,ny,v4); 
      do j=1,ny
         y(j)=v4(ny-j+1,1,1)
      enddo
      deallocate(v4)
      call ncclos (ncid, iret)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      ! time vector
      t1 =  set_date(start_y, 1, 1, 0, 0, 0)
      toffset=t1%days + (nfraction_y-1)*(full_y/fraction_y)
      write(file,'("data/ecmwf_T511/smooth/EWSS_",i4,".cdf")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'EWSS',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      call ncclos (ncid, iret)
      do n=1,nt
       time(n) =  365.0/(full_y*1.0)
      enddo

      print*,' reading year ',start_y
      call sub_flush(6)
      write(file,'("data/ecmwf_T511/smooth/EWSS_",i4,".cdf")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'EWSS',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nn)
      n1=1+(nfraction_y-1)*(full_y/fraction_y)
      nn=nfraction_y*(full_y/fraction_y)
     
      corner=(/1,1,n1,1/); edges = (/nx,ny,nt,1/)
      allocate( v4(nx,ny,nt) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      iret= nf_get_att_real  (ncid,varid ,'missing_value',sp)
c    convert N/m**2 to dyn/cm**2 
      do n=1,nt; do j=1,ny; do i=1,nx
          if (v4(i,j,n) /= sp ) then
           tau(i,j,n,1) =v4(i,j,n)*10.0
           mask(i,j)=0.
          else
           tau(i,j,n,1) = spval
           mask(i,j)=1.
          endif
      enddo; enddo; enddo
      call ncclos (ncid, iret)

      write(file,'("data/ecmwf_T511/smooth/NSSS_",i4,".cdf")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'NSSS',varid)
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      iret= nf_get_att_real  (ncid,varid ,'missing_value',sp)
c    convert N/m**2 to dyn/cm**2 
      do n=1,nt; do j=1,ny; do i=1,nx
          if (v4(i,j,n) /= sp ) then
           tau(i,j,n,2) =v4(i,j,n)*10.0
          else
           tau(i,j,n,2) = spval
          endif
      enddo; enddo; enddo
      call ncclos (ncid, iret)
      deallocate(v4)

      i=indp(290.2,x,nx); j=indp(11.51,y,ny)
      tau(i,j,:,:) = 0.; mask(i,j)=0.
      i=indp(276.5,x,nx); j=indp(12.73,y,ny)
      tau(i,j,:,:) = 0.; mask(i,j)=0.
      i=indp(354.3,x,nx); j=indp(56.38,y,ny)
      tau(i,j,:,:) = 0.; mask(i,j)=0.
      i=indp(354.1,x,nx); j=indp(57.43,y,ny)
      tau(i,j,:,:) = 0.; mask(i,j)=0.
      i=indp(319.1,x,nx); j=indp(64.83,y,ny)
      tau(i,j,:,:) = 0.; mask(i,j)=0.
      i=indp(296.6,x,nx); j=indp(65.39,y,ny)
      tau(i,j,:,:) = 0.; mask(i,j)=0.

      end subroutine read_tau


#elif defined read_tau_ncep_clim
c
c--------------------------------------------------------
c     use climatological monthly mean wind stress data from the
c     NCEP/NCAR reanalysis project
c--------------------------------------------------------

      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4)
      character (len=80) :: file
      print*,' reading clim. monthly wind stress data from NCEP '
      file='data/ncep_for_pml/uflx.mon.ltm.nc'
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'uflx',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt)
      call ncclos (ncid, iret)
      end subroutine read_dims_tau


      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      use time_type_module
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt)
#include "netcdf.inc"
      integer :: iret,ncid,varid,corner(4),edges(4),k,dims(4),nn
      integer :: n1,i,j,n,xid,yid
      character (len=80) :: file
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:,:,:)
      real (kind=4) :: add_offset,scale_factor
#else
      real , allocatable :: v4(:,:,:)
      real  :: add_offset,scale_factor
#endif

c   read land mask
      iret=nf_open('data/ncep_daily/lsmask.19294.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'lsmask',varid)
      corner=(/1,1,1,1/); edges = (/nx,ny,1,1/)
      allocate( v4(nx,ny,1) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      do j=1,ny; do i=1,nx
         mask(i,j)=v4(i,ny-j+1,1)
      enddo; enddo
      deallocate(v4)
      ! read grid info
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),file) 
      iret=nf_inq_varid(ncid,file,xid)
      allocate(v4(nx,1,1) )
      iret= nf_get_vara_real (ncid,xid ,1,nx,v4); x=v4(1:nx,1,1)
      deallocate(v4)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),file) 
      iret=nf_inq_varid(ncid,file,yid)
      allocate(v4(ny,1,1) )
      iret= nf_get_vara_real (ncid,yid ,1,ny,v4); 
      do j=1,ny
         y(j)=v4(ny-j+1,1,1)
      enddo
      deallocate(v4)
      call ncclos (ncid, iret)

      file='data/ncep_for_pml/uflx.mon.ltm.nc'
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'uflx',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      corner=(/1,1,1,1/); edges = (/nx,ny,nt,1/)
      allocate( v4(nx,ny,nt) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      add_offset=0.0
      iret= nf_get_att_real (ncid,varid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,varid ,'scale_factor',scale_factor)
c    convert N/m**2 to dyn/cm**2 and change sign
      do j=1,ny
       do i=1,nx
          if (mask(i,ny-j+1) == 0.0 ) then
           tau(i,ny-j+1,:,1) = 
     &       (v4(i,j,:)*scale_factor+add_offset)*(-10.0)
          else
           tau(i,ny-j+1,:,1) = spval
          endif
       enddo
      enddo
      call ncclos (ncid, iret)
      deallocate(v4)


      file='data/ncep_for_pml/vflx.mon.ltm.nc'
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'vflx',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      corner=(/1,1,1,1/); edges = (/nx,ny,nt,1/)
      allocate( v4(nx,ny,nt) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      add_offset=0.0
      iret= nf_get_att_real (ncid,varid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,varid ,'scale_factor',scale_factor)
c    convert N/m**2 to dyn/cm**2 and change sign
      do j=1,ny
       do i=1,nx
          if (mask(i,ny-j+1) == 0.0 ) then
           tau(i,ny-j+1,:,2) = 
     &       (v4(i,j,:)*scale_factor+add_offset)*(-10.0)
          else
           tau(i,ny-j+1,:,2) = spval
          endif
       enddo
      enddo
      call ncclos (ncid, iret)
      deallocate(v4)
c   time
      toffset=0.0
      time = (/31,28,31,30,31,30,31,31,30,31,30,31/)
      end subroutine read_tau

#elif defined read_tau_ncep
c
c--------------------------------------------------------
c     use monthly mean wind stress data from the
c     NCEP/NCAR reanalysis project
c--------------------------------------------------------

      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4)
      character (len=80) :: file
      integer :: start_y, nyear, year0
      common / flx_ncep/ start_y,nyear, year0

      start_y=1999
      nyear = 5
      year0=1948

      print*,' reading clim. monthly wind stress data from NCEP '
      print*,' start year = ',start_y
      print*,' number of years = ',nyear

      file='data/ncep/uflx.mon.mean.nc'
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'uflx',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt)
      nt=12*nyear
      call ncclos (ncid, iret)
      end subroutine read_dims_tau


      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      use time_type_module
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt)
#include "netcdf.inc"
      integer :: iret,ncid,varid,corner(4),edges(4),k,dims(4),nn
      integer :: n1,i,j,n,xid,yid
      character (len=80) :: file
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:,:,:)
      real (kind=4) :: add_offset,scale_factor
#else
      real , allocatable :: v4(:,:,:)
      real  :: add_offset,scale_factor
#endif
      integer :: start_y, nyear, year0
      common / flx_ncep/ start_y,nyear, year0
      type(time_type) :: t1,t2,t3
      real :: dpm(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)

c   read land mask
      iret=nf_open('data/ncep/lsmask.19294.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'lsmask',varid)
      corner=(/1,1,1,1/); edges = (/nx,ny,1,1/)
      allocate( v4(nx,ny,1) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      do j=1,ny; do i=1,nx
         mask(i,j)=v4(i,ny-j+1,1)
      enddo; enddo
      deallocate(v4)
      ! read grid info
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),file) 
      iret=nf_inq_varid(ncid,file,xid)
      allocate(v4(nx,1,1) )
      iret= nf_get_vara_real (ncid,xid ,1,nx,v4); x=v4(1:nx,1,1)
      deallocate(v4)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),file) 
      iret=nf_inq_varid(ncid,file,yid)
      allocate(v4(ny,1,1) )
      iret= nf_get_vara_real (ncid,yid ,1,ny,v4); 
      do j=1,ny
         y(j)=v4(ny-j+1,1,1)
      enddo
      deallocate(v4)
      call ncclos (ncid, iret)

      file='data/ncep/uflx.mon.mean.nc'
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'uflx',varid)
      n=(start_y-year0)*12
      corner=(/1,1,n+1,1/); edges = (/nx,ny,nt,1/)
      allocate( v4(nx,ny,nt) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      add_offset=0.0
      iret= nf_get_att_real (ncid,varid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,varid ,'scale_factor',scale_factor)
c    convert N/m**2 to dyn/cm**2 and change sign
      do j=1,ny
       do i=1,nx
          if (mask(i,ny-j+1) == 0.0 ) then
           tau(i,ny-j+1,:,1) = 
     &       (v4(i,j,:)*scale_factor+add_offset)*(-10.0)
          else
           tau(i,ny-j+1,:,1) = spval
          endif
       enddo
      enddo
      call ncclos (ncid, iret)
      deallocate(v4)


      file='data/ncep/vflx.mon.mean.nc'
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'vflx',varid)
      n=(start_y-year0)*12
      corner=(/1,1,n+1,1/); edges = (/nx,ny,nt,1/)
      allocate( v4(nx,ny,nt) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      add_offset=0.0
      iret= nf_get_att_real (ncid,varid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,varid ,'scale_factor',scale_factor)
c    convert N/m**2 to dyn/cm**2 and change sign
      do j=1,ny
       do i=1,nx
          if (mask(i,ny-j+1) == 0.0 ) then
           tau(i,ny-j+1,:,2) = 
     &       (v4(i,j,:)*scale_factor+add_offset)*(-10.0)
          else
           tau(i,ny-j+1,:,2) = spval
          endif
       enddo
      enddo
      call ncclos (ncid, iret)
      deallocate(v4)
c   time
      t1 =  set_date(start_y, 1, 1, 0, 0, 0)
      toffset=t1%days 
      do n=1,nyear
       do j=1,12
        time((n-1)*12+j )=dpm(j)
       enddo
      enddo
      end subroutine read_tau

#elif defined read_tau_omip
c
c--------------------------------------------------------
c    use a monthly mean climatology from the OMIP II
c    data (derived from ECMWF reanalysis data)
c--------------------------------------------------------
c
      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4)
      print*,' reading monthly wind stress data from OMIP'
      iret=nf_open('data/OMIP/mean_zonal_stress.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'stress_x',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt); nt=12
      call ncclos (ncid, iret)
      end subroutine read_dims_tau

      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt)
#include "netcdf.inc"
      integer :: iret,ncid,txid,tyid,dims(4),aid,n
      character (len=80) :: name
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:,:,:)
#else
      real, allocatable :: v4(:,:,:)
#endif
      allocate(v4(nx,ny,nt))
c   read taux
      iret=nf_open('data/OMIP/mean_zonal_stress.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'stress_x',txid)
      iret= nf_get_vara_real(ncid,txid,(/1,1,1/),(/nx,ny,nt/),v4)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      tau(:,:,:,1)=v4
      ! read grid info
      iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
      iret=nf_inq_varid(ncid,name,aid)
      iret= nf_get_vara_real(ncid,aid ,1,nx,v4(:,1,1))
      x=v4(:,1,1)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
      iret=nf_inq_varid(ncid,name,aid)
      iret= nf_get_vara_real (ncid,aid ,1,ny,v4(1,:,1)) 
      y=v4(1,:,1)
      call ncclos (ncid, iret)
c   read tauy
      iret=nf_open('data/OMIP/mean_merid_stress.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'stress_y',txid)
      iret= nf_get_vara_real(ncid,txid,(/1,1,1/),(/nx,ny,nt/),v4)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      tau(:,:,:,2)= v4
      call ncclos (ncid, iret)
c   read land mask
      iret=nf_open('data/OMIP/land_sea_mask.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'lsm',txid)
      iret= nf_get_vara_real(ncid,txid,(/1,1/),(/nx,ny/),v4(1,1,1))
      mask=v4(:,:,1)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      call ncclos (ncid, iret)
c   apply land mask
      do n=1,nt
        where ( mask == 0.0 ) tau(:,:,n,1)=spval
        where ( mask == 0.0 ) tau(:,:,n,2)=spval
      enddo
      deallocate(v4)
c   flip second dimension up/down
      y=y(ny:1:-1)
      tau=tau(:,ny:1:-1,:,:)
      mask=mask(:,ny:1:-1)
c   convert N/m^2 to dyn/cm^2
      do n=1,nt
        where ( mask /= 0.0 ) tau(:,:,n,1)=tau(:,:,n,1)*10.0
        where ( mask /= 0.0 ) tau(:,:,n,2)=tau(:,:,n,2)*10.0
      enddo
c   time
      toffset=0.0
      time = (/31,28,31,30,31,30,31,31,30,31,30,31/)
      end subroutine read_tau
c
#elif defined read_tau_dasilva
c
c--------------------------------------------------------
c    use a monthly mean climatology from Da Silva 1994
c    (derived from COADS data)
c--------------------------------------------------------
c
      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4)
      print*,' reading monthly wind stress data from DaSilva94'
      iret=nf_open('data/DaSilva94/taux.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'taux3',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt); nt=12
      call ncclos (ncid, iret)
      end subroutine read_dims_tau

      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt)
#include "netcdf.inc"
      integer :: iret,ncid,txid,tyid,dims(4),aid,n
      character (len=80) :: name
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:,:,:)
      real (kind=4) :: sp,add_offset,scale_factor
#else
      real, allocatable :: v4(:,:,:)
      real  :: sp,add_offset,scale_factor
#endif
      allocate(v4(nx,ny,nt))
c   read taux
      iret=nf_open('data/DaSilva94/taux.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'taux3',txid)
      iret= nf_get_vara_real(ncid,txid,(/1,1,1/),(/nx,ny,nt/),v4)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      tau(:,:,:,1)=v4
      iret= nf_get_att_real (ncid,txid ,'missing_value',sp)
      where (tau(:,:,:,1) == sp ) tau(:,:,:,1)=spval
      add_offset=0.0
      iret= nf_get_att_real (ncid,txid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,txid ,'scale_factor',scale_factor)
      where (tau(:,:,:,1)/=spval) 
     &     tau(:,:,:,1)=tau(:,:,:,1)*scale_factor+add_offset
      ! read grid info
      iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
      iret=nf_inq_varid(ncid,name,aid)
      iret= nf_get_vara_real(ncid,aid ,1,nx,v4(:,1,1))
      x=v4(:,1,1)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
      iret=nf_inq_varid(ncid,name,aid)
      iret= nf_get_vara_real (ncid,aid ,1,ny,v4(1,:,1)) 
      y=v4(1,:,1)
      call ncclos (ncid, iret)
c   read tauy
      iret=nf_open('data/DaSilva94/tauy.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'tauy3',txid)
      iret= nf_get_vara_real(ncid,txid,(/1,1,1/),(/nx,ny,nt/),v4)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      tau(:,:,:,2)=v4
      iret= nf_get_att_real (ncid,txid ,'missing_value',sp)
      where (tau(:,:,:,2) == sp ) tau(:,:,:,2)=spval
      add_offset=0.0
      iret= nf_get_att_real (ncid,txid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,txid ,'scale_factor',scale_factor)
      where (tau(:,:,:,2)/=spval) 
     &     tau(:,:,:,2)=tau(:,:,:,2)*scale_factor+add_offset
      call ncclos (ncid, iret)
      deallocate(v4)
c   read land mask
      mask=1.0
      where ( tau(:,:,1,1)==spval)mask=0.0
c   convert N/m^2 to dyn/cm^2
      do n=1,nt
        where ( mask /= 0.0 ) tau(:,:,n,1)=tau(:,:,n,1)*10.0
        where ( mask /= 0.0 ) tau(:,:,n,2)=tau(:,:,n,2)*10.0
      enddo
c   time
      toffset=0.0
      time = (/31,28,31,30,31,30,31,31,30,31,30,31/)
      end subroutine read_tau

#elif defined read_tau_ecco
c
c--------------------------------------------------------
c    use a monthly mean climatology from ECCO
c    (assimilated data into 2x2 deg ocean model. 
c     Control parameter are initial and surface boundary 
c     conditions. Assimilation period is 1992-1997.
c     Climatology is taken as average over this period.)
c--------------------------------------------------------
c
      subroutine read_dims_tau(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4)
      print*,' reading monthly wind stress data from ECCO'
      iret=nf_open('data/ECCO/tau_mean.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'taux',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), nx) ! axis are swaped 
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), ny) ! somehow
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt)  ! ordering is T,Y,Z
      call ncclos (ncid, iret)
      end subroutine read_dims_tau

      subroutine read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
      implicit none
      integer :: nx,ny,nt
      real :: spval, tau(nx,ny,nt,2) ,x(nx), y(ny),mask(nx,ny)
      real :: toffset,time(nt)
#include "netcdf.inc"
      integer :: iret,ncid,txid,tyid,dims(4),aid,n,j
      character (len=80) :: name
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:,:,:)
      real (kind=4) :: sp,add_offset,scale_factor
#else
      real, allocatable :: v4(:,:,:)
      real  :: sp,add_offset,scale_factor
#endif
      allocate(v4(ny,nx,nt)) ! axes are somehow swaped in file
c   read taux
      iret=nf_open('data/ECCO/tau_mean.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'taux',txid)
      iret= nf_get_vara_real(ncid,txid,(/1,1,1/),(/ny,nx,nt/),v4)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      do n=1,nt
      do j=1,ny
        tau(:,j,n,1)=v4(j,:,n)
      enddo
      enddo
      iret= nf_get_att_real (ncid,txid ,'missing_value',sp)
      where (tau(:,:,:,1) == sp ) tau(:,:,:,1)=spval
      ! read grid info
      iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
      iret=nf_inq_varid(ncid,name,aid)
      iret= nf_get_vara_real(ncid,aid ,1,ny,v4(:,1,1))
      y=v4(:,1,1)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
      iret=nf_inq_varid(ncid,name,aid)
      iret= nf_get_vara_real (ncid,aid ,1,nx,v4(1,:,1)) 
      x=v4(1,:,1)
      call ncclos (ncid, iret)
c   read tauy
      iret=nf_open('data/ECCO/tau_mean.cdf',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'tauy',txid)
      iret= nf_get_vara_real(ncid,txid,(/1,1,1/),(/ny,nx,nt/),v4)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      do n=1,nt
      do j=1,ny
        tau(:,j,n,2)=v4(j,:,n)
      enddo
      enddo
      iret= nf_get_att_real (ncid,txid ,'missing_value',sp)
      where (tau(:,:,:,2) == sp ) tau(:,:,:,2)=spval
      call ncclos (ncid, iret)
      deallocate(v4)
c   read land mask
      mask=1.0
      where ( tau(:,:,1,1)==spval) mask=0.0
c   convert N/m^2 to dyn/cm^2
      do n=1,nt
        where ( mask /= 0.0 ) tau(:,:,n,1)=tau(:,:,n,1)*10.0
        where ( mask /= 0.0 ) tau(:,:,n,2)=tau(:,:,n,2)*10.0
      enddo
c   time
      toffset=0.0
      time = (/31,28,31,30,31,30,31,31,30,31,30,31/)
      end subroutine read_tau

#endif





#ifdef read_ustar_ecmwf
c
c--------------------------------------------------------
c    use a monthly mean climatology from the ECMWF
c    DYNAMO data for friction velocity.
c--------------------------------------------------------
c
       subroutine read_dims_ustar(nx,ny,nt)
       implicit none
       integer nx,ny,nt,io
       character*32 stamp
       real (kind=4) :: dpm
       integer (kind=4) nx_,ny_
       print*,' reading frcition velocity data from ECMWF'
       open(10,file='data/ecmwf/ecmwf_BB.ustar',form='unformatted')
       read (10); read (10) stamp, dpm, nx_, ny_
       nx=nx_;ny=ny_; nt=12;
       close(10)
       end subroutine read_dims_ustar

       subroutine read_ustar(ustar,x,y,nx,ny,nt,spval,time,toffset)
       implicit none
       integer :: nx,ny,nt
       real :: spval, ustar(nx,ny,nt) ,x(nx), y(ny),time(nt),toffset
       character*32 s
       real (kind=4) u4(nx,ny) 
       real (kind=4) x4(nx), y4(ny)
       real (kind=4) dpm,   mask(nx,ny)
       integer (kind=4) k,n
       real :: factor,rho_air = 1.29, rho_wat = 1026.
       open(10,file='data/ecmwf/ecmwf_BB.ustar',form='unformatted')
       do n=1,nt
        read (10); read (10) s, dpm, k, k, k, x4,y4,u4(:,:)
        x=x4;y=y4; ustar(:,:,n)=u4
       enddo
       close(10)
c    convert ustar
c    the ustar**3 field is referenced to air-density (1.29 kg/m**3), 
c    therefore re-reference to density of water (1035 kg/m**3)
c    i.e. rho_wat*(vstar**2)=rho_air*(ustar**2), when vstar is the 
c    friction velocity in water and ustar the one in air.
!    reference to dens. in water and convert (m/s)**3 to (cm/s)**3
       factor = 1.e+06 * (rho_air/rho_wat)**1.5
       ustar = factor*ustar
!   apply land mask
       open(10,file='data/ecmwf/ecmwf_BB.mask',form='unformatted')
       read (10) ; read (10); read (10) mask; close(10)
       do n=1,nt
        where (mask==0.) ustar(:,:,n) = spval
       enddo
       toffset=0.0
       time = (/31,28,31,30,31,30,31,31,30,31,30,31/)
       end subroutine read_ustar

#elif defined read_ustar_ncep_monthly
c
c--------------------------------------------------------
c    use a monthly mean climatology from the ECMWF
c    DYNAMO data for friction velocity.
c--------------------------------------------------------
c
       subroutine read_dims_ustar(nx,ny,nt)
       implicit none
       integer nx,ny,nt,io
       character*32 stamp
       real (kind=4) :: dpm
       integer (kind=4) nx_,ny_
#include "netcdf.inc"
       integer :: iret, ncid,dims(4),vid

       print*,' reading frcition velocity data from ECMWF '
       print*,' plus anomalies from NCEP/NCAR'
       open(10,file='data/ecmwf/ecmwf_BB.ustar',form='unformatted')
       read (10); read (10) stamp, dpm, nx_, ny_
       nx=nx_;ny=ny_; 
       close(10)
       iret=nf_open('data/ncep/ustar_anom.cdf',NF_nowrite,ncid)
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=nf_inq_varid(ncid,'ustar',vid)
       iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
       iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt)
       call ncclos (ncid, iret)
       end subroutine read_dims_ustar



       subroutine read_ustar(ustar,x,y,nx,ny,nt,spval,time,toffset)
       use time_type_module
       implicit none
       integer :: nx,ny,nt
       real :: spval, ustar(nx,ny,nt) ,x(nx), y(ny),time(nt),toffset
       character*32 s
       real (kind=4) u4(nx,ny) , mtime(12), mustar(nx,ny,12)
       real (kind=4) x4(nx), y4(ny)
       real (kind=4) dpm,   mask(nx,ny)
       integer (kind=4) k,n
       real :: factor,rho_air = 1.29, rho_wat = 1026.
#include "netcdf.inc"
       integer :: ncid,txid,dims(4),lonid,latid,iret
       character*80 name
       integer :: nx2,ny2
       real (kind=4) ::  spv
       real (kind=4), allocatable :: x2(:),y2(:),ustar2(:,:),tx(:)
       real (kind=4), allocatable :: ustar3(:,:),nmask(:,:)
       real, allocatable   :: work_t(:)
       integer :: init_t=1, lenw_t,nn
       type(time_type) :: t1

       open(10,file='data/ecmwf/ecmwf_BB.ustar',form='unformatted')
       do n=1,12
        read (10); read (10) s, dpm, k, k, k, x4,y4,u4(:,:)
        x=x4;y=y4; mustar(:,:,n)=u4
       enddo
       close(10)
c    convert ustar
       factor = 1.e+06 * (rho_air/rho_wat)**1.5
       mustar = factor*mustar
!    apply land mask
       open(10,file='data/ecmwf/ecmwf_BB.mask',form='unformatted')
       read (10) ; read (10); read (10) mask; close(10)
       do n=1,12
        where (mask==0.) mustar(:,:,n) = spval
       enddo
       mtime = (/31,28,31,30,31,30,31,31,30,31,30,31/)
c
c   read land/sea mask for NCEP data
c
       iret=nf_open('data/ncep/lsmask.19294.nc',NF_nowrite,ncid)
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=nf_inq_varid(ncid,'lsmask',txid)
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
       iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx2)
       iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny2)
       allocate( nmask(nx2,ny2) )
       iret= nf_get_vara_real (ncid,txid ,(/1,1,1/),(/nx2,ny2,1/),nmask)
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       call ncclos (ncid, iret)
c
c    read NCEP anomalies
c
       iret=nf_open('data/ncep/ustar_anom.cdf',NF_nowrite,ncid)
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=nf_inq_varid(ncid,'ustar',txid)
       iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
       iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx2)
       iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny2)
       nx2=nx2+2 ! append some rows
       ny2=ny2+4

       allocate(x2(nx2),y2(ny2), ustar2(nx2,ny2),tx(nx2),ustar3(nx,ny))
       name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
       iret=NF_INQ_VARID (NCID, name, lonid)
       name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
       iret=NF_INQ_VARID (NCID, name, latid)

       iret= nf_get_vara_real (ncid,lonid ,1,nx2-2,x2(2:nx2-1))
       x2(nx2)=2*x2(nx2-1)-x2(nx2-2)
       x2(1)=x2(2)-(x2(3)-x2(2))
       iret= nf_get_vara_real (ncid,latid ,1,ny2-4,y2(3:ny2-2)) 
       y2(ny2-1)=2*y2(ny2-2)-y2(ny2-3)
       y2(2)=y2(3)-(y2(4)-y2(3))
       y2(ny2)=2*y2(ny2-1)-y2(ny2-2)
       y2(1)=y2(2)-(y2(3)-y2(2))
       y2=y2(ny2:1:-1)
       iret = nf_get_att_real(ncid, txid, 'missing_value', spv)

c    read each month and interpolate on ECMWF grid
       print*,' interpolating data '
       lenw_t = 10*max(nx,ny)
       allocate( work_t(lenw_t) )
       t1 =  set_date(1948, 1, 1, 0, 0, 0)
       toffset=t1%days
       do n=1,nt
        ustar2=spv
        iret= nf_get_vara_real (ncid,txid ,
     &         (/1,1,n/),(/nx2-2,ny2-4,1/),ustar2(2:nx2-1,3:ny2-2))

        where( nmask == -1) ustar2(2:nx2-1,3:ny2-2) = spv
c     acount for appended rows and switch latitudinal ordering
        ustar2(nx2,:)=ustar2(2,:)
        ustar2(1,:)  =ustar2(nx2-1,:)
        ustar2(:,:) = ustar2(:,ny2:1:-1)
c    convert special values
        where( ustar2 == spv) ustar2=spval
c     interpolate
        call tranlon (ustar2(1,1), nx2, nx2-2, ny2, x2, x,nx, tx)
        call ctf (ustar2(1,1), nx2, ny2, tx, y2, 
     &                  ustar3(1,1),nx, ny, 1, nx, 1, ny,
     &                  x, y, init_t, work_t, lenw_t,spval)
        init_t=0 ! we can use the same interpolation weights
c     add mean
        nn = mod(n-1,12)+1
        where( ustar3 == spval) ustar3=0.
        ustar(:,:,n)=mustar(:,:,nn)+ustar3(:,:)
        where (mask==0.) ustar(:,:,n) = spval
        where (mask/=0. .and. ustar(:,:,n) < 0.0) ustar(:,:,n) = 0.0
        time(n)=mtime(nn)
       enddo
       call ncclos (ncid, iret)
       print*,' done '
       deallocate(x2,y2,tx,ustar2,ustar3,work_t,nmask)
       end subroutine read_ustar

#elif defined read_ustar_omip
c
c--------------------------------------------------------
c    use a daily climatology from the OMIP
c    wind stress data for friction velocity.
c--------------------------------------------------------
      subroutine read_dims_ustar(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4)
      print*,' reading daily friction velocity data from OMIP'
      iret=nf_open('data/OMIP365/north_south_stress.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'stress_y',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)  
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt)  
      call ncclos (ncid, iret)
      end subroutine read_dims_ustar



      subroutine read_ustar(ustar,x,y,nx,ny,nt,spval,time,toffset)
      implicit none
      integer :: nx,ny,nt
      real :: spval, ustar(nx,ny,nt) ,x(nx), y(ny),time(nt),toffset
#include "netcdf.inc"
      integer :: iret,ncid,txid,tyid,dims(4),aid,n,i,j
      character (len=80) :: name
#ifdef netcdf_real4
      real (kind=4),allocatable :: taux(:,:,:),tauy(:,:,:),mask(:,:)
      real (kind=4) :: sp
#else
      real ,allocatable :: taux(:,:,:),tauy(:,:,:),mask(:,:)
      real  :: sp
#endif
      real :: rho_air = 1.29, rho_wat = 1026.
      allocate( taux(nx,ny,nt), tauy(nx,ny,nt), mask(nx,ny) )
c   read land mask
      mask=1.0
      iret=nf_open('data/OMIP/land_sea_mask.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'lsm',txid)
      iret= nf_get_vara_real(ncid,txid,(/1,1/),(/nx,ny/),mask)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      call ncclos (ncid, iret)
c   read taux
      iret=nf_open('data/OMIP365/east_west_stress.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'stress_x',txid)
      iret= nf_get_vara_real(ncid,txid,(/1,1,1/),(/nx,ny,nt/),taux)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret= nf_get_att_real (ncid,txid ,'missing_value',sp)
      do n=1,nt
        where ( taux(:,:,n) == sp ) taux(:,:,n)=spval
        where ( mask == 0.0  )      taux(:,:,n)=spval
      enddo
      ! read grid info
      iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
      iret=nf_inq_varid(ncid,name,aid)
      iret= nf_get_vara_real(ncid,aid ,1,nx,tauy(:,1,1))
      x=tauy(:,1,1)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
      iret=nf_inq_varid(ncid,name,aid)
      iret= nf_get_vara_real (ncid,aid ,1,ny,tauy(1,:,1)) 
      y=tauy(1,:,1)
      call ncclos (ncid, iret)
c   read tauy
      iret=nf_open('data/OMIP365/north_south_stress.nc',NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'stress_y',txid)
      iret= nf_get_vara_real(ncid,txid,(/1,1,1/),(/nx,ny,nt/),tauy)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret= nf_get_att_real (ncid,txid ,'missing_value',sp)
      do n=1,nt
        where ( tauy(:,:,n) == sp ) tauy(:,:,n)=spval
        where ( mask == 0.0  )      tauy(:,:,n)=spval
      enddo
      call ncclos (ncid, iret)
c   ustar**2 = tau/rho_w    ustar**3 = (tau/rho_w)**3/2
!    reference to dens. in water and convert (m/s)**3 to (cm/s)**3
      do n=1,nt
        do i=1,nx
        do j=1,ny
         if ( mask(i,j) /= 0.0 .and. taux(i,j,n)/=spval
     &              .and.tauy(i,j,n)/=spval )   then
          taux(i,j,n)= (taux(i,j,n)**2 + tauy(i,j,n)**2 )**0.5 
          ustar(i,j,n) = 1.e6*(taux(i,j,n)/rho_wat)**1.5
         else 
          ustar(i,j,n) =  spval
         end if
       enddo
       enddo
      enddo
c   flip second dimension up/down
      toffset=0.0
      y=y(ny:1:-1)
      do n=1,nt
       ustar(:,:,n)=ustar(:,ny:1:-1,n)
       time(n)=1.0
      enddo
      deallocate(taux,tauy,mask)
      end subroutine read_ustar

#elif defined read_ustar_T511
c
c--------------------------------------------------------
c    use daily data from ECMWF/T511
c    wind stress data for friction velocity.
c--------------------------------------------------------
      subroutine read_dims_ustar(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4),nn
      character (len=80) :: file
      integer :: start_y, fraction_y, nfraction_y,full_y
      common / flx_ncep/ start_y,fraction_y,nfraction_y,full_y
      namelist /T511_time/ start_y,fraction_y,nfraction_y

      call getunit(ncid,'namelist.T511_time','fsr')
      read(ncid, nml = T511_time) 
      close(ncid)

      print*,' reading daily friction velocity data from ECMWF/T511'
      print*,' year ',start_y,' in ',fraction_y,' fractions'
      print*,' current fraction ',nfraction_y
      write(file,'("data/ecmwf_T511/EWSS_",i4,".cdf")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'EWSS',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)  
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nn)  
      call ncclos (ncid, iret)
      full_y = nn
      print*,' days in year:',full_y
      if (mod(full_y,fraction_y)/=0) then
          print*,'fraction must be integer divider of 365'
          print*,'but ', mod(full_y,fraction_y)
          stop
      endif
      nt=full_y/fraction_y ! must be an integer
      end subroutine read_dims_ustar



      subroutine read_ustar(ustar,x,y,nx,ny,nt,spval,time,toffset)
      use time_type_module
      implicit none
      integer :: nx,ny,nt
      real :: spval, ustar(nx,ny,nt) ,x(nx), y(ny),time(nt),toffset
#include "netcdf.inc"
      integer :: iret,ncid,txid,tyid,dims(4),aid,n,i,j,varid
      integer :: corner(4),edges(4),xid,yid,nn,n1,n2
      character (len=80) :: name,file
#ifdef netcdf_real4
      real (kind=4),allocatable :: taux(:,:,:),tauy(:,:,:),mask(:,:)
      real (kind=4) :: sp
      real (kind=4),allocatable :: v4(:,:,:)
      real (kind=4) :: add_offset,scale_factor
#else
      real ,allocatable :: taux(:,:,:),tauy(:,:,:),mask(:,:)
      real  :: sp
      real ,allocatable :: v4(:,:,:)
      real  :: add_offset,scale_factor
#endif
      real :: mask2(nx,ny)
      type(time_type) :: t1,t2,t3
      real :: rho_air = 1.29, rho_wat = 1026., fxa
      integer :: start_y, fraction_y, nfraction_y,full_y
      common / flx_ncep/ start_y,fraction_y,nfraction_y,full_y

      allocate( taux(nx,ny,nt), tauy(nx,ny,nt), mask(nx,ny) )

c   read land mask
      iret=nf_open('data/ecmwf_T511/land-sea-mask.nc.new',
     &              NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'tcw',varid)
      corner=(/1,1,1,1/); edges = (/nx,ny,1,1/)
      allocate( v4(nx,ny,1) )
      iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
      add_offset=0.0
      iret= nf_get_att_real (ncid,varid ,'add_offset',add_offset)
      scale_factor=1.0
      iret= nf_get_att_real (ncid,varid ,'scale_factor',scale_factor)
      do j=1,ny; do i=1,nx
        mask(i,j)=0.0
        fxa = v4(i,ny-j+1,1)*scale_factor + add_offset
        if (fxa.gt. 0.1)  mask(i,j)=1.
        mask2(i,j)=0.0
        if (fxa.gt. 0.5)  mask2(i,j)=1.
      enddo; enddo
      deallocate(v4)
      call ncclos (ncid, iret)

c     we have to manipulate the land mask a bit
      mask(705:712,288:300)=mask2(705:712,288:300) ! this is Scotland
      mask(635:645,306:312)=mask2(635:645,306:312) ! East greenland
      mask(585:600,310:318)=mask2(585:600,310:318)

      ! read grid info
      iret=nf_open('data/ecmwf_T511/land-sea-mask.nc',NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'lsm',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),file) 
      iret=nf_inq_varid(ncid,file,xid)
      allocate(v4(nx,1,1) )
      iret= nf_get_vara_real (ncid,xid ,1,nx,v4); x=v4(1:nx,1,1)
      deallocate(v4)
      file='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),file) 
      iret=nf_inq_varid(ncid,file,yid)
      allocate(v4(ny,1,1) )
      iret= nf_get_vara_real (ncid,yid ,1,ny,v4); 
      do j=1,ny
           y(j)=v4(ny-j+1,1,1)
      enddo
      deallocate(v4)
      call ncclos (ncid, iret)
      ! time vector
      t1 =  set_date(start_y, 1, 1, 0, 0, 0)
      toffset=t1%days + (nfraction_y-1)*(full_y/fraction_y)
      write(file,'("data/ecmwf_T511/EWSS_",i4,".cdf")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'EWSS',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      call ncclos (ncid, iret)
      do n=1,nt
        time(n) =  365.0/(full_y*1.0)
      enddo

c   read taux and tauy

      print*,' reading year ',start_y
      call sub_flush(6)
      write(file,'("data/ecmwf_T511/EWSS_",i4,".cdf")') start_y
      iret=nf_open(file,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'EWSS',varid)
      iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nn)
      n1=1+(nfraction_y-1)*(full_y/fraction_y)
      nn=nfraction_y*(full_y/fraction_y)
      corner=(/1,1,n1,1/); edges = (/nx,ny,nt,1/)
      iret= nf_get_vara_real (ncid,varid ,corner,edges,taux)
      call ncclos (ncid, iret)

      write(file,'("data/ecmwf_T511/NSSS_",i4,".cdf")') start_y
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_open(file,NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'NSSS',varid)
      iret= nf_get_vara_real (ncid,varid ,corner,edges,tauy)
      call ncclos (ncid, iret)

c   ustar**2 = tau/rho_w    ustar**3 = (tau/rho_w)**3/2
!    reference to dens. in water and convert (m/s)**3 to (cm/s)**3
      do n=1,nt
        do i=1,nx
        do j=1,ny
         if ( mask(i,j) /= 1.0) then 
          taux(i,j,n)= (taux(i,j,n)**2 + tauy(i,j,n)**2 )**0.5 
          ustar(i,j,n) = 1.e6*(taux(i,j,n)/rho_wat)**1.5
         else 
          ustar(i,j,n) =  spval
         end if
       enddo
       enddo
      enddo
      deallocate(taux,tauy,mask)
      end subroutine read_ustar


#endif
