#include "options.inc"
c
c ------------------------------------------------------
c     templates to read climatology of temperature/salinity
c     and other tracers from external files used for initial
c     conditions, open boundaries and sponge layers.
c
c     Notes:
c      temperature should be potential temperature in
c      deg C, salinity in (psu-35)/1000.
c      Two subroutines are defined:
c       "read_dims_ic" and "read_ic" both called
c      in prep_ic (in prep_main.F).
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_ic_from_LEV94
c
c ------------------------------------------------------
c    read from WOA94 (Levitus 94)
c    (given as monthly means for T and S)
c ------------------------------------------------------
c
      subroutine read_dims_ic(nx,ny,nz,nt)
      implicit none
#include "netcdf.inc"
      integer nx,ny,nz,nt
      integer  ncid,iret,rid,ndims,latid,lonid
      integer dims(4)
      character (len=80) :: infile ='data/WOA94/lev94.cdf'
      print*,' reading T/S climatology from file ',
     &      infile(1:len_trim(infile))
      iret=nf_open(infile, NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'temp',rid)
      iret=NF_INQ_VARNDIMS(NCID, rid, ndims)  
      if (ndims/=4) then
       print*,' number of dimensions for variable not equal 4'
       call halt_stop(' in read_dims_ic (LEV94)')
      endif
      iret=NF_INQ_VARDIMID (NCID,rid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nz)
      iret=NF_INQ_DIMLEN (NCID, DIMs(4), nt)
      call ncclos (ncid, iret)
      end subroutine read_dims_ic


      subroutine read_ic(nx,ny,nz,ntr,nt,x,y,z,tr,spval)
      implicit none
      integer nx,ny,nz,nt,ntr
      real x(nx),y(ny),tr(nx,ny,nz,ntr),z(nz)
      real spval
#include "netcdf.inc"
      integer  ncid,iret,latid,lonid,tid,sid,zid,n
      integer  dims(4),corner(4),edges(4)
#ifdef netcdf_real4
      real (kind=4) :: v4(nx,ny,nz) ,x4(nx),y4(ny),z4(nz),sp
#else
      real :: v4(nx,ny,nz) ,x4(nx),y4(ny),z4(nz),sp
#endif
      integer :: indp,i,j,k,ii,jj,kk
      real :: fxa, fxb
      character*80 dname
      character (len=80) :: infile = 'data/WOA94/lev94.cdf'
      iret=nf_open(infile, NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'temp',tid)
      iret=nf_inq_varid(ncid,'salt',sid)
      iret=NF_INQ_VARDIMID (NCID,tid,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)
      dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(3),dname) 
      iret=nf_inq_varid(ncid,dname,zid)
      corner=1; edges(1)=nx; 
      iret= nf_get_vara_real (ncid,lonid ,corner,edges,x4); x=x4
      corner=1; edges(1)=ny; 
      iret= nf_get_vara_real (ncid,latid ,corner,edges,y4); y=y4
      corner=1; edges(1)=nz; 
      iret= nf_get_vara_real (ncid,zid ,corner,edges,z4); z=z4
      tr=0. ! setting all other tracers to zero
      corner=(/1,1,1,nt/); edges=(/nx,ny,nz,1/)
      iret= nf_get_vara_real (ncid,tid ,corner,edges,v4)
      iret= nf_get_att_real (ncid,tid ,'missing_value',sp)
      where( v4 == sp) v4 = -1000.0
      tr(:,:,:,1)=v4
      where( tr(:,:,:,1) <= -1000.0) tr(:,:,:,1) = spval
      iret= nf_get_vara_real (ncid,sid ,corner,edges,v4)
      iret= nf_get_att_real (ncid,sid ,'missing_value',sp)
c   apply special values, convert to cgs units
      where( v4 == sp) v4 = -1000.0
      tr(:,:,:,2)=v4
      where( tr(:,:,:,2) <= -1000.0 ) tr(:,:,:,2) = spval
      where( tr(:,:,:,2) /= spval) 
     &             tr(:,:,:,2) = (tr(:,:,:,2)-35.)/1000
      call ncclos (ncid, iret)

