#include "options.inc"

c
c--------------------------------------------------------
c  Templates to read surface tracer fluxes from dataset
c  stored in external files.
c
c  Notes:
c   heat flux should be in cal/s/cm**2
c   restoring temperature should be in deg C
c   restoring coeff. for temp. should be in cal/s/cm**2/K
c   salt flux should be in (psu/1000-35)/s/cm**2
c   restoring salinity should be in psu/1000-35
c   restoring coeff. for sal. should be in cm/s
c
c   Two subroutines are defined:
c   "read_dims_tr" and "read_tr" both called in prep_sbc (in prep_main.F) 
c   Each tracer can be treated different, i.e. could
c   have different grids and is  referenced by "ntr"
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_density_flx_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_tr(nx,ny,nt,ntr)
       use prep_module
       implicit none
       integer nx,ny,nt,ntr
       character*32 stamp
       real :: dpm
       integer :: io
#include "netcdf.inc"
       integer :: iret,ncid,vid,ndims,dims(4)
       if (ntr == 1 ) then
        print*,' reading heat flux boundary conditions from'
        print*,' Barniers ECMWF analysis, which were also used'
        print*,' in the DYNAMO project'
        nt=12
        call getunit(io,'data/ecmwf/ecmwf_BB.qnet','usr ieee')
        read (io); read (io) stamp, dpm, nx, ny
        close(io)
       elseif (ntr == 2) then
        print*,' reading salt flux boundary condition from'
        print*,' monthly climatological SSS (file ts.dta) '
        print*,' and using a fixed restoring time scale '
#ifdef read_salt_fl_from_E5
        print*,' reading salt flux data from E5 exp. '
#endif
        nt=12; nx=imt; ny=jmt
#if defined  prep_for_LEDWELL || defined prep_for_OBC_tracer
       elseif (ntr >= 3) then
        print*,' for tracer #',ntr,' there is no flux'
c     however read some dummy data in any case
        iret=nf_open('data/U10/u10.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'u10',vid)
        iret=NF_INQ_VARNDIMS(NCID, vid, ndims)  
        iret=NF_INQ_VARDIMID (NCID,vid,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)
#endif

#ifdef  prep_for_CFC 
       elseif (ntr == 3) then
        print*,' for tracer 3 we are using a climatology'
        print*,' of U(10m) for the restoring time scale '
        iret=nf_open('data/U10/u10.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'u10',vid)
        iret=NF_INQ_VARNDIMS(NCID, vid, ndims)  
        iret=NF_INQ_VARDIMID (NCID,vid,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)
#endif

#ifdef  prep_for_NPZD
       elseif (ntr >= 3 .and. ntr<=6 ) then
        print*,' for passive tracer in the NPZD model ',
     &         'we are using a climatology of solar radiation'
        nt=12
        call getunit(io,'data/ecmwf/ecmwf_BB.qsol','usr ieee')
        read (io); read (io) stamp, dpm, nx, ny
        close(io)
       elseif (ntr == 7 .or. ntr==8 
#ifdef  prep_for_ACO2
     &       .or. ntr >= 9  
#endif
     &               ) then
        print*,' for DIC/O2 we are using a climatology'
        print*,' of U(10m) for the restoring time scale '
#ifdef read_piston_velocity_U10
        print*,' reading from file data/U10/u10.cdf'
        iret=nf_open('data/U10/u10.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'u10',vid)
        iret=NF_INQ_VARNDIMS(NCID, vid, ndims)  
        if (ndims/=3) then
         print*,' number of dimensions for variable not equal 3'
         print*,ndims
         call halt_stop(' in read_dims_tr (NPZD) ')
        endif
        iret=NF_INQ_VARDIMID (NCID,vid,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)
#endif
#ifdef read_piston_velocity_WSPD
        print*,' reading from file data/U10/wspd_mcycle.cdf'
        iret=nf_open('data/U10/wspd_mcycle.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'wspd',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx);nx=nx+2
        iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
        iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt)
        call ncclos (ncid, iret)
        print*,' this is converted to piston velocity during model run'
        print*,' and time series of pressure for DIC boundary condition'
#endif
#ifndef  prep_for_ACO2
       elseif (ntr ==9 ) then
c     this is for dissolved organic carbon 
        print*,' for passive tracer in the NPZD model ',
     &         'we are using a climatology of solar radiation'
        nt=12
        call getunit(io,'data/ecmwf/ecmwf_BB.qsol','usr ieee')
        read (io); read (io) stamp, dpm, nx, ny
        close(io)
#endif
#endif
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_dims_tr




       subroutine read_tr(tr_clim,tr_rest,tr_flux,
     &                    x,y,time,toffset,nx,ny,nt,ntr)
c
       ! returns inverse restoring time scale for tracers 
       ! number ntr in tr_rest
       ! and climatology to restore against in tr_clim
       ! nx,ny are horizontal dimensions of the data
       ! nt the time steps,
       ! x/y the longitude/latitude of the grid
       ! should mask data with spval
       ! return also net fluxes in Barniers heat flux formulation
       ! or alternatively the shortwave solar radiation in the net flux array
       ! if available 
       ! Units: mom units for tr_clim and tr_rest
c       
       use prep_module
       implicit none
       integer nx,ny,nt,ntr,m
       real tr_clim(nx,ny,nt) ,x(nx), y(ny)
       real tr_rest(nx,ny,nt) 
       real tr_flux(nx,ny,nt) ,time(nt),toffset
       character*32 stamp
       real, allocatable, dimension(:,:,:) :: qnet,qnec,airt,qsol
       real, allocatable, dimension(:,:)   :: mask
       real :: dpm
       integer  k,n,io,i,j
       real,parameter ::rhonull=1.026, Cpw=0.955, dimfac=0.2389e-4
       real,parameter ::qfac = dimfac/(rhonull*Cpw)
#include "netcdf.inc"
       integer :: iret,ncid,vid,dims(4),lonid,latid,corner(4),edges(4)
       character (len=80) :: dname
#ifdef netcdf_real4
      real (kind=4) :: v4(nx,ny) ,x4(nx),y4(ny),sp
#else
      real :: v4(nx,ny) ,x4(nx),y4(ny),sp
#endif
       time = (/31,28,31,30,31,30,31,31,30,31,30,31/)
       toffset=0.0

       if (ntr == 1) then
c     SBC for heat comes from Barniers formulation
        allocate( qnet(nx,ny,nt), qnec(nx,ny,nt), airt(nx,ny,nt))
        allocate( mask(nx,ny), qsol(nx,ny,nt))
        call getunit(io,'data/ecmwf/ecmwf_BB.qnet','usr ieee')
        do n=1,nt
         read (io); read (io) stamp, dpm, k, k, k, x,y,qnet(:,:,n)
         time(n)=dpm
        enddo
        close(io)
        call getunit(io,'data/ecmwf/ecmwf_BB.qnec','usr ieee')
        do n=1,nt
         read (io); read (io) stamp, dpm, k, k, k, x,y,qnec(:,:,n)
        enddo
        close(io)
        call getunit(io,'data/ecmwf/ecmwf_BB.airt','usr ieee')
        do n=1,nt
         read (io); read (io) stamp, dpm, k, k, k, x,y,airt(:,:,n)
        enddo
        close(io)
        call getunit(io,'data/ecmwf/ecmwf_BB.qsol','usr ieee')
        do n=1,nt
         read (io); read (io) stamp, dpm, k, k, k, x,y,qsol(:,:,n)
        enddo
        close(io)
        where (qnec/=0.) 
             tr_clim(:,:,:) = airt + qnet/qnec
c             tr_flux(:,:,:) = qnet  ! store net heat flux
             tr_flux(:,:,:) = qsol/41868.0  ! no store shortwave radiation in the array
        elsewhere
             tr_clim(:,:,:) = spval
        end where
        tr_rest(:,:,:) = qnec*qfac
!     land apply mask
        call getunit(io,'data/ecmwf/ecmwf_BB.mask','usr ieee')
        read (io) ; read (io); read (io) mask; close(io)
        do n=1,nt
         where (mask==0.) tr_clim(:,:,n) = spval
         where (mask==0.) tr_rest(:,:,n) = spval
         where (mask==0.) tr_flux(:,:,n) = spval
        enddo
        deallocate( qnet, qnec, airt,mask,qsol)
       elseif (ntr == 2) then
c    surface salt flux for salinity
#ifdef read_salt_fl_from_E5
c      read salt flux from previous model solution
c      (to implement mixed boundary conditions)
        iret=nf_open(
     &'/users/model/ceden/dov3/N_atl/E5_bbl/salt_fl_mmean_ufo.cdf'
     &              , NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'salt_fl',vid)
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,n,1/); edges=(/nx,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4)
         where( v4 /= sp) v4 = v4/10.  ! convert to cgs
         where( v4 == sp) v4 = -1000.0
         tr_flux(:,:,n)=v4
        enddo
        where( tr_flux == -1000.0 ) tr_flux = spval
        call ncclos (ncid, iret)
#else
        tr_flux=0. ! unknown otherwise 
c      (flux is not used for restoring formulation)
#endif
        tr_rest=1./15.0*dzt(1)/86400.
c        tr_rest=1./15.0*dzt(1)/86400.*10.0
c     read restoring salinity from previously interpolated data
c     which are also used for the initial conditions, sponge layers,  etc.
        call getunit(io,'ts.dta','usr ieee')
        do n=1,nt
         read(io) ; read(io) tr_clim(:,:,n) ! spval is already correct
         where( kmt == 0) tr_clim(:,:,n)=spval
         do m=3,number_tr
          read(io) ! tracer
         enddo
         do k=2,km; do m=1,number_tr
           read(io) ! rest of the levels
        enddo; enddo; enddo
        close(io)
        x=xt;y=yt
#ifdef limAmazoneSSS
c     this is a bug fix for spurious dispersion error in the Amazon river region
        do j = 1,ny
         if (yt(j).ge.-10. .and. yt(j).le.10.) then
          tr_clim = max(tr_clim,-0.005)
         endif
        enddo
#endif
#if defined  prep_for_LEDWELL || defined prep_for_OBC_tracer
       elseif (ntr >= 3) then
c    SBC for tracer 3 is no flux but read dummy data
        tr_clim=0.; tr_flux=0.; tr_rest=0.
        iret=nf_open('data/U10/u10.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'u10',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),dname) 
        iret=nf_inq_varid(ncid,dname,lonid)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),dname) 
        iret=nf_inq_varid(ncid,dname,latid)
        iret= nf_get_vara_real (ncid,lonid ,1,nx,x4); x=x4+360.
        iret= nf_get_vara_real (ncid,latid ,1,ny,y4); y=y4
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        call ncclos (ncid, iret)
#endif

#ifdef  prep_for_CFC 
       elseif (ntr == 3) then
c     SBC for CFC11 and CFC12 (wind in 10m)
        tr_clim=0.; tr_flux=0.
        tr_rest=1./15.0*dzt(1)/86400.
        iret=nf_open('data/U10/u10.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'u10',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),dname) 
        iret=nf_inq_varid(ncid,dname,lonid)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),dname) 
        iret=nf_inq_varid(ncid,dname,latid)
        iret= nf_get_vara_real (ncid,lonid ,1,nx,x4); x=x4+360.
        iret= nf_get_vara_real (ncid,latid ,1,ny,y4); y=y4
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,n,1/); edges=(/nx,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4)
         where( v4 == sp) v4 = -1000.0
         tr_rest(:,:,n)=v4
        enddo
        where( tr_rest== -1000.0 ) tr_rest = spval
        call ncclos (ncid, iret)