#ifdef  prep_for_CFC
c    initial condition is just zero
#endif
#ifdef  prep_for_LEDWELL
      if (ntr>=3) then
       print*,' constructing artificial initial conditions for tracer3'
       print*,' not yet implimented for LEV94'
       stop
      endif
#endif
#ifdef prep_for_OBC_tracer
      if (ntr>=3) then
       print*,' constructing initial conditions for tracer3'
       print*,' not yet implimented for LEV94'
       stop
      endif
#endif
#ifdef  prep_for_NPZD
      if (ntr>=3) then
       print*,' constructing initial conditions for tracer3'
       tr(:,:,:,4)=1e-2  ! phyto plankton
       tr(:,:,:,5)=1e-2  ! zoo
       tr(:,:,:,6)=1e-2  ! detritus
       print*,' reading nitrate from file data/WOA98/no3_annual.nc'
       iret=nf_open('data/WOA98/no3_annual.nc',NF_nowrite,ncid)
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=nf_inq_varid(ncid,'NO3',tid)
       v4=-99.9
       corner=(/1,1,1,1/); edges=(/nx-2,ny,nz,1/)
       iret= nf_get_vara_real (ncid,tid ,corner,edges,v4(2:nx-1,:,:))
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       tr(:,:,:,3)=v4
       where( tr(:,:,:,3) <= -99.9) tr(:,:,:,3) = spval
       call ncclos (ncid, iret)
c    limit nitrate to prevent spurious biological production
c    in restoring zones like the open boundary
       do k=1,nz
        do j=1,ny
          do i=1,nx
            if ( tr(i,j,k,3) <= 5.0 .and. tr(i,j,k,3)>=-1.0 ) then
                tr(i,j,k,3)=0.0
            endif
          enddo
        enddo
       enddo

       if (ntr >= 7) then
c       dissolved inorganic carbon
        print*,' reading DIC from data/GLODAP/glodap_dic.cdf'
        iret=nf_open('data/GLODAP/glodap_dic.cdf',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'ALK',tid)
        iret= nf_get_att_real (ncid,tid ,'missing_value',sp)
        v4=sp
        corner=(/1,1,1,1/); edges=(/nx-181+1,ny,nz,1/)
        iret= nf_get_vara_real (ncid,tid ,corner,edges,v4(181:nx,:,:))
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        corner=(/181,1,1,1/); edges=(/180,ny,nz,1/)
        iret= nf_get_vara_real (ncid,tid ,corner,edges,
     &                   v4(1:180,:,:))
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        where( v4 == sp) v4 = -1000.0
        tr(:,:,:,7)=v4
        where( tr(:,:,:,7) <= -1000.0) tr(:,:,:,7) = spval
        call ncclos (ncid, iret)
       endif

       if (ntr >= 8) then
c       total dissolved oxygen
        print*,' reading O2 from WOA98/oxygen_annual.nc'
        iret=nf_open('data/WOA98/oxygen_annual.nc',NF_nowrite,ncid)
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        iret=nf_inq_varid(ncid,'O2',tid)
        v4=-99.9
        corner=(/1,1,1,1/); edges=(/nx-2,ny,nz,1/)
        iret= nf_get_vara_real (ncid,tid ,corner,edges,v4(2:nx-1,:,:))
        if (iret/=0) then; print*,nf_strerror(iret);stop;endif
        tr(:,:,:,8)=v4
        where( tr(:,:,:,8) <= -99.9) tr(:,:,:,8) = spval
        call ncclos (ncid, iret)
       endif

       if (ntr >= 9) then
c        DOC, not yet implemented
       endif
      endif
#endif

      end subroutine read_ic

#elif defined read_ic_from_LEV98

c ------------------------------------------------------
c   read from a combined WOA94 + WOA98 dataset.
c   the high resolution (1/4 deg) seasonal mean data of 
c   WOA98 was combined with the monthly mean temporal
c   resolution of WOA94.
c ------------------------------------------------------

      subroutine read_dims_ic(nx,ny,nz,nt)
      implicit none
      integer nx,ny,nz,nt
      character*32 stamp
      real  dpm
      integer io
      print*,' reading T/S climatology from WOA9497'
      call getunit(io,'data/WOA9497/ann.temp','usr ieee')
      read (io); read (io) stamp, dpm, nx, ny
      close(io)
      nt=12; nz=33
      end subroutine read_dims_ic


      subroutine read_ic(nx,ny,nz,ntr,nt,x,y,z,tr,spval)
      implicit none
      integer nx,ny,nz,nt,ntr,indp
      real x(nx),y(ny),tr(nx,ny,nz,ntr),z(nz)
      real spval,dpm
      integer n,k,l,io,i,j,is,ie,js,je,ks,ke
      character(len=3) :: month(12)
      character(len=32) :: stamp
      data month/'jan','feb','mar','apr',
     &  'may','jun','jul','aug','sep','oct','nov','dec'/
      print*,' reading T/S climatology from WOA9497'
       tr=0. ! setting all other tracers to zero
       call getunit(io,'data/WOA9497/ann.temp','usr ieee')
       do k=1,19
        read(io); read(io)
       enddo
       do k=20,nz
        read(io); read(io) stamp, dpm, l,l,l,l,l,x,y,z,tr(:,:,k,1)
       enddo
       close(io)
       call getunit(io,'data/WOA9497/ann.salt','usr ieee')
       do k=1,19
        read(io); read(io)
       enddo
       do k=20,nz
        read(io); read(io) stamp, dpm, l,l,l,l,l,x,y,z,tr(:,:,k,2)
       enddo
       close(io)
       call getunit(io,'data/WOA9497/'//month(nt)//'.temp','usr ieee')
       do k=1,19
        read(io); read(io) stamp, dpm, l,l,l,l,l,x,y,z,tr(:,:,k,1)
       enddo
       close(io)
       call getunit(io,'data/WOA9497/'//month(nt)//'.salt','usr ieee')
       do k=1,19
        read(io); read(io) stamp, dpm, l,l,l,l,l,x,y,z,tr(:,:,k,2)
       enddo
       close(io)
c    apply special values, convert to cgs units
       do k=1,nz; do j=1,ny; do i=1,nx
         if (tr(i,j,k,2) > -99.9) tr(i,j,k,2)=(tr(i,j,k,2)-35.0)/1000.0
       enddo; enddo; enddo
       z=z/100
       do k=1,nz; do j=1,ny; do i=1,nx
         if (tr(i,j,k,1) <= -99.9) then
           tr(i,j,k,1)=spval; tr(i,j,k,2)=spval
         endif
       enddo; enddo; enddo

c    other tracers:

#ifdef  prep_for_ACO2
c    initial condition is just zero
#endif
#ifdef  prep_for_CFC
c    initial condition is just zero
#endif
#ifdef  prep_for_LEDWELL
      if (ntr>=3) then
       print*,' constructing artificial initial conditions for tracer3'
       print*,' not yet implimented for WOA9498'
      endif
#endif
#ifdef prep_for_OBC_tracer
      if (ntr>=3) then
       print*,' constructing initial conditions for tracer3'
       print*,' just zeros, except for '
       tr(:,:,:,3:ntr)=0.0
c    at y = 11.9N:12.2N from 60W to 56W    tr=1
       is=indp(360.0-60.0,x,nx)
       ie=indp(360.0-56.0,x,nx)
       js=indp(11.9,y,ny)
       je=indp(15.0,y,ny)
       ks=indp(1000.0,z,nz)
       print*,'is,ie,js,je,ks = ',is,ie,js,je,ks
       tr(is:ie,js:je,ks:nz,3)=1.0
       print*,' okay'
      endif
#endif
#ifdef  prep_for_NPZD
      if (ntr>=3) then
       print*,' reading Nitrate from WOA98/ann.no3'
       call getunit(io,'data/WOA98/ann.no3','usr ieee')
       do k=1,nz
        read(io); read(io) stamp, dpm, l,l,l,l,l,x,y,z,tr(:,:,k,3)
       enddo
       close(io)
       tr(:,:,:,4)=1e-2  ! phyto plankton
       tr(:,:,:,5)=1e-2  ! zoo
       tr(:,:,:,6)=1e-2  ! detritus

c    limit nitrate to prevent spurious biological production
c    in restoring zones like the open boundary
       do k=1,nz
        do j=1,ny
          do i=1,nx
            if ( tr(i,j,k,3) <= 5.0 .and. tr(i,j,k,3)>=-1.0 ) then
                tr(i,j,k,3)=0.0
            endif
          enddo
        enddo
       enddo

c    at southern boundary, set biology to zero
c    otherwise it will leak into model domain by restoring
       print*,' special !!!!'
       tr(:,200:230,:,4:6) = 0.0

       if (ntr >= 7) then
c       dissolved inorganic carbon
c        print*,' reading DIC from FRIIS/ann.pre_dic'
c        call getunit(io,'data/FRIIS/ann.pre_dic','usr ieee')
c        print*,' reading DIC from GOYET/ann.tco2'
c        call getunit(io,'data/GOYET/ann.tco2','usr ieee')
        print*,' reading DIC from GLODAP/ann.pre_dic'
        call getunit(io,'data/GLODAP/ann.pre_dic','usr ieee')
        do k=1,nz
         read(io); read(io) stamp, dpm, l,l,l,l,l,x,y,z,
     &     ((tr(i,j,k,7),i=1,nx),j=1,ny)
         do j=1,ny
          do i=1,nx
            if (tr(i,j,k,7)<0.0) then
                tr(i,j,k,7)=-99.9
            elseif (tr(i,j,k,7)> 5000.0) then
                tr(i,j,k,7)=-99.9
            endif
          enddo
         enddo
        enddo
        close(io)
       endif
       if (ntr >= 8) then
c       total dissolved oxygen
        print*,' reading O2 from WOA98/ann.oxygen'
        call getunit(io,'data/WOA98/ann.oxygen','usr ieee')
        do k=1,nz
         read(io); read(io) stamp, dpm, l,l,l,l,l,x,y,z,tr(:,:,k,8)
        enddo
        close(io)
       endif
#ifdef  prep_for_ACO2
       if (ntr >= 9) then
c       dissolved inorganic carbon from fixed pco2_air
        print*,' reading CO2_C from GLODAP/ann.pre_dic'
        call getunit(io,'data/GLODAP/ann.pre_dic','usr ieee')
c        print*,' reading DIC from FRIIS/ann.pre_dic'
c        call getunit(io,'data/FRIIS/ann.pre_dic','usr ieee')
c        print*,' reading DIC from GOYET/ann.tco2'
c        call getunit(io,'data/GOYET/ann.tco2','usr ieee')
        do k=1,nz
         read(io); read(io) stamp, dpm, l,l,l,l,l,x,y,z,
     &     ((tr(i,j,k,9),i=1,nx),j=1,ny)
         do j=1,ny
          do i=1,nx
            if (tr(i,j,k,9)<0.0) then
                tr(i,j,k,9)=-99.9
            elseif (tr(i,j,k,9)> 5000.0) then
                tr(i,j,k,9)=-99.9
            endif
          enddo
         enddo
        enddo
        close(io)
       endif
#else
       if (ntr >= 9) then
c      dissolved organic carbon
         tr(:,:,:,9)=0.0
       endif
#endif
c      convert to cgs units
       z=z/100
c      apply special values
       do n=3,ntr; do k=1,nz; do j=1,ny; do i=1,nx
         if (tr(i,j,k,n) <= -99.9) tr(i,j,k,n)=spval
       enddo; enddo; enddo; enddo
      endif
#endif

      end subroutine read_ic




#elif defined read_ic_from_SAC

c ------------------------------------------------------
c   read from the SAC (Woce) climatology 
c    (annual mean only)
c    (note that the data comes as in situ temperature)
c ------------------------------------------------------

      subroutine read_dims_ic(nx,ny,nz,nt)
      implicit none
#include "netcdf.inc"
      integer nx,ny,nz,nt
      integer  ncid,iret,rid,ndims,latid,lonid
      integer dims(4)
      character (len=80) :: infile ='data/SAC/sac_wtso.cdf'
      print*,' reading T/S climatology from file ',
     &      infile(1:len_trim(infile))
      nt=12
      iret=nf_open(infile,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'TEMP',rid)
      iret=NF_INQ_VARNDIMS(NCID, rid, ndims)  
      if (ndims/=3) then
       print*,' number of dimensions for variable not equal 3'
       print*,ndims
       call halt_stop(' in read_dims_ic (SAC)')
      endif
      iret=NF_INQ_VARDIMID (NCID,rid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), ndims); nx=ndims
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ndims); ny=ndims
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), ndims); nz=ndims
      call ncclos (ncid, iret)
      end subroutine read_dims_ic


      subroutine read_ic(nx,ny,nz,ntr,nt,x,y,z,tr,spval)
      implicit none
      integer nx,ny,nz,nt,ntr
      real x(nx),y(ny),tr(nx,ny,nz,ntr),z(nz)
      real spval, fxa