#endif

#ifdef  prep_for_NPZD
       elseif (ntr == 3) then
c    SBC for nutrient growing function (solar radiation)
        tr_clim=0.; tr_flux=0.; tr_rest=1./15.0*dzt(1)/86400.
        call getunit(io,'data/ecmwf/ecmwf_BB.qsol','usr ieee')
        do n=1,nt
         read (io); read (io) stamp, dpm, k, k, k, x,y,tr_clim(:,:,n)
        enddo
        close(io)
        allocate( mask(nx,ny))
        call getunit(io,'data/ecmwf/ecmwf_BB.mask','usr ieee')
        read (io) ; read (io); read (io) mask; close(io)
        do n=1,nt
         where (mask==0.) tr_clim(:,:,n) = spval
         where (mask==0.) tr_rest(:,:,n) = spval
         where (mask==0.) tr_flux(:,:,n) = spval
        enddo
        deallocate( mask)
       elseif (ntr > 3 .and. ntr <=6) then
c       phyto, zoo plankton and detritus, no surface fluxes
        tr_clim=0.0; tr_flux=0.0; tr_rest=0.0
c     read dummy values, only for grid
        call getunit(io,'data/ecmwf/ecmwf_BB.qsol','usr ieee')
        read (io); read (io) stamp, dpm, k, k, k, x,y
        close(io)
       elseif (ntr ==7 .or. ntr == 8
#ifdef  prep_for_ACO2
     &       .or. ntr >= 9  
#endif
     &             ) then
c       SBC for DIC: piston velocity  depends on u10 (same as for CFC)
c       SBC for O2: same as for DIC
        tr_clim=0.; tr_flux=0.
#ifdef read_piston_velocity_U10
        iret=nf_open('data/U10/u10.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'u10',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),dname) 
        iret=nf_inq_varid(ncid,dname,lonid)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),dname) 
        iret=nf_inq_varid(ncid,dname,latid)
        iret= nf_get_vara_real (ncid,lonid ,1,nx,x4); x=x4+360.
        iret= nf_get_vara_real (ncid,latid ,1,ny,y4); y=y4
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,n,1/); edges=(/nx,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4)
         do i=1,nx; do j=1,ny
            if( v4(i,j) == sp) v4(i,j) = -1000.0
            tr_rest(i,j,n)=v4(i,j)
            if( tr_rest(i,j,n)== -1000.0 ) tr_rest(i,j,n) = spval
         enddo; enddo
        enddo
        call ncclos (ncid, iret)
#endif
#ifdef read_piston_velocity_WSPD
        iret=nf_open('data/U10/wspd_mcycle.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'wspd',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),dname) 
        iret=nf_inq_varid(ncid,dname,lonid)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),dname) 
        iret=nf_inq_varid(ncid,dname,latid)
        iret= nf_get_vara_real (ncid,lonid ,1,nx-2,x4(2:nx-1)); x=x4
        x(nx)=2*x(nx-1)-x(nx-2)
        x(1)=x(2)-(x(3)-x(2))
        iret= nf_get_vara_real (ncid,latid ,1,ny,y4); y=y4
        y=y(ny:1:-1)
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,n,1/); edges=(/nx-2,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4(2:nx-1,:))
         v4=v4(:,ny:1:-1)
         v4(nx,:)=v4(2,:)
         v4(1,:)=v4(nx-1,:)
         do i=1,nx; do j=1,ny
            if( v4(i,j) == sp) v4(i,j) = -1000.0
            if( v4(i,j) /= sp) v4(i,j) = v4(i,j)*0.01+202.65
            tr_rest(i,j,n)=v4(i,j)
            if( tr_rest(i,j,n)== -1000.0 ) tr_rest(i,j,n) = spval
         enddo; enddo
        enddo
        call ncclos (ncid, iret)