#include "netcdf.inc"
      integer  ncid,iret,latid,lonid,tid,sid,zid,n
      integer  dims(4),corner(4),edges(4),i,j,k
#ifdef netcdf_real4
      real (kind=4) :: v4(nx,ny,nz) ,x4(nx),y4(ny),z4(nz),sp
#else
      real :: v4(nx,ny,nz) ,x4(nx),y4(ny),z4(nz),sp
#endif
      character*80 dname
      character (len=80) :: infile = 'data/SAC/sac_wtso.cdf'
      iret=nf_open(infile, NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'TEMP',tid)
      iret=nf_inq_varid(ncid,'SAL',sid)
      iret=NF_INQ_VARDIMID (NCID,tid,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)
      dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(3),dname) 
      iret=nf_inq_varid(ncid,dname,zid)
      corner=1; edges(1)=nx; 
      iret= nf_get_vara_real (ncid,lonid ,corner,edges,x4); x=x4
      corner=1; edges(1)=ny; 
      iret= nf_get_vara_real (ncid,latid ,corner,edges,y4); y=y4
      corner=1; edges(1)=nz; 
      iret= nf_get_vara_real (ncid,zid ,corner,edges,z4); z=z4
      tr=0. ! setting all other tracers to zero
      corner=1; edges(1)=nx; edges(2)=ny; edges(3)=nz
      iret= nf_get_vara_real (ncid,tid ,corner,edges,v4)
      iret= nf_get_att_real (ncid,tid ,'missing_value',sp)
      where( v4 == sp) v4 = -1000.0
      tr(:,:,:,1)=v4
      where( tr(:,:,:,1) <= -1000.0) tr(:,:,:,1) = spval
      iret= nf_get_vara_real (ncid,sid ,corner,edges,v4)
      iret= nf_get_att_real (ncid,sid ,'missing_value',sp)
      where( v4 == sp) v4 = -1000.0
      tr(:,:,:,2)=v4
      where( tr(:,:,:,2) <= -1000.0 ) tr(:,:,:,2) = spval