c     now read slp and store in tr_clim
        iret=nf_open('data/ncep/slp.mon.ltm.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'slp',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,n,1/); edges=(/nx-2,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4(2:nx-1,:))
         v4=v4(:,ny:1:-1)
         v4(nx,:)=v4(2,:)
         v4(1,:)=v4(nx-1,:)
         do i=1,nx; do j=1,ny
            if( v4(i,j) == sp) v4(i,j) = -1000.0
            if( v4(i,j) /= sp) v4(i,j) = v4(i,j)
            tr_clim(i,j,n)=v4(i,j)
            if( tr_clim(i,j,n)== -1000.0 ) tr_clim(i,j,n) = spval
         enddo; enddo
        enddo
        call ncclos (ncid, iret)
#endif

#ifndef  prep_for_ACO2
       elseif (ntr ==9 ) then ! no flux condition for DOC
        tr_clim=0.; tr_flux=0.; tr_rest=0.0
c     read dummy values, only for grid
        call getunit(io,'data/ecmwf/ecmwf_BB.qsol','usr ieee')
        read (io); read (io) stamp, dpm, k, k, k, x,y
        close(io)
#endif
#endif
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_tr



#elif defined read_density_flx_pml_ncep
c
c--------------------------------------------------------
c     Density flux for PML formulation from mean NCEP
c--------------------------------------------------------
c
       subroutine read_dims_tr(nx,ny,nt,ntr)
       use prep_module
       implicit none
       integer nx,ny,nt,ntr
       character*32 stamp
       real :: dpm
       integer :: io
#include "netcdf.inc"
       integer :: iret,ncid,vid,ndims,dims(4)
       if (ntr == 1 ) then
        print*,' for tracer 1 we are computing variables'
        print*,' necessary for Seagers PML '
        iret=nf_open('data/ncep_for_pml/shum.mon.ltm.nc',
     &                NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'shum',vid)
        iret=NF_INQ_VARNDIMS(NCID, vid, ndims)  
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
        iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
        iret=NF_INQ_DIMLEN (NCID, DIMs(4), nt)
        call ncclos (ncid, iret)
       elseif (ntr == 2) then
        print*,' reading salt flux boundary condition from'
        print*,' monthly climatological SSS (file ts.dta) '
        print*,' and using a fixed restoring time scale '
        nt=12; nx=imt; ny=jmt
       elseif (ntr == 3) then
        print*,' for tracer 3 we are computing variables'
        print*,' necessary for Seagers PML '
        iret=nf_open('data/ncep_for_pml/uwnd.mon.ltm.nc',
     &                NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'uwnd',vid)
        iret=NF_INQ_VARNDIMS(NCID, vid, ndims)  
        iret=NF_INQ_VARDIMID (NCID,vid,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)
       elseif (ntr == 4) then
        print*,' for tracer 4 we are also computing variables'
        print*,' necessary for Seagers PML '
        iret=nf_open('data/ncep_for_pml/nswrs.mon.ltm.nc',
     &                NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'nswrs',vid)
        iret=NF_INQ_VARNDIMS(NCID, vid, ndims)  
        iret=NF_INQ_VARDIMID (NCID,vid,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)
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_dims_tr




       subroutine read_tr(tr_clim,tr_rest,tr_flux,
     &                    x,y,time,toffset,nx,ny,nt,ntr)
c
       ! returns inverse restoring time scale for tracers 
       ! number ntr in tr_rest
       ! and climatology to restore against in tr_clim
       ! nx,ny are horizontal dimensions of the data
       ! nt the time steps,
       ! x/y the longitude/latitude of the grid
       ! should mask data with spval
       ! return also net fluxes in Barniers heat flux formulation
       ! or alternatively the shortwave solar radiation in the net flux array
       ! if available 
       ! Units: mom units for tr_clim and tr_rest
c       
       use prep_module
       implicit none
       integer nx,ny,nt,ntr,m
       real tr_clim(nx,ny,nt) ,x(nx), y(ny)
       real tr_rest(nx,ny,nt) 
       real tr_flux(nx,ny,nt) ,time(nt),toffset
       character*32 stamp
       real, allocatable, dimension(:,:,:) :: qnet,qnec,airt,qsol
       real, allocatable, dimension(:,:)   :: mask
       real :: dpm
       integer  k,n,io,i,j
       real,parameter ::rhonull=1.026, Cpw=0.955, dimfac=0.2389e-4
       real,parameter ::qfac = dimfac/(rhonull*Cpw)
#include "netcdf.inc"
       integer :: iret,ncid,vid,dims(4),lonid,latid,corner(4),edges(4)
       character (len=80) :: dname
#ifdef netcdf_real4
      real (kind=4) :: v4(nx,ny) ,x4(nx),y4(ny),sp
#else
      real :: v4(nx,ny) ,x4(nx),y4(ny),sp
#endif
       time = (/31,28,31,30,31,30,31,31,30,31,30,31/)
       toffset=0.0

       if (ntr == 1) then

        tr_clim=0.; tr_flux=0.; tr_rest=0.
c
c       specific humidity to  sst_rest
c
        iret=nf_open('data/ncep_for_pml/shum.mon.ltm.nc',
     &                NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'shum',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),dname) 
        iret=nf_inq_varid(ncid,dname,lonid)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),dname) 
        iret=nf_inq_varid(ncid,dname,latid)
        iret= nf_get_vara_real (ncid,lonid ,1,nx,x4); x=x4
        iret= nf_get_vara_real (ncid,latid ,1,ny,y4); y=y4
        y=y(ny:1:-1)
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,1,n/); edges=(/nx,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4)
         v4=v4(:,ny:1:-1)
         where( v4 == sp) v4 = -1000.0
         tr_rest(:,:,n)=v4
        enddo
        where( tr_rest/= -1000.0 ) tr_rest = tr_rest/1000.0
        where( tr_rest== -1000.0 ) tr_rest = spval
        call ncclos (ncid, iret)

c
c        air temperature to tr_clim
c
        iret=nf_open('data/ncep_for_pml/air.mon.ltm.nc',
     &                NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'air',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,1,n/); edges=(/nx,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4)
         v4=v4(:,ny:1:-1)
         where( v4 == sp) v4 = -1000.0
         tr_clim(:,:,n)=v4
        enddo
        where( tr_clim/= -1000.0 ) tr_clim = tr_clim+273.15
        where( tr_clim== -1000.0 ) tr_clim = spval
        call ncclos (ncid, iret)


        apply_landmask = .false.


       elseif (ntr == 2) then
c    surface salt flux for salinity
        tr_flux=0. ! unknown 
c      (flux is not used for restoring formulation)
        tr_rest=1./15.0*dzt(1)/86400.
c        tr_rest=1./15.0*dzt(1)/86400.*10.0
c     read restoring salinity from previously interpolated data
c     which are also used for the initial conditions, sponge layers,  etc.
        call getunit(io,'ts.dta','usr ieee')
        do n=1,nt
         read(io) ; read(io) tr_clim(:,:,n) ! spval is already correct
         where( kmt == 0) tr_clim(:,:,n)=spval
         do m=3,number_tr
          read(io) ! tracer
         enddo
         do k=2,km; do m=1,number_tr
           read(io) ! rest of the levels
        enddo; enddo; enddo
        close(io)
        x=xt;y=yt
#ifdef limAmazoneSSS
c     this is a bug fix for spurious dispersion error in the Amazon river region
        do j = 1,ny
         if (yt(j).ge.-10. .and. yt(j).le.10.) then
          tr_clim = max(tr_clim,-0.005)
         endif
        enddo
#endif

c    SBC for Seagers PML 
c      they are stored in str03_flux/clim/rest
       elseif (ntr == 3) then

        tr_clim=0.; tr_flux=0.; tr_rest=0.
c
c      uwind to tr_rest
c
        iret=nf_open('data/ncep_for_pml/uwnd.mon.ltm.nc',
     &                NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'uwnd',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),dname) 
        iret=nf_inq_varid(ncid,dname,lonid)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),dname) 
        iret=nf_inq_varid(ncid,dname,latid)
        iret= nf_get_vara_real (ncid,lonid ,1,nx,x4); x=x4
        iret= nf_get_vara_real (ncid,latid ,1,ny,y4); y=y4
        y=y(ny:1:-1)
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,n,1/); edges=(/nx,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4)
         v4=v4(:,ny:1:-1)
         where( v4 == sp) v4 = -1000.0
         tr_rest(:,:,n)=v4
        enddo
        where( tr_rest== -1000.0 ) tr_rest = spval
        call ncclos (ncid, iret)
c
c      vwind to tr_clim
c
        iret=nf_open('data/ncep_for_pml/vwnd.mon.ltm.nc',
     &                NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'vwnd',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,n,1/); edges=(/nx,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4)
         v4=v4(:,ny:1:-1)
         where( v4 == sp) v4 = -1000.0
         tr_clim(:,:,n)=v4
        enddo
        where( tr_clim== -1000.0 ) tr_clim = spval
        call ncclos (ncid, iret)
c
c      wspd to tr_flux
c
        iret=nf_open('data/ncep_for_pml/wspd.mon.ltm.nc',
     &                NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'wspd',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,n,1/); edges=(/nx,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4)
         v4=v4(:,ny:1:-1)
         where( v4 == sp) v4 = -1000.0
         tr_flux(:,:,n)=v4
        enddo
        where( tr_flux/= -1000.0 ) tr_flux = 
     &                       tr_flux*0.0099999998+202.64999
        where( tr_flux== -1000.0 ) tr_flux = spval
        call ncclos (ncid, iret)

        apply_landmask = .false.

       elseif (ntr == 4) then
        tr_clim=0.; tr_flux=0.; tr_rest=0.
c
c      qsol to tr_flux 4
c
        iret=nf_open('data/ncep_for_pml/nswrs.mon.ltm.nc',
     &                NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'nswrs',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),dname) 
        iret=nf_inq_varid(ncid,dname,lonid)
        dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),dname) 
        iret=nf_inq_varid(ncid,dname,latid)
        iret= nf_get_vara_real (ncid,lonid ,1,nx,x4); x=x4
        iret= nf_get_vara_real (ncid,latid ,1,ny,y4); y=y4
        y=y(ny:1:-1)
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,n,1/); edges=(/nx,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4)
         v4=v4(:,ny:1:-1)
         where( v4 == sp) v4 = -1000.0
         tr_flux(:,:,n)=v4
        enddo
        where( tr_flux/= -1000.0 ) tr_flux = -tr_flux/41868.0
        where( tr_flux== -1000.0 ) tr_flux = spval
        call ncclos (ncid, iret)