c   convert in situ to potential temperature
      do i=1,nx; do j=1,ny; do k=1,nz
          if (tr(i,j,k,1) /=spval ) then
           call potem(tr(i,j,k,1),tr(i,j,k,2),z(k),fxa)
           tr(i,j,k,1)=fxa
          endif
      enddo; enddo; enddo;
      where( tr(:,:,:,2) /= spval) 
     &             tr(:,:,:,2) = (tr(:,:,:,2)-35.)/1000
      call ncclos (ncid, iret)
      end subroutine read_ic


      subroutine potem (t, s, p, theta)
c     this subroutine calculates potential temperature as a function
c     of in-situ temperature, salinity, and pressure.
c     input [units]:
c       in-situ temperature (t): [degrees centigrade]
c       salinity (s): [per mil] psu ??
c       pressure (p): [decibars, approx. as meters of depth]
c     output [units]:
c       potential temperature (theta): [degrees centigrade]
c     references:
c        based on Fofonoff and Froese (1958) as shown in ...
c        Fofonoff, N., The Sea: Vol 1, (ed. M. Hill). Interscience,
c          New York, 1962, page 17, table iv.
      implicit real (a-h,o-z)
      b1    = -1.60d-5*p
      b2    = 1.014d-5*p*t
      t2    = t*t
      t3    = t2*t
      b3    = -1.27d-7*p*t2
      b4    = 2.7d-9*p*t3
      b5    = 1.322d-6*p*s
      b6    = -2.62d-8*p*s*t
      s2    = s*s
      p2    = p*p
      b7    = 4.1d-9*p*s2
      b8    = 9.14d-9*p2
      b9    = -2.77d-10*p2*t
      b10   = 9.5d-13*p2*t2
      b11   = -1.557d-13*p2*p
      potmp = b1+b2+b3+b4+b5+b6+b7+b8+b9+b10+b11
      theta = t-potmp
      end subroutine potem