c
c      cloud cover to tr_clim
c
        iret=nf_open('data/ncep_for_pml/tcdc_mcycle.cdf',
     &                NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'tcdc',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         corner=(/1,1,n,1/); edges=(/nx,ny,1,1/)
         iret= nf_get_vara_real (ncid,vid ,corner,edges,v4)
         v4=v4(:,ny:1:-1)
         where( v4 == sp) v4 = -1000.0
         tr_clim(:,:,n)=v4
        enddo
        where( tr_clim/= -1000.0 ) tr_clim =(tr_clim*0.1+3276.5)/100
        where( tr_clim== -1000.0 ) tr_clim =spval
        call ncclos (ncid, iret)

        apply_landmask = .false.

       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_tr


#elif defined read_density_flx_ncep_monthly
c
c--------------------------------------------------------
c    use a 3 year monthly mean climatology from the european
c    ECMWF model (same data used in the DYNAMO project)
c    plus anomalies from NCEP/NCAR data for heat flux
c--------------------------------------------------------
c

       module density_flx_ncep_monthly
       integer           :: end_nt = 300*12   ! start and end index of months to be read
       integer,parameter :: start_nt = 25*12+1   ! start and end index of months to be read
       integer,parameter :: start_year = 1948 +(start_nt)/12! start year to be read
       end module density_flx_ncep_monthly

       subroutine read_dims_tr(nx,ny,nt,ntr)
       use prep_module
       use density_flx_ncep_monthly
       implicit none
       integer nx,ny,nt,ntr
       character*32 stamp
       real :: dpm
       integer :: io
#include "netcdf.inc"
       integer :: iret,ncid,vid,ndims,dims(4)
       if (ntr == 1 ) then
        print*,' reading heat flux boundary conditions from'
        print*,' Barniers ECMWF analysis, which were also used'
        print*,' in the DYNAMO project'
        print*,' plus heat flux anomalies from the NCEP/NCAR data'
        call getunit(io,'data/ecmwf/ecmwf_BB.qnet','usr ieee')
        read (io); read (io) stamp, dpm, nx, ny
        close(io)
        iret=nf_open('data/ncep/net_hflx_anom.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'hfl',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,DIMS)
        iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt)
        end_nt=min(end_nt,int(nt/12)*12 )
        nt=end_nt - start_nt + 1
        call ncclos (ncid, iret)

       elseif (ntr == 2) then
        print*,' reading salt flux boundary condition from'
        print*,' monthly climatological SSS (file ts.dta) '
        print*,' and using a fixed restoring time scale '
        nt=12; nx=imt; ny=jmt

#ifdef  prep_for_NPZD

       elseif (ntr == 3 ) then
        print*,' for Nitrate in the NPZD model ',
     &         'we are using NCEP/NCAR solar radiation'
        iret=nf_open('data/ncep/nswrs.mon.mean.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'nswrs',vid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=NF_INQ_VARDIMID (NCID,vid,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)
        nx=nx+2 ! append some rows
        ny=ny+4
        end_nt=min(end_nt,int(nt/12)*12 )
        nt=end_nt - start_nt + 1
        call ncclos (ncid, iret)

       elseif (ntr > 3 .and. ntr<=6 ) then
        print*,' for P/Z/D in the NPZD model ',
     &         'we are using a dummy climatology of solar radiation'
        nt=12
        call getunit(io,'data/ecmwf/ecmwf_BB.qsol','usr ieee')
        read (io); read (io) stamp, dpm, nx, ny
        close(io)

       elseif (ntr == 7 .or. ntr==8 ) then
        print*,' for DIC/O2 we are using a NCEP/NCAR'
        print*,' time series of U(10m) for the restoring time scale '
        print*,' this is converted to piston velocity during model run'
        print*,' and time series of pressure for DIC boundary condition'

        iret=nf_open('data/ncep/wspd.mon.mean.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'wspd',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,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)
        nx=nx+2 ! append some rows
        ny=ny+4
        end_nt=min(end_nt,int(nt/12)*12 )
        nt=end_nt - start_nt + 1
        call ncclos (ncid, iret)

       elseif (ntr ==9 ) then
        print*,' for DOC in the NPZD model ',
     &         'we are using a climatology of solar radiation'
        nt=12
        call getunit(io,'data/ecmwf/ecmwf_BB.qsol','usr ieee')
        read (io); read (io) stamp, dpm, nx, ny
        close(io)

#endif
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_dims_tr




       subroutine read_tr(tr_clim,tr_rest,tr_flux,
     &                    x,y,time,toffset,nx,ny,nt,ntr)
c
       ! returns inverse restoring time scale for tracers 
       ! number ntr in tr_rest
       ! and climatology to restore against in tr_clim
       ! nx,ny are horizontal dimensions of the data
       ! nt the time steps,
       ! x/y the longitude/latitude of the grid
       ! should mask data with spval
       ! return also net fluxes in Barniers heat flux formulation
       ! if available 
       ! Units: mom units for tr_clim and tr_rest
c       
       use time_type_module
       use prep_module
       use density_flx_ncep_monthly
       implicit none
       integer nx,ny,nt,ntr,m
       real tr_clim(nx,ny,nt) ,x(nx), y(ny)
       real tr_rest(nx,ny,nt) 
       real tr_flux(nx,ny,nt) ,time(nt),toffset
       character*32 stamp
       real, allocatable, dimension(:,:,:) :: qnet,qnec,airt
       real, allocatable, dimension(:,:)   :: mask
       real :: dpm
       integer  k,n,io,i,j
       real,parameter ::rhonull=1.026, Cpw=0.955, dimfac=0.2389e-4
       real,parameter ::qfac = dimfac/(rhonull*Cpw)
#include "netcdf.inc"
       integer :: iret,ncid,vid,dims(4),lonid,latid,corner(4),edges(4)
       character (len=80) :: dname
#ifdef netcdf_real4
      real (kind=4) :: v4(nx,ny) ,x4(nx),y4(ny),sp
#else
      real :: v4(nx,ny) ,x4(nx),y4(ny),sp
#endif
      real (kind=4) :: mtime(12),spv,add_offset,scale_factor,spv2
      real (kind=4), allocatable :: x2(:),y2(:),hfl2(:,:),tx(:)
      real (kind=4), allocatable :: hfl3(:,:),nmask(:,:),slp(:,:)
       character*80 name
      integer :: nx2,ny2,varid,ncid2
      real, allocatable   :: work_t(:)
      integer :: init_t=1, lenw_t,nn
      type(time_type) :: t1

       time(1:12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
       mtime = (/31,28,31,30,31,30,31,31,30,31,30,31/)
       toffset=0.0

       if (ntr == 1) then
c     SBC for heat comes from Barniers formulation
        allocate( qnet(nx,ny,12), qnec(nx,ny,12), airt(nx,ny,12))
        allocate( mask(nx,ny))
        call getunit(io,'data/ecmwf/ecmwf_BB.qnet','usr ieee')
        do n=1,12
         read (io); read (io) stamp, dpm, k, k, k, x,y,qnet(:,:,n)
         mtime(n)=dpm
        enddo
        close(io)
        call getunit(io,'data/ecmwf/ecmwf_BB.qnec','usr ieee')
        do n=1,12
         read (io); read (io) stamp, dpm, k, k, k, x,y,qnec(:,:,n)
        enddo
        close(io)
        call getunit(io,'data/ecmwf/ecmwf_BB.airt','usr ieee')
        do n=1,12
         read (io); read (io) stamp, dpm, k, k, k, x,y,airt(:,:,n)
        enddo
        close(io)
        call getunit(io,'data/ecmwf/ecmwf_BB.mask','usr ieee')
        read (io) ; read (io); read (io) mask; close(io)
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',varid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=NF_INQ_VARDIMID (NCID,varid,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,varid,(/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 data
c
        iret=nf_open('data/ncep/net_hflx_anom.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'hfl',varid)
        iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
        nx2=nx2+2 ! append some rows
        ny2=ny2+4
        allocate(x2(nx2),y2(ny2), hfl2(nx2,ny2),tx(nx2),hfl3(nx,ny))
        name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
        iret=NF_INQ_VARID (NCID, name, lonid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
        iret=NF_INQ_VARID (NCID, name, latid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif

        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, varid, 'missing_value', spv)
c
c    read each month and interpolate on ECMWF grid
c
        print*,' interpolating data '
        lenw_t = 10*max(nx,ny)
        allocate( work_t(lenw_t) )
        t1 =  set_date(start_year, 1, 1, 0, 0, 0)
        toffset=t1%days
        do n=1,nt
         hfl2=spv
         iret= nf_get_vara_real (ncid,varid ,
     &     (/1,1,n+start_nt-1/),(/nx2-2,ny2-4,1/),hfl2(2:nx2-1,3:ny2-2))

         where( nmask == -1) hfl2(2:nx2-1,3:ny2-2) = spv
c     acount for appended rows and switch latitudinal ordering
         hfl2(nx2,:)=hfl2(2,:)
         hfl2(1,:)  =hfl2(nx2-1,:)
         hfl2(:,:) = hfl2(:,ny2:1:-1)
c    convert special values
         where( hfl2 == spv) hfl2=spval
c     interpolate
         call tranlon (hfl2(1,1), nx2, nx2-2, ny2, x2, x,nx, tx)
         call ctf (hfl2(1,1), nx2, ny2, tx, y2, 
     &                  hfl3(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( hfl3 == spval) hfl3=0.
         where (qnec(:,:,nn)/=0.) 
          tr_clim(:,:,n)=airt(:,:,nn)+(qnet(:,:,nn)-hfl3)/qnec(:,:,nn)
          tr_flux(:,:,n)=qnet(:,:,nn)-hfl3
         elsewhere
             tr_clim(:,:,n) = spval
         end where
         tr_rest(:,:,n) = qnec(:,:,nn)*qfac
         where (mask==0.) tr_clim(:,:,n) = spval
         where (mask==0.) tr_rest(:,:,n) = spval
         where (mask==0.) tr_flux(:,:,n) = spval
         time(n)=mtime(nn)
        enddo
        call ncclos (ncid, iret)
        print*,' done '
        deallocate(x2,y2,tx,hfl2,hfl3,work_t,nmask)
        deallocate( qnet, qnec, airt,mask)

       elseif (ntr == 2) then
c    surface salt flux for salinity
        tr_flux=0. ! unknown otherwise 
c      (flux is not used for restoring formulation)
        tr_rest=1./15.0*dzt(1)/86400.
c        tr_rest=1./15.0*dzt(1)/86400.*10.0
c     read restoring salinity from previously interpolated data
c     which are also used for the initial conditions, sponge layers,  etc.
        call getunit(io,'ts.dta','usr ieee')
        do n=1,nt
         read(io) ; read(io) tr_clim(:,:,n) ! spval is already correct
         where( kmt == 0) tr_clim(:,:,n)=spval
         do m=3,number_tr
          read(io) ! tracer
         enddo
         do k=2,km; do m=1,number_tr
           read(io) ! rest of the levels
        enddo; enddo; enddo
        close(io)
        x=xt;y=yt

#ifdef  prep_for_NPZD

       elseif (ntr == 3) then
c    SBC for nutrient growing function (solar radiation)
        tr_clim=0.; tr_flux=0.; tr_rest=1./15.0*dzt(1)/86400.
c        call getunit(io,'data/ecmwf/ecmwf_BB.qsol','usr ieee')
c        do n=1,nt
c         read (io); read (io) stamp, dpm, k, k, k, x,y,tr_clim(:,:,n)
c        enddo
c        close(io)
c        allocate( mask(nx,ny))
c        call getunit(io,'data/ecmwf/ecmwf_BB.mask','usr ieee')
c        read (io) ; read (io); read (io) mask; close(io)
c        do n=1,nt
c         where (mask==0.) tr_clim(:,:,n) = spval
c         where (mask==0.) tr_rest(:,:,n) = spval
c         where (mask==0.) tr_flux(:,:,n) = spval
c        enddo
c        deallocate( mask)

c
c   read land/sea mask for NCEP data
c
        print*,' reading NCEP land/sea 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)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=NF_INQ_VARDIMID (NCID,varid,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,varid,(/1,1,1/),(/nx2,ny2,1/),nmask)
        call ncclos (ncid, iret)
c
c      read NCEP data
c
        iret=nf_open('data/ncep/nswrs.mon.mean.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'nswrs',varid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
        nx2=nx2+2 ! append some rows
        ny2=ny2+4
        allocate(x2(nx2),y2(ny2),hfl2(nx2,ny2) )
        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, varid, 'missing_value', spv)
        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)
        x=x2;y=y2
c
c    read each month 
c
        t1 =  set_date(start_year, 1, 1, 0, 0, 0)
        toffset=t1%days
        do n=1,nt
         hfl2=spv
         iret=nf_inq_varid(ncid,'nswrs',varid)
         iret= nf_get_vara_real (ncid,varid ,
     &     (/1,1,n+start_nt-1/),(/nx2-2,ny2-4,1/),hfl2(2:nx2-1,3:ny2-2))
         where( nmask == -1) hfl2(2:nx2-1,3:ny2-2) = spv
c      change sign and account for scaling
         where( hfl2 /= spv) hfl2=-hfl2*scale_factor+add_offset 
c     acount for appended rows and switch latitudinal ordering
         hfl2(nx2,:)=hfl2(2,:)
         hfl2(1,:)  =hfl2(nx2-1,:)
         hfl2(:,:) = hfl2(:,ny2:1:-1)
c    convert special values and units
         where( hfl2 == spv) hfl2=spval
         tr_clim(:,:,n) = hfl2
c     add mean
         nn = mod(n-1,12)+1
         time(n)=mtime(nn)
        enddo
        print*,' done '
        deallocate(x2,y2,hfl2,nmask)
        call ncclos (ncid, iret)

       elseif (ntr > 3 .and. ntr <=6) then
c       phyto, zoo plankton and detritus, no surface fluxes
        tr_clim=0.0; tr_flux=0.0; tr_rest=0.0
        call getunit(io,'data/ecmwf/ecmwf_BB.qsol','usr ieee')
        read (io); read (io) stamp, dpm, k, k, k, x,y
        close(io)

       elseif (ntr ==7 .or. ntr == 8) then
c       SBC for DIC: piston velocity  depends on u10 (same as for CFC)
c       SBC for O2: same as for DIC
        tr_clim=0.; tr_flux=0.
c
c   read land/sea mask for NCEP data
c
        iret=nf_open('data/ncep/land.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'land',varid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=NF_INQ_VARDIMID (NCID,varid,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,varid,(/1,1,1/),(/nx2,ny2,1/),nmask)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        call ncclos (ncid, iret)
c
c      read wind speed data
c
        print*,' reading wind speed data';call sub_flush(6)
        iret=nf_open('data/ncep/wspd.mon.mean.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'wspd',varid)
        iret=NF_INQ_VARDIMID (NCID,varid,DIMS)
        nx2=nx2+2 ! append some rows
        ny2=ny2+4
        allocate(x2(nx2),y2(ny2),hfl2(nx2,ny2),slp(nx2,ny2))
        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, varid, 'missing_value', spv)
        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)
        x=x2;y=y2
c
c      read SLP data
c
        print*,' reading SLP data';call sub_flush(6)
        iret=nf_open('data/ncep/slp.mon.mean.nc',NF_nowrite,ncid2)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid2,'slp',varid)
        iret = nf_get_att_real(ncid2, varid, 'missing_value', spv2)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
c
c    read each month and interpolate on ECMWF grid
c
        t1 =  set_date(start_year, 1, 1, 0, 0, 0)
        toffset=t1%days
        do n=1,nt
         hfl2=spv
         iret=nf_inq_varid(ncid,'wspd',varid)
         iret= nf_get_vara_real (ncid,varid ,
     &    (/1,1,n+start_nt-1/),(/nx2-2,ny2-4,1/),hfl2(2:nx2-1,3:ny2-2))
         where( nmask == -1) hfl2(2:nx2-1,3:ny2-2) = spv
         where( hfl2 /= spv) hfl2=hfl2*scale_factor+add_offset
c     acount for appended rows and switch latitudinal ordering
         hfl2(nx2,:)=hfl2(2,:)
         hfl2(1,:)  =hfl2(nx2-1,:)
         hfl2(:,:) = hfl2(:,ny2:1:-1)
c    convert special values and units
         where( hfl2 == spv) hfl2=spval
         tr_rest(:,:,n) = hfl2

         slp=spv2
         iret=nf_inq_varid(ncid2,'slp',varid)
         iret= nf_get_vara_real (ncid2,varid ,
     &    (/1,1,n+start_nt-1/),(/nx2-2,ny2-4,1/), slp(2:nx2-1,3:ny2-2))
         where( nmask == -1) slp(2:nx2-1,3:ny2-2) = spv2
         where( slp /= spv) slp=slp*1.0
c     acount for appended rows and switch latitudinal ordering
         slp(nx2,:)=slp(2,:)
         slp(1,:)  =slp(nx2-1,:)
         slp(:,:) = slp(:,ny2:1:-1)
c    convert special values and units
         where( slp == spv2) slp=spval
         tr_clim(:,:,n) = slp
c     the time
         nn = mod(n-1,12)+1
         time(n)=mtime(nn)
        enddo
        print*,' done '
        deallocate(x2,y2,hfl2,nmask,slp)
        call ncclos (ncid, iret)


       elseif (ntr ==9 ) then ! no flux condition for DOC
        tr_clim=0.; tr_flux=0.; tr_rest=0.0
        call getunit(io,'data/ecmwf/ecmwf_BB.qsol','usr ieee')
        read (io); read (io) stamp, dpm, k, k, k, x,y
        close(io)
#endif
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_tr



#elif defined read_density_flx_nao_proxy
c
c--------------------------------------------------------
c    use the 3 year monthly mean ECMWF climatology plus
c    NAO related time series of heat flux anomalies
c    (which is given by the product of the monthly NAO-index
c    and the monthly regression patterns of the 
c    NCEP net heat flux anomalies with 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 hflx_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 hflx_nao_reg
c    See also prep_template_tau.F 
c--------------------------------------------------------
c
       subroutine read_dims_tr(nx,ny,nt,ntr)
       use prep_module
       implicit none
       integer nx,ny,nt,ntr
       character*32 stamp
       real :: dpm
       integer :: io
#include "netcdf.inc"
       integer :: iret,ncid,vid,ndims,dims(4)
       integer :: start_y,end_y
       common / hflx_nao_reg/ start_y,end_y

       if (ntr == 1 ) then
        print*,' reading heat flux boundary conditions from'
        print*,' Barniers ECMWF analysis, which were also used'
        print*,' in the DYNAMO project'
        print*,' plus NAO related heat flux anomalies from NCEP/NCAR'
        call getunit(io,'data/ecmwf/ecmwf_BB.qnet','usr ieee')
        read (io); read (io) stamp, dpm, nx, ny
        close(io)
        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',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,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
       elseif (ntr == 2) then
        print*,' reading salt flux boundary condition from'
        print*,' monthly climatological SSS (file ts.dta) '
        print*,' and using a fixed restoring time scale '
#ifdef read_salt_fl_from_E5
        print*,' option read_salt_fl_from_E5 is not working here'
        stop
#endif
        nt=12; nx=imt; ny=jmt
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_dims_tr


       subroutine read_tr(tr_clim,tr_rest,tr_flux,
     &                    x,y,time,toffset,nx,ny,nt,ntr)
c
       ! returns inverse restoring time scale for tracers 
       ! number ntr in tr_rest
       ! and climatology to restore against in tr_clim
       ! nx,ny are horizontal dimensions of the data
       ! nt the time steps,
       ! x/y the longitude/latitude of the grid
       ! should mask data with spval
       ! return also net fluxes in Barniers heat flux formulation
       ! if available 
       ! Units: mom units for tr_clim and tr_rest
c       
       use prep_module
       use time_type_module
       implicit none
       integer nx,ny,nt,ntr,m,nn
       real tr_clim(nx,ny,nt) ,x(nx), y(ny)
       real tr_rest(nx,ny,nt) 
       real tr_flux(nx,ny,nt) ,time(nt),toffset
       character*32 stamp
       real, allocatable, dimension(:,:,:) :: qnet,qnec,airt
       real, allocatable, dimension(:,:)   :: mask
       real :: dpm
       integer  k,n,io,i,j
       real,parameter ::rhonull=1.026, Cpw=0.955, dimfac=0.2389e-4
       real,parameter ::qfac = dimfac/(rhonull*Cpw)
#include "netcdf.inc"
       integer :: iret,ncid,vid,dims(4),lonid,latid,corner(4),edges(4)
       character (len=80) :: dname
#ifdef netcdf_real4
      real (kind=4) :: v4(nx,ny) ,x4(nx),y4(ny),sp
#else
      real :: v4(nx,ny) ,x4(nx),y4(ny),sp
#endif
      integer :: init_t=1, lenw_t,nx2,ny2
      real (kind=4), allocatable :: x2(:),y2(:),hflx(:,:,:),tx(:),
     &               hflx2(:,:)
      character*80 name
      real (kind=4) ::  spv,add_offset,scale_factor
      real, allocatable   :: work_t(:)
      type(time_type) :: t1
      real :: mtime(12)
      integer :: start_y,end_y
      common / hflx_nao_reg/ start_y,end_y

       if (ntr == 1) then
c     read mean data from Barniers formulation
        allocate( qnet(nx,ny,12), qnec(nx,ny,12), airt(nx,ny,12))
        allocate( mask(nx,ny))
        call getunit(io,'data/ecmwf/ecmwf_BB.qnet','usr ieee')
        do n=1,12
         read (io); read (io) stamp, dpm, k, k, k, x,y,qnet(:,:,n)
         mtime(n)=dpm
        enddo
        close(io)
        call getunit(io,'data/ecmwf/ecmwf_BB.qnec','usr ieee')
        do n=1,12
         read (io); read (io) stamp, dpm, k, k, k, x,y,qnec(:,:,n)
        enddo
        close(io)
        call getunit(io,'data/ecmwf/ecmwf_BB.airt','usr ieee')
        do n=1,12
         read (io); read (io) stamp, dpm, k, k, k, x,y,airt(:,:,n)
        enddo
        close(io)
c
c     read NAO regression patterns
c
        iret=nf_open('data/NAO_reg/reg_hflx.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'hflx_reg',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,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), hflx2(nx2,ny2), tx(nx2),
     &           hflx(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, vid, 'missing_value', spv)
        add_offset=0.0
        iret= nf_get_att_real (ncid,vid ,'add_offset',add_offset)
        scale_factor=1.0
        iret= nf_get_att_real (ncid,vid ,'scale_factor',scale_factor)
        print*,' missing_value of hflx_reg = ',spv
        print*,' add_offset for   hflx_reg = ',add_offset
        print*,' scale_factor for hflx_reg = ',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
         hflx2=spv
         iret= nf_get_vara_real (ncid,vid ,
     &          (/1,1,n/),(/nx2-2,ny2-4,1/),hflx2(2:nx2-1,3:ny2-2))
         hflx2(nx2,:)=hflx2(2,:)
         hflx2(1,:)  =hflx2(nx2-1,:)
         hflx2(:,:)  =hflx2(:,ny2:1:-1)
         where( hflx2 == spv) hflx2=spval
         where( hflx2 /= spval) hflx2=(hflx2*scale_factor+add_offset)
         call tranlon (hflx2(1,1), nx2, nx2-2, ny2, x2, x,nx, tx)
         call ctf (hflx2(1,1), nx2, ny2, tx, y2, 
     &                  hflx(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
        enddo
        call ncclos (ncid, iret)
        print*,' done '
        deallocate(x2,y2,tx,work_t,hflx2)

c     read NAO index, multiply with regression pattern and add mean
c     change sign of regression patterns here !!!!
        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',vid)
        t1 =  set_date(start_y, 1, 1, 0, 0, 0)
        toffset=t1%days

        allocate(hflx2(nx,ny))
        do n=1,nt
         nn = mod(n-1,12)+1
         time(n)=mtime(nn)
         iret= nf_get_vara_real (ncid,vid ,n+(start_y-1865)*12,1,spv)
c         print*,' n=',n,nn,' nao(n)=',spv
         where (qnec(:,:,nn)/=0.) 
             hflx2(:,:)     = qnet(:,:,nn)+ hflx(:,:,nn)*spv*(-1)
             tr_clim(:,:,n) = airt(:,:,nn) + hflx2(:,:)/qnec(:,:,nn)
             tr_flux(:,:,n) = hflx2(:,:)
             tr_rest(:,:,n) = qnec(:,:,nn)*qfac
         elsewhere
             tr_clim(:,:,n) = spval
             tr_rest(:,:,n) = spval
             tr_flux(:,:,n) = spval
         end where
        enddo
        call ncclos (ncid, iret)
        deallocate(hflx2,hflx)
!     apply land mask
        call getunit(io,'data/ecmwf/ecmwf_BB.mask','usr ieee')
        read (io) ; read (io); read (io) mask; close(io)
        do n=1,nt
         where (mask==0.) tr_clim(:,:,n) = spval
         where (mask==0.) tr_rest(:,:,n) = spval
         where (mask==0.) tr_flux(:,:,n) = spval
        enddo
        deallocate( qnet, qnec, airt,mask)

       elseif (ntr == 2) then
c    surface salt flux for salinity
        tr_flux=0. ! unknown otherwise 
c      (flux is not used for restoring formulation)
        tr_rest=1./15.0*dzt(1)/86400.
c     read restoring salinity from previously interpolated data
c     which are also used for the initial conditions, sponge layers,  etc.
        call getunit(io,'ts.dta','usr ieee')
        do n=1,nt
         read(io) ; read(io) tr_clim(:,:,n) ! spval is already correct
         where( kmt == 0) tr_clim(:,:,n)=spval
         do m=3,number_tr
          read(io) ! tracer
         enddo
         do k=2,km; do m=1,number_tr
           read(io) ! rest of the levels
        enddo; enddo; enddo
        close(io)
        x=xt;y=yt
        time = (/31,28,31,30,31,30,31,31,30,31,30,31/)
        toffset=0.0
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_tr


#elif defined read_density_flx_nao_ideal
c
c--------------------------------------------------------
c    use the 3 year monthly mean ECMWF climatology plus
c    NAO related heat flux anomalies
c    (which is given by 
c    the monthly regression patterns of the 
c    NCEP net heat flux anomalies with the NAO-index). 
c    Same forcing was used in Eden+Willebrand 01, J.Clim.
c    See also prep_template_tau.F 
c--------------------------------------------------------
c
       subroutine read_dims_tr(nx,ny,nt,ntr)
       use prep_module
       implicit none
       integer nx,ny,nt,ntr
       character*32 stamp
       real :: dpm
       integer :: io

       if (ntr == 1 ) then
        print*,' reading heat flux boundary conditions from'
        print*,' Barniers ECMWF analysis, which were also used'
        print*,' in the DYNAMO project'
        print*,' plus NAO related heat flux anomalies from NCEP/NCAR'
        call getunit(io,'data/ecmwf/ecmwf_BB.qnet','usr ieee')
        read (io); read (io) stamp, dpm, nx, ny
        close(io)
        nt=12
       elseif (ntr == 2) then
        print*,' reading salt flux boundary condition from'
        print*,' monthly climatological SSS (file ts.dta) '
        print*,' and using a fixed restoring time scale '
#ifdef read_salt_fl_from_E5
        print*,' option read_salt_fl_from_E5 is not working here'
        stop
#endif
        nt=12; nx=imt; ny=jmt
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_dims_tr


       subroutine read_tr(tr_clim,tr_rest,tr_flux,
     &                    x,y,time,toffset,nx,ny,nt,ntr)
c
       ! returns inverse restoring time scale for tracers 
       ! number ntr in tr_rest
       ! and climatology to restore against in tr_clim
       ! nx,ny are horizontal dimensions of the data
       ! nt the time steps,
       ! x/y the longitude/latitude of the grid
       ! should mask data with spval
       ! return also net fluxes in Barniers heat flux formulation
       ! if available 
       ! Units: mom units for tr_clim and tr_rest
c       
       use prep_module
       use time_type_module
       implicit none
       integer nx,ny,nt,ntr,m,nn
       real tr_clim(nx,ny,nt) ,x(nx), y(ny)
       real tr_rest(nx,ny,nt) 
       real tr_flux(nx,ny,nt) ,time(nt),toffset
       character*32 stamp
       real, allocatable, dimension(:,:,:) :: qnet,qnec,airt
       real, allocatable, dimension(:,:)   :: mask
       real :: dpm
       integer  k,n,io,i,j
       real,parameter ::rhonull=1.026, Cpw=0.955, dimfac=0.2389e-4
       real,parameter ::qfac = dimfac/(rhonull*Cpw)
#include "netcdf.inc"
       integer :: iret,ncid,vid,dims(4),lonid,latid,corner(4),edges(4)
       character (len=80) :: dname
#ifdef netcdf_real4
      real (kind=4) :: v4(nx,ny) ,x4(nx),y4(ny),sp
#else
      real :: v4(nx,ny) ,x4(nx),y4(ny),sp
#endif
      integer :: init_t=1, lenw_t,nx2,ny2
      real (kind=4), allocatable :: x2(:),y2(:),hflx(:,:,:),tx(:),
     &               hflx2(:,:)
      character*80 name
      real (kind=4) ::  spv,add_offset,scale_factor
      real, allocatable   :: work_t(:)
      type(time_type) :: t1
      real :: mtime(12)

       if (ntr == 1) then
c     read mean data from Barniers formulation
        allocate( qnet(nx,ny,12), qnec(nx,ny,12), airt(nx,ny,12))
        allocate( mask(nx,ny))
        call getunit(io,'data/ecmwf/ecmwf_BB.qnet','usr ieee')
        do n=1,12
         read (io); read (io) stamp, dpm, k, k, k, x,y,qnet(:,:,n)
         mtime(n)=dpm
        enddo
        close(io)
        call getunit(io,'data/ecmwf/ecmwf_BB.qnec','usr ieee')
        do n=1,12
         read (io); read (io) stamp, dpm, k, k, k, x,y,qnec(:,:,n)
        enddo
        close(io)
        call getunit(io,'data/ecmwf/ecmwf_BB.airt','usr ieee')
        do n=1,12
         read (io); read (io) stamp, dpm, k, k, k, x,y,airt(:,:,n)
        enddo
        close(io)
c
c     read NAO regression patterns
c
        iret=nf_open('data/NAO_reg/reg_hflx.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'hflx_reg',vid)
        iret=NF_INQ_VARDIMID (NCID,vid,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), hflx2(nx2,ny2), tx(nx2),
     &           hflx(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, vid, 'missing_value', spv)
        add_offset=0.0
        iret= nf_get_att_real (ncid,vid ,'add_offset',add_offset)
        scale_factor=1.0
        iret= nf_get_att_real (ncid,vid ,'scale_factor',scale_factor)
        print*,' missing_value of hflx_reg = ',spv
        print*,' add_offset for   hflx_reg = ',add_offset
        print*,' scale_factor for hflx_reg = ',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
         hflx2=spv
         iret= nf_get_vara_real (ncid,vid ,
     &          (/1,1,n/),(/nx2-2,ny2-4,1/),hflx2(2:nx2-1,3:ny2-2))
         hflx2(nx2,:)=hflx2(2,:)
         hflx2(1,:)  =hflx2(nx2-1,:)
         hflx2(:,:)  =hflx2(:,ny2:1:-1)
         where( hflx2 == spv) hflx2=spval
         where( hflx2 /= spval) hflx2=(hflx2*scale_factor+add_offset)
         call tranlon (hflx2(1,1), nx2, nx2-2, ny2, x2, x,nx, tx)
         call ctf (hflx2(1,1), nx2, ny2, tx, y2, 
     &                  hflx(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
        enddo
        call ncclos (ncid, iret)
        print*,' done '
        deallocate(x2,y2,tx,work_t,hflx2)

c     read NAO index, multiply with regression pattern and add mean
c     change sign of regression patterns here !!!!
        t1 =  set_date(1999, 1, 1, 0, 0, 0)
        toffset=t1%days

        allocate(hflx2(nx,ny))
        do n=1,nt
         nn = mod(n-1,12)+1
         time(n)=mtime(nn)
         where (qnec(:,:,nn)/=0.) 
             hflx2(:,:)     = qnet(:,:,nn)+ hflx(:,:,nn)*3.0*(-1)
             tr_clim(:,:,n) = airt(:,:,nn) + hflx2(:,:)/qnec(:,:,nn)
             tr_flux(:,:,n) = hflx2(:,:)
             tr_rest(:,:,n) = qnec(:,:,nn)*qfac
         elsewhere
             tr_clim(:,:,n) = spval
             tr_rest(:,:,n) = spval
             tr_flux(:,:,n) = spval
         end where
        enddo
        deallocate(hflx2,hflx)
!     apply land mask
        call getunit(io,'data/ecmwf/ecmwf_BB.mask','usr ieee')
        read (io) ; read (io); read (io) mask; close(io)
        do n=1,nt
         where (mask==0.) tr_clim(:,:,n) = spval
         where (mask==0.) tr_rest(:,:,n) = spval
         where (mask==0.) tr_flux(:,:,n) = spval
        enddo
        deallocate( qnet, qnec, airt,mask)
       elseif (ntr == 2) then
c    surface salt flux for salinity
        tr_flux=0. ! unknown otherwise 
c      (flux is not used for restoring formulation)
        tr_rest=1./15.0*dzt(1)/86400.
c     read restoring salinity from previously interpolated data
c     which are also used for the initial conditions, sponge layers,  etc.
        call getunit(io,'ts.dta','usr ieee')
        do n=1,nt
         read(io) ; read(io) tr_clim(:,:,n) ! spval is already correct
         where( kmt == 0) tr_clim(:,:,n)=spval
         do m=3,number_tr
          read(io) ! tracer
         enddo
         do k=2,km; do m=1,number_tr
           read(io) ! rest of the levels
        enddo; enddo; enddo
        close(io)
        x=xt;y=yt
        time = (/31,28,31,30,31,30,31,31,30,31,30,31/)
        toffset=0.0
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_tr

#elif defined read_density_flx_omip
c
c--------------------------------------------------------
c    use the monthly mean OMIP climatology 
c--------------------------------------------------------
c
       subroutine read_dims_tr(nx,ny,nt,ntr)
       use prep_module
       implicit none
       integer nx,ny,nt,ntr
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4)

       if (ntr == 1 ) then
        print*,' reading heat flux boundary condition from'
        print*,' monthly mean climatology from OMIP '
        iret=nf_open('data/OMIP/mean_net_heatflux.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'net_heat',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)
       elseif (ntr == 2) then
        print*,' reading salt flux boundary condition from'
        print*,' monthly climatological SSS (file ts.dta) '
        print*,' and using a fixed restoring time scale '
#ifdef read_salt_fl_from_E5
        print*,' option read_salt_fl_from_E5 is not working here'
        stop
#endif
        nt=12; nx=imt; ny=jmt
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_dims_tr


       subroutine read_tr(tr_clim,tr_rest,tr_flux,
     &                    x,y,time,toffset,nx,ny,nt,ntr)
c
       ! returns inverse restoring time scale for tracers 
       ! number ntr in tr_rest
       ! and climatology to restore against in tr_clim
       ! nx,ny are horizontal dimensions of the data
       ! nt the time steps,
       ! x/y the longitude/latitude of the grid
       ! should mask data with spval
       ! return also net fluxes in Barniers heat flux formulation
       ! if available 
       ! Units: mom units for tr_clim and tr_rest
c       
       use prep_module
       implicit none
       integer nx,ny,nt,ntr
       real tr_clim(nx,ny,nt) ,x(nx), y(ny), tr_rest(nx,ny,nt) 
       real tr_flux(nx,ny,nt) ,time(nt),toffset
#include "netcdf.inc"
      integer :: iret,ncid,vid,dims(4),aid,n,m,k,io
      character (len=80) :: name
#ifdef netcdf_real4
      real (kind=4),allocatable :: v4(:,:,:),mask(:,:)
      real (kind=4) :: sp
#else
      real ,allocatable :: v4(:,:,:),mask(:,:)
      real  :: sp
#endif
       real,parameter ::rhonull=1.026, Cpw=0.955, dimfac=0.2389e-4
       real,parameter ::qfac = dimfac/(rhonull*Cpw)

       if (ntr == 1) then
        allocate( v4(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',vid)
        iret= nf_get_vara_real(ncid,vid,(/1,1/),(/nx,ny/),mask)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        call ncclos (ncid, iret)
c     read net heat flux
        iret=nf_open('data/OMIP/mean_net_heatflux.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'net_heat',vid)
        iret= nf_get_vara_real(ncid,vid,(/1,1,1/),(/nx,ny,nt/),v4)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         where ( v4(:,:,n) == sp ) v4(:,:,n)=spval
         where ( mask == 0.0  )    v4(:,:,n)=spval
         tr_flux(:,:,n)=v4(:,:,n)
        enddo
        ! read grid info
        iret=NF_INQ_VARDIMID (NCID,vid,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 dq/dt
        iret=nf_open('data/OMIP/mean_dq_dt.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'dq_dt',vid)
        iret= nf_get_vara_real(ncid,vid,(/1,1,1/),(/nx,ny,nt/),v4)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         where ( v4(:,:,n) == sp ) v4(:,:,n)=spval
         where ( mask == 0.0  )    v4(:,:,n)=spval
         tr_rest(:,:,n)=max(5.0,v4(:,:,n)*-1)
        enddo
        call ncclos (ncid, iret)
c     read SST
        iret=nf_open('data/OMIP/mean_sst.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'SST',vid)
        iret= nf_get_vara_real(ncid,vid,(/1,1,1/),(/nx,ny,nt/),v4)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
        do n=1,nt
         where ( v4(:,:,n) == sp ) v4(:,:,n)=spval
         where ( mask == 0.0  )    v4(:,:,n)=spval
         tr_clim(:,:,n)=v4(:,:,n)-273.15
        enddo
        call ncclos (ncid, iret)
c    construct Haney forcing, flip second dimension up/down
c    and convert units to MOM units
        do n=1,nt
         where (tr_rest(:,:,n)/=spval .and. tr_clim(:,:,n)/=spval 
     &           .and. tr_flux(:,:,n) /=spval  ) 
          tr_clim(:,:,n) = tr_clim(:,:,n) 
     &                      + tr_flux(:,:,n)/tr_rest(:,:,n)
          tr_rest(:,:,n) = tr_rest(:,:,n)*qfac
         elsewhere
          tr_clim(:,:,n) = spval
          tr_rest(:,:,n) = spval
          tr_flux(:,:,n) = spval
         end where
         tr_clim(:,:,n)=tr_clim(:,ny:1:-1,n)
         tr_flux(:,:,n)=tr_flux(:,ny:1:-1,n)
         tr_rest(:,:,n)=tr_rest(:,ny:1:-1,n)
        enddo
        y=y(ny:1:-1)
        deallocate(v4,mask)
        toffset=0.0
        time = (/31,28,31,30,31,30,31,31,30,31,30,31/)
       elseif (ntr == 2) then
c    surface salt flux for salinity
        tr_flux=0. ! unknown otherwise 
c      (flux is not used for restoring formulation)
        tr_rest=1./15.0*dzt(1)/86400.
c     read restoring salinity from previously interpolated data
c     which are also used for the initial conditions, sponge layers,  etc.
        call getunit(io,'ts.dta','usr ieee')
        do n=1,nt
         read(io) ; read(io) tr_clim(:,:,n) ! spval is already correct
         where( kmt == 0) tr_clim(:,:,n)=spval
         do m=3,number_tr
          read(io) ! tracer
         enddo
         do k=2,km; do m=1,number_tr
           read(io) ! rest of the levels
        enddo; enddo; enddo
        close(io)
        x=xt;y=yt
        time = (/31,28,31,30,31,30,31,31,30,31,30,31/)
        toffset=0.0
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_tr


#elif defined read_density_flx_T511
c
c--------------------------------------------------------
c    use the daily data from ECMWF/T511
c--------------------------------------------------------
c
       subroutine read_dims_tr(nx,ny,nt,ntr)
       use prep_module
       implicit none
       integer nx,ny,nt,ntr
#include "netcdf.inc"
      integer :: iret,ncid,varid,dims(4),nn
      integer :: start_y, fraction_y, nfraction_y,full_y
      common / flx_ncep/ start_y,fraction_y,nfraction_y,full_y
      character (len=80) :: file
      namelist /T511_time/ start_y,fraction_y,nfraction_y

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

       if (ntr == 1 ) then
        print*,' reading heat flux boundary condition from'
        print*,' daily data from ECMWF/T511 '
        print*,' year ',start_y,' in ',fraction_y,' fractions'
        print*,' current fraction ',nfraction_y
        write(file,'("data/ecmwf_T511/SSHF_",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,'SSHF',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
       elseif (ntr == 2) then
        print*,' no salt flux here'
        iret=nf_open('data/ecmwf_T511/SSHF_2001.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'SSHF',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=1
        call ncclos (ncid, iret)
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_dims_tr


       subroutine read_tr(tr_clim,tr_rest,tr_flux,
     &                    x,y,time,toffset,nx,ny,nt,ntr)
c
       ! returns inverse restoring time scale for tracers 
       ! number ntr in tr_rest
       ! and climatology to restore against in tr_clim
       ! nx,ny are horizontal dimensions of the data
       ! nt the time steps,
       ! x/y the longitude/latitude of the grid
       ! should mask data with spval
       ! return also net fluxes in Barniers heat flux formulation
       ! if available 
       ! Units: mom units for tr_clim and tr_rest
c       
       use prep_module
       use time_type_module
       implicit none
       integer nx,ny,nt,ntr
       real tr_clim(nx,ny,nt) ,x(nx), y(ny), tr_rest(nx,ny,nt) 
       real tr_flux(nx,ny,nt) ,time(nt),toffset
#include "netcdf.inc"
      integer :: iret,ncid,varid,xid,yid,i,j,n,n1,nn
      integer :: corner(4),edges(4),dims(4)
      character (len=80) :: name,file
#ifdef netcdf_real4
      real (kind=4),allocatable :: v4(:,:,:)
      real (kind=4) :: sp
      real (kind=4) :: add_offset,scale_factor
#else
      real ,allocatable :: v4(:,:,:)
      real  :: sp
      real  :: add_offset,scale_factor
#endif
        
       real :: mask(nx,ny),sp2,mask2(nx,ny),fxa
       type(time_type) :: t1,t2,t3
       real,parameter ::rhonull=1.026, Cpw=0.955, dimfac=0.2389e-4
       real,parameter ::qfac = dimfac/(rhonull*Cpw)
       integer :: start_y, fraction_y, nfraction_y,full_y
       common / flx_ncep/ start_y,fraction_y,nfraction_y,full_y

       
       if (ntr == 1) then

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/SSHF_",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,'SSHF',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/SSHF_",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,'SSHF',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)

        do n=1,nt; do j=1,ny; do i=1,nx
          if (mask(i,j) /= 1.0 ) then
           tr_flux(i,j,n) =v4(i,j,n)
          else
           tr_flux(i,j,n) = spval
          endif
        enddo; enddo; enddo
        call ncclos (ncid, iret)

        if (full_y /= 365 ) then
         iret=nf_open('data/ecmwf_T511/tair_leap.cdf',NF_nowrite,ncid)
        else
         iret=nf_open('data/ecmwf_T511/tair_2001.cdf',NF_nowrite,ncid)
        endif

        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'TAIR',varid)
        corner=(/1,1,n1,1/); edges = (/nx,ny,nt,1/)
        iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
        iret= nf_get_att_real (ncid,varid ,'missing_value',sp)
        do n=1,nt; do j=1,ny; do i=1,nx
          if (v4(i,j,n) /= sp ) then
           tr_clim(i,j,n) =v4(i,j,n)
          else
           tr_clim(i,j,n) = spval
          endif
        enddo; enddo; enddo
        call ncclos (ncid, iret)

        if (full_y /= 365 ) then
         iret=nf_open('data/ecmwf_T511/q2_leap.cdf',NF_nowrite,ncid)
        else
         iret=nf_open('data/ecmwf_T511/q2_2001.cdf',NF_nowrite,ncid)
        endif

        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'Q2',varid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        corner=(/1,1,n1,1/); edges = (/nx,ny,nt,1/)
        iret= nf_get_vara_real (ncid,varid ,corner,edges,v4)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret= nf_get_att_real (ncid,varid ,'missing_value',sp)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        do n=1,nt; do j=1,ny; do i=1,nx
          if (v4(i,j,n) /= sp ) then
           tr_rest(i,j,n) =v4(i,j,n)
          else
           tr_rest(i,j,n) = spval
          endif
        enddo; enddo; enddo
        call ncclos (ncid, iret)

        deallocate(v4)

        do n=1,nt; do j=1,ny; do i=1,nx
         if (tr_rest(i,j,n)/=spval .and. tr_clim(i,j,n)/=spval 
     &    .and. tr_flux(i,j,n) /=spval.and.tr_rest(i,j,n)/=0.0  ) then
          tr_clim(i,j,n) = tr_clim(i,j,n) 
     &                      + tr_flux(i,j,n)/tr_rest(i,j,n)
          tr_rest(i,j,n) = tr_rest(i,j,n)*qfac
         else
          tr_clim(i,j,n) = spval
          tr_rest(i,j,n) = spval
          tr_flux(i,j,n) = spval
         endif
        enddo; enddo; enddo

       elseif (ntr == 2) then
c no flux here
        tr_flux=0. 
        tr_rest=0.
        tr_clim=0.
        iret=nf_open('data/ecmwf_T511/land-sea-mask.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'lsm',varid)
        ! 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)
       else
        print*,' do not know SBC for tracer ',ntr
        call halt_stop(' in read_tr')
       endif
       end subroutine read_tr
#endif