#elif defined read_ic_from_E5

c ------------------------------------------------------
c   read from model climatology
c ------------------------------------------------------

      subroutine read_dims_ic(nx,ny,nz,nt)
      implicit none
#include "netcdf.inc"
      integer nx,ny,nz,nt
      integer  ncid,iret,rid,ndims,latid,lonid
      integer dims(4)
      character (len=80) :: infile = 
     &  '/users/model/ceden/dov3/N_atl/E5_bbl/ts_mmean_ufo.cdf'
      print*,' reading T/S climatology from file ',
     &      infile(1:len_trim(infile))
      iret=nf_open(infile, NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'temp',rid)
      iret=NF_INQ_VARNDIMS(NCID, rid, ndims)  
      if (ndims/=4) then
       print*,' number of dimensions for variable not equal 4'
       call halt_stop(' in read_dims_ic (E5)')
      endif
      iret=NF_INQ_VARDIMID (NCID,rid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nz)
      iret=NF_INQ_DIMLEN (NCID, DIMs(4), nt)
      call ncclos (ncid, iret)
      end subroutine read_dims_ic



      subroutine read_ic(nx,ny,nz,ntr,nt,x,y,z,tr,spval)
      implicit none
      integer nx,ny,nz,nt,ntr
      real x(nx),y(ny),tr(nx,ny,nz,ntr),z(nz)
      real spval
#include "netcdf.inc"
      integer  ncid,iret,latid,lonid,tid,sid,zid,n
      integer  dims(4),corner(4),edges(4)
#ifdef netcdf_real4
      real (kind=4) :: v4(nx,ny,nz) ,x4(nx),y4(ny),z4(nz),sp
#else
      real :: v4(nx,ny,nz) ,x4(nx),y4(ny),z4(nz),sp
#endif
      character*80 dname
      character (len=80) :: infile = 
     &  '/users/model/ceden/dov3/N_atl/E5_bbl/ts_mmean_ufo.cdf'
      iret=nf_open(infile, NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,'temp',tid)
      iret=nf_inq_varid(ncid,'salt',sid)
      iret=NF_INQ_VARDIMID (NCID,tid,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)
      dname='';iret=NF_INQ_DIMNAME (NCID, DIMs(3),dname) 
      iret=nf_inq_varid(ncid,dname,zid)
      corner=1; edges(1)=nx; 
      iret= nf_get_vara_real (ncid,lonid ,corner,edges,x4); x=x4
      corner=1; edges(1)=ny; 
      iret= nf_get_vara_real (ncid,latid ,corner,edges,y4); y=y4
      corner=1; edges(1)=nz; 
      iret= nf_get_vara_real (ncid,zid ,corner,edges,z4); z=z4
      tr=0. ! setting all other tracers to zero
      corner=(/1,1,1,nt/); edges=(/nx,ny,nz,1/)
      iret= nf_get_vara_real (ncid,tid ,corner,edges,v4)
      iret= nf_get_att_real (ncid,tid ,'missing_value',sp)
      where( v4 == sp) v4 = -1000.0
      tr(:,:,:,1)=v4
      where( tr(:,:,:,1) <= -1000.0) tr(:,:,:,1) = spval
      iret= nf_get_vara_real (ncid,sid ,corner,edges,v4)
      iret= nf_get_att_real (ncid,sid ,'missing_value',sp)
      where( v4 == sp) v4 = -1000.0
      tr(:,:,:,2)=v4
      where( tr(:,:,:,2) <= -1000.0 ) tr(:,:,:,2) = spval
      where( tr(:,:,:,2) /= spval) 
     &             tr(:,:,:,2) = (tr(:,:,:,2)-35.)/1000
      call ncclos (ncid, iret)
      end subroutine read_ic
#endif
