module interp 
    use const
    use comm
    use var
!TS    use io
implicit real*8 (a-h,o-z)
#   include "pnetcdf.inc" !dog
contains

subroutine inmets
    ! ======================================================================
    ! extract regional subset of etopo5 bathymetry
    ! user specified lat-long window and resolution
    ! by david dietrich, may 1, 1995
    !
    ! reads one latitude at a time in a simple do loop.
    ! reads full etopo5 data sets into memory, then extracts desired subset
    ! simple and versatile.
    ! most workstations have the required storage capacity.
    ! etopo5 has 360 x 12 integer*2 data elements per (latitudinal) record.
    !   this gives 8640 bytes or 2160 sgi (or vax) "words". thus, under sun
    !   fortran, recl=8640, while under sgi or vax fortran recl=2160.
    ! ======================================================================
    !nres is the etopo resolution =5 (5min)
    parameter(nres=2)
    parameter(nx=21600/nres,ny=10800/nres,nxpd=nx/360,nypd=ny/180)
    integer*2 meters,in

    ! ========================================================
    ! tlev, slev is levitus climatology data with landfills
    ! t,s is tlev, slev interpolated to model horizontal grid
    ! t1, s1 is t,s vertically interpolated to model z-levels
    ! ========================================================
    dimension rmets(i0,j0),meters(nx,ny),in(i0,j0)
    dimension depth(nx,ny)
    dimension istart(2),icount(2),dt(nx,ny)
    ! read nc data, which is single precision
    real,dimension(nx,ny) :: dep_tmp
    ! ======================================================================
    ! i0= number of grid points in the x-direction including ghost zones
    ! j0= number of grid points in the y-direction including ghost zones
    ! westdeg is western most longitude degrees (west cell face of i=2 zone)
    ! longitude coordinate increases eastward
    ! 0 < westdeg < 360
    ! longitudinal increment in minutes,dxmnut, is read from yzgrid
    ! latitudes are read from yzgrid
    ! ======================================================================

!    open(9,file='data/westdeg'//grd)
!    rewind 9
!    read(9,'(f9.4)') westdeg
!    close(9)
    call rdnc_westdeg !dog
    ! ================
    ! input data files
    ! ================
!    open(10,file='data/yzgrid'//grd,form='unformatted')

    ! new, single record integer*2 data
    ! no record length bullshit!
    if (myid .eq. 0)  then
	if(nres==1)then
	    istatus=nf_open('../worldata/etopo/etopo1/etopo1_bed_c_gmt4.grd',nf_nowrite,ncid)
	    print*,'istatus1',istatus,ncid,nx,ny
	    elseif(nres==2) then
	    istatus=nf_open('../worldata/etopo/etopo2/etopo2v2g_f4.nc',nf_nowrite,ncid)
	    print*,'istatus2',istatus,ncid,nx,ny,NF_STRERROR(istatus)
	    elseif(nres==5) then
	    open(12,file='../worldata/etopo/etopo5/itopo5',form='unformatted')
	endif
    endif
!    rewind 10
    rewind 12
    rewind 35

    ! =================
    ! output data files
    ! =================
!    open(8,file='data/depth'//grd,form='unformatted')
!    open(14,file='data/inview'//grd)
!    rewind 8
!    rewind 14

    ! ==============================================
    ! read metrics data file yzgrid (from metgen.f)
    ! ==============================================
    call rdnc_yzgrid(j,k,dxmnut)  !dog
 !   read(10) j,k
!    if (j0.eq.j.and.k0.eq.k) go to 9
    if (j0.ne.j.and.k0.ne.k) then
    write(*,6)
    6    format('stop due to resolution mismatch. set j0,k0 to values from metgen.f.')
    stop
    endif
!    9    read(10) tlz,z,odz,odzw,dxmnut,xdeg,y,yv,yvdeg,ydeg,cs,csv,ocs,ocsv,dx,dxv,odx,odxv,dy,dyv,ody,odyv
        
    ! new single record input
    if (myid .eq. 0) then
	if(nres.le.2)then
	    istart(1) = 1
	    istart(2) = 1
	    icount(1) = nx
	    icount(2) = ny
	    istatus=nf_get_vara_real(ncid,3,istart,icount,dep_tmp)
            depth = dep_tmp
	    do j=1,ny
		do i=1,nx/2 
		    depth(i,j)=depth(i,j)+depth(i+nx/2,j)
		    depth(i+nx/2,j)=depth(i,j)-depth(i+nx/2,j)
		    depth(i,j)=depth(i,j)-depth(i+nx/2,j)
		enddo
	    enddo
	    print*,'icount,',icount
	    do j=1,ny
		do i=1,nx
		    dt(i,j)=depth(i,ny-j+1)
		enddo
	    enddo
	    do j=1,ny
		do i=1,nx
		    depth(i,j)=dt(i,j)
		enddo
	    enddo
	    !ts   will add this later
	else if(nres.eq.5) then
	    read(12) meters
	    do j=1,ny
		do i=1,nx
		    depth(i,j)=meters(i,j)
		enddo
	    enddo
	endif
    endif

    call mpi_bcast(depth,nx*ny,mpi_real8,0,mpi_comm_world,ierr)

    ! ==============================================================
    ! determine depths at model grid points and store in rmets array
    ! ==============================================================
    ! first element is at 0.0 degrees east (or 360 degrees east).
    ! 
    if(nres.eq.1)then
	xinc=dxmnut/1.
    else if(nres.eq.2)then
	xinc=dxmnut/2.
    else if(nres.eq.5)then
	xinc=dxmnut/5.
    endif

    if(myid.eq.0) print*,'xinc',xinc

    ! xin! is the grid distance of etop2. therefore, for 15' grid size
    ! each grid has 15/2 = 7.5 grid point in etop2 data set.
    do 100 jj=1,j0
	! ylat is measured from north pole.  first point is at pole.
	! xlon is measured units 1/12 deg east, from grenwich, with first point
	! at 1/12 deg east (meters(1,j) is at 1/12 deg east; meters(i,1) is
	! at north pole, or ydeg=90)
	ylat=nypd*(90-ydeg(jj))+1
	! average depths in lat-long window
	j=ylat
	jp=j-1
	xlon=nxpd*westdeg-.5*xinc
	!      print*,'xinc',xinc,xlon,ydeg(524)
      
	do 100 ii=2,i1
	    xlon=xlon+xinc
	    i=xlon
	    tmp1=xlon-i
	    tmp2=1.-tmp1
	    ip=mod(i+1,nx)+1
	    i=mod(i,nx)+1
	    !choke point depths are critical in modeling straits. when they are not
	    ! adequately resolved they are always too shallow when evaluating them from
	    ! world data base. thus, at a given control volume location, it is
	    ! physically better to assign the deepest adjacent point in the control
	    ! volume region than to use the actual local value.

	    ! special for 1/4 deg global grid in critical gom and gibraltar regions
	    if (xdeg(ii).lt.269.625.or.xdeg(ii).gt.359.125) go to 80
	    ! avoid cape hatteras and northern north atlanti! region
	    if (xdeg(ii).lt.299.625.and.ydeg(jj).gt.29.91505) go to 80
	    ! this may make some single-point islands wet, but that's ok!
	    ! note: negative meters values mean positive ocean depths


	    rmets(ii,jj)=min(depth(i,j),depth(ip,j),depth(i,jp),depth(ip,jp))

	    rmets(ii,jj)=-rmets(ii,jj)
	    rmets(ii,jj)=max(0.,rmets(ii,jj))
	    go to 100

80  dsw=depth(i,jp)
    dse=depth(ip,jp)
    dnw=depth(i,j)
    dne=depth(ip,j)
    dn=tmp1*dne+tmp2*dnw
    ds=tmp1*dse+tmp2*dsw
    rmets(ii,jj)=(ylat-j)*ds+(1-ylat+j)*dn
    rmets(ii,jj)=-rmets(ii,jj)
    rmets(ii,jj)=max(0.,rmets(ii,jj))
100 continue

!    open(333,file='data/rmets'//grd,form='unformatted')
!    write(333) rmets
!    close(333)

!    write(8) rmets
    call wr8nc_depth(rmets)  !dog

    m=0
    n=0
    tmp=0.
    temp=.01*tlz
    do j=2,j1
	do i=2,i1
	    m=m+in(i,j)
	    tmp=max(tmp,rmets(i,j))
	    if (rmets(i,j).gt.temp) n=n+1
	enddo  
    enddo  
    write(*,103) m,n,temp,tmp
    103  format(i8,' total wet points,',i6,' deeper than ',f5.0,' m, deepest= ',f6.0,' m')

end subroutine

subroutine indata
! ======================================================================
! read bathymetry file from inmets.f
! read salinity, temperature and wind data sets
! interpolate data to model grid
! ======================================================================
! lscowf: .true.  scow wind climatological forcing
! lscowf: .false. hr   wind climatological forcing
! lclim: 1 levitus 94 climatology
! lclim: 2 phc3 climatology
    parameter(lscowf=.true.,lclim=2)
    parameter(nxscow=1440,nyscow=560,nxhel=360,nyhel=180)
    integer*2 nzdat
!    real*4 slon,slat,ilon,ilat,tlevs,slevs

! real*4 tmp variables for file reading
    real*4,dimension(360,180,33) :: swoa_tmp,twoa_tmp
    real*4,dimension(360,180,33) :: tlevs_tmp,slevs_tmp
    real*4,dimension(nxhel,nyhel) :: tx_tmp,ty_tmp
    real*4,dimension(nxscow,nyscow) :: sctx_tmp,scty_tmp
! ========================================================
! tlev, slev is levitus climatology data with landfills
! t,s is tlev, slev interpolated to model horizontal grid
! t1_data, s1_data is t_data,s_data vertically interpolated to model z-levels
! ========================================================
      dimension tlev(360,180,33),slev(360,180,33),rho_data(360,180,33), &
      t_data(i0,j0,33),s_data(i0,j0,33),t1_data(i0,j0,k1),s1_data(i0,j0,k1),taux_data(i2,j2), &
      tauy_data(i2,j2),tauxs(i2,j2),tauys(i2,j2),nzdat(33)
    
      dimension tlevs(360,180,33),slevs(360,180,33)

      dimension twoa(360,180,33),swoa(360,180,33)
      dimension twoal(0:361,0:181),swoal(0:361,0:181)
      dimension tx_data(nxhel,nyhel),ty_data(nxhel,nyhel)
      dimension scowtx_data(nxscow,nyscow),scowty_data(nxscow,nyscow),istart(2),icount(2)
      
      character*59 txc(12)
      character*59 tyc(12)
      character*59 season(4)
      character*60 twoaf(12),swoaf(12),seasonswoa(4),seasontwoa(4)
      character*80 scowtaux,scowtauy

      integer*2 maxj(i2,j2)
      dimension scr_data(i0,j0,2)
      real*8,dimension(i0,j0,k1,12) :: clidata_t1,clidata_s1
      real*8,dimension(i2,j2,12) :: cli_tauxs,cli_tauys
      real*8,dimension(i2,j2,12) :: clidata_tauy,clidata_taux

      data txc/				    &
       '../worldata/hellerman/tx-jan.dat',  &
       '../worldata/hellerman/tx-feb.dat',  &
       '../worldata/hellerman/tx-mar.dat',  &
       '../worldata/hellerman/tx-apr.dat',  &
       '../worldata/hellerman/tx-may.dat',  &
       '../worldata/hellerman/tx-jun.dat',  &
       '../worldata/hellerman/tx-jul.dat',  &
       '../worldata/hellerman/tx-aug.dat',  &
       '../worldata/hellerman/tx-sep.dat',  &
       '../worldata/hellerman/tx-oct.dat',  &
       '../worldata/hellerman/tx-nov.dat',  &
       '../worldata/hellerman/tx-dec.dat'/
      data tyc/				    &
       '../worldata/hellerman/ty-jan.dat',  &
       '../worldata/hellerman/ty-feb.dat',  &
       '../worldata/hellerman/ty-mar.dat',  &
       '../worldata/hellerman/ty-apr.dat',  &
       '../worldata/hellerman/ty-may.dat',  &
       '../worldata/hellerman/ty-jun.dat',  &
       '../worldata/hellerman/ty-jul.dat',  &
       '../worldata/hellerman/ty-aug.dat',  &
       '../worldata/hellerman/ty-sep.dat',  &
       '../worldata/hellerman/ty-oct.dat',  &
       '../worldata/hellerman/ty-nov.dat',  &
       '../worldata/hellerman/ty-dec.dat'/
      data season/				&
       '../worldata/levitus/modelev_win.dat',   &
       '../worldata/levitus/modelev_spr.dat',   &
       '../worldata/levitus/modelev_sum.dat',   &
       '../worldata/levitus/modelev_aut.dat'/
      data swoaf/			    &
       '../worldata/phc3/salt01_p3',        &
       '../worldata/phc3/salt02_p3',        &
       '../worldata/phc3/salt03_p3',        &
       '../worldata/phc3/salt04_p3',        &
       '../worldata/phc3/salt05_p3',        &
       '../worldata/phc3/salt06_p3',        &
       '../worldata/phc3/salt07_p3',        &
       '../worldata/phc3/salt08_p3',        &
       '../worldata/phc3/salt09_p3',        &
       '../worldata/phc3/salt10_p3',        &
       '../worldata/phc3/salt11_p3',        &
       '../worldata/phc3/salt12_p3'/
      data twoaf/			    &
       '../worldata/phc3/temp01_p3',        &
       '../worldata/phc3/temp02_p3',        &
       '../worldata/phc3/temp03_p3',        &
       '../worldata/phc3/temp04_p3',        &
       '../worldata/phc3/temp05_p3',        &
       '../worldata/phc3/temp06_p3',        &
       '../worldata/phc3/temp07_p3',        &
       '../worldata/phc3/temp08_p3',        &
       '../worldata/phc3/temp09_p3',        &
       '../worldata/phc3/temp10_p3',        &
       '../worldata/phc3/temp11_p3',        &
       '../worldata/phc3/temp12_p3'/

      data seasonswoa/			    & 		
       '../worldata/phc3/salt13_p3',        &
       '../worldata/phc3/salt15_p3',        &
       '../worldata/phc3/salt15_p3',        &
       '../worldata/phc3/salt13_p3'/
      data seasontwoa/			    &
       '../worldata/phc3/temp13_p3',        &
       '../worldata/phc3/temp15_p3',        &
       '../worldata/phc3/temp15_p3',        &
       '../worldata/phc3/temp13_p3'/

      data scowtauy/			    & 
       '../worldata/scow/ws_v.nc'/          
      data scowtaux/                        &
       '../worldata/scow/ws_u.nc'/          

! -------------------------------------
! standard levitus data depths (m)
! -------------------------------------
    data nzdat/   0,  10,  20,  30,  50,  75, 100, 125, 150, 200, 250,  &
      300, 400, 500, 600, 700, 800, 900,1000,1100,1200,1300,1400,1500, &
      1750,2000,2500,3000,3500,4000,4500,5000,5500/
!lat: /lat (degree_north) ordered (89.5s) to (89.5n) by 1.0 n= 180 pts :grid
!lon: /lon (degree_east) periodi! (0.5e) to (0.5w) by 1.0 n= 360 pts :grid
! ======================================================================
! i0= number of grid points in the x-direction including ghost zones
! j0= number of grid points in the y-direction including ghost zones
! westdeg is western most longitude degrees (west cell face of i=2 zone)
! longitude coordinate increases eastward
! 0 < westdeg < 360
! longitudinal increment in minutes,dxmnut, is read from yzgrid
! latitudes are read from yzgrid
! ======================================================================
! ================
! input data files
! ================
!      open(10,file='data/yzgrid'//grd,form='unformatted')
!      open(14,file='data/invu'//grd)

!      rewind 10
!      rewind 14

! =================
! output data files
! =================
      if (myid .eq. 0) open(41,file='data/fixedlev',form='unformatted')

! annual cycle levitus climatology on model grid
!      open(70,file='data/annualevitus'//grd,form='unformatted')

! s,t at surface and lateral boundaries
!      open(72,file='data/winds'//grd,form='unformatted')

      rewind 41
!      rewind 70
!      rewind 72

! ==============================================
! read metrics data file yzgrid (from metrics.f)
! ==============================================
      call rdnc_yzgrid(j,k,dxmnut) !dog
!      read(10) j,k
      if (j0.ne.j.and.k0.ne.k) then
      write(*,6)
 6    format('stop due to resolution mismatch. set j0,k0 to values from metgen.f.')
      stop
      endif
! 9    read(10) tlz,z,odz,odzw,dxmnut,xdeg,y,yv,yvdeg,ydeg,cs,csv,ocs,ocsv,dx,dxv,odx,odxv,dy,dyv,ody,odyv


!     scow wind forcing 
      if (myid.eq.0.and.lscowf.eq..true.)  then
      ist=nf_open(scowtaux,nf_nowrite,nctx)
      ist=nf_open(scowtauy,nf_nowrite,ncty)
      istart(1)=1
      istart(2)=1
      icount(1)=nxscow
      icount(2)=nyscow
      endif
! ==================================================
! prepare monthly winds and seasonal s,t climatology
! ==================================================
      do m=1,12
      if (myid .eq. 0) then

      rewind 41
      write(*,111) m
 111  format(i3)
! hellerman
      open(15,file=txc(m))
      open(16,file=tyc(m))
      rewind 15
      rewind 16
      
      if (mod(m,3).eq.1) then
      nm=(m+2)/3
! ===================================================
! read modelfied (with landfills) levitus climatology
! and fix negative stratification areas
! ===================================================

      open(36,file=season(nm),form='unformatted')
      rewind 36

      if(lclim.ne.1)then
      open(37,file=seasonswoa(nm))
      open(38,file=seasontwoa(nm))
      rewind 37
      rewind 38
      endif
! ===================================================
! read modelfied (with landfills) levitus climatology
! and fix negative stratification areas
! ===================================================
!      read(40)
!     1 ((tlev(i,j,1),i=1,360),j=1,180),((slev(i,j,1),i=1,360),j=1,180)

      print*,'read seasonal data'
      do k=1,33
       read(36) ((tlevs_tmp(i,j,k),i=1,360),j=1,180),((slevs_tmp(i,j,k),i=1,360),j=1,180) 
      enddo
      tlevs = tlevs_tmp     
      slevs = slevs_tmp 

      if(lclim.ne.1)then
      do k=1,33
      do j=1,180
      do i=1,36
       read(37,'(10f8.4)') (swoa_tmp(ii,j,k),ii=(i-1)*10+1,i*10)
       read(38,'(10f8.4)') (twoa_tmp(ii,j,k),ii=(i-1)*10+1,i*10)
      enddo
      enddo
      enddo
      swoa = swoa_tmp
      twoa = twoa_tmp
      endif

      endif

!      read(40)
!     1 ((tlev(i,j,1),i=1,360),j=1,180),((slev(i,j,1),i=1,360),j=1,180)
      if(lclim.ne.1)then
      open(39,file=swoaf(m))
      open(40,file=twoaf(m))
      rewind 39
      rewind 40
! ===================================================
! read modelfied (with landfills) levitus climatology
! and fix negative stratification areas
! ===================================================
      do 124 k=1,24
      do 124 j=1,180
      do 124 i=1,36
      read(40,'(10f8.4)') (twoa_tmp(ii,j,k),ii=(i-1)*10+1,i*10)
 124  read(39,'(10f8.4)') (swoa_tmp(ii,j,k),ii=(i-1)*10+1,i*10)
      swoa = swoa_tmp
      twoa = twoa_tmp
      endif
!      print*,'month twoa',(twoa(180,j,1),j=1,180)

 127  continue

      if(lclim.eq.1)then
      do k=1,33
      do j=1,180 
      do i=1,360
      swoa(i,j,k)=slevs(i,j,k)
      twoa(i,j,k)=tlevs(i,j,k)
      enddo
      enddo
      enddo
      else
      do 128 k=1,33
      do 128 j=1,180 
      do 128 i=1,360
      if (swoa(i,j,k) .le. -10.) swoa(i,j,k)=slevs(i,j,k)
 128  if (twoa(i,j,k) .le. -10.) twoa(i,j,k)=tlevs(i,j,k)
      endif

      write(41) ((twoa(i,j,1),i=1,360),j=1,180),((swoa(i,j,1),i=1,360),j=1,180)


! nzdat are levitus levels
      zlev = dble(nzdat(1))
      do 142 j=1,180
      do 142 i=1,360
      if ((twoa(i,j,1) .lt. -10.0)) twoa(i,j,1) = 0.0
      if ((swoa(i,j,1) .lt. -10.0)) swoa(i,j,1) = 0.0
      twoa(i,j,1)=theta(swoa(i,j,1),twoa(i,j,1),zlev,dble(0.))
 142  rho_data(i,j,1)=r(twoa(i,j,1),swoa(i,j,1),zlev)

! eliminate unstable stratification data

      ssum=0.
      n=0
      lb=1
      lt=2

      do 150 k=2,33

!      read(40)
!     1 ((tlev(i,j,lt),i=1,360),j=1,180),((slev(i,j,lt),i=1,360),j=1,180)
! for density referenced to in situ pressure
      zlev = dble(nzdat(k))
! for density referenced to mean pressure level
!     zlev = float(nzdat(k-1)+nzdat(k))/2.0
!     write(*,*) zlev
      do 145 j=1,180
      do 145 i=1,360
      if ((twoa(i,j,k) .lt. -10.0)) twoa(i,j,k) = 0.0
      if ((swoa(i,j,k) .lt. -10.0)) swoa(i,j,k) = 0.0
      twoa(i,j,k)=theta(swoa(i,j,k),twoa(i,j,k),zlev,dble(0.))
 145  rho_data(i,j,k)=r(twoa(i,j,k),swoa(i,j,k),zlev)
      do 148 j=1,180
      do 148 i=1,360
      rhok=rho_data(i,j,k)
      rhokm=rho_data(i,j,k-1)
      if (rhok.ge.rhokm) go to 148
      n=n+1
      twoa(i,j,k)=twoa(i,j,k-1)
      swoa(i,j,k)=swoa(i,j,k-1)
      rho_data(i,j,k)=rhokm
 148  continue
 150  write(41) ((twoa(i,j,k),i=1,360),j=1,180),((swoa(i,j,k),i=1,360),j=1,180)
      insum=32*180*360

! =================================================
! interpolate levitus data to model horizontal grid
! =================================================
      close(40)
      endif 
      rewind 41     

      do 230 k=1,33
      if (myid .eq. 0) then
      read(41) ((tlev(i,j,1),i=1,360),j=1,180),((slev(i,j,1),i=1,360),j=1,180)
      endif
      call mpi_bcast(tlev,360*180,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(slev,360*180,mpi_real8,0,mpi_comm_world,ierr)


!      if(k.eq.1.or.k.eq.2) then
!      write(*,'(i2,a6,10f10.7)'),myid,' slev=',slev(340,161:170,1)
!      write(*,'(i2,a6,10f10.7)'),myid,' slev=',slev(320,101:110,1)
!      endif


      slon=.5
      slat=-89.5
      ilon=1.
      ilat=1.


      do 230 jj=1,j0
! xlon is model data point, measured in degrees east of grenwich
! ylat is model data point, measured in degrees north of south pole
! ydeg is measured from -90 at south pole to 90 at north pole
! tlev(1,1) is at 0.5 deg e, -89.5 deg lat
! tlev(360,180,) is at 359.5 deg e, 89.5 deg lat
! allowed model latitudes: -89.5.lt.ydeg.lt.89.5
! allowed model longitudes: 0.le.xlon.le.720

!      ylat=90.+ydeg(jj)
!      j=ylat+0.5
      ylat=(ydeg(jj)-slat)*ilat+1.
      j=ylat
      jp=j+1
      do 230 ii=1,i0
!      xlon=xdeg(ii)
      xlon=(xdeg(ii)-slon)*ilon+1.
      if (xlon.gt.360) xlon=xlon-360. 

! xplus is model data point distance east of 0.5 deg west(!)
! note that xplus is always positive due to xlon restrictions

!      xplus=xlon+0.5
      xplus=xlon
      i=xplus      
      ip=i+1
! tmp1 is weighting of east point
      tmp1=xplus-dble(i)
      tmp2=1.-tmp1
 
! special treatments when 359.5 < xlon < 360 and 0 < xlon < 0.5
      if (ip.eq.361) ip=1
!      if (i.eq.0) i=359
      if (i.eq.0) i=360
 
 
      dnw=tlev(i,jp,1)
      dne=tlev(ip,jp,1)
      dsw=tlev(i,j,1)
      dse=tlev(ip,j,1)

      dn=tmp1*dne+tmp2*dnw
      ds=tmp1*dse+tmp2*dsw
      t_data(ii,jj,k)=(ylat-j)*dn+(jp-ylat)*ds


      dnw=slev(i,jp,1)
      dne=slev(ip,jp,1)
      dsw=slev(i,j,1)
      dse=slev(ip,j,1)
      dn=tmp1*dne+tmp2*dnw
      ds=tmp1*dse+tmp2*dsw
 230  s_data(ii,jj,k)=(ylat-j)*dn+(jp-ylat)*ds
      tmax=-100.
      tmin=100.
      do 240 j=2,j1
      do 240 i=2,i1
      tmax=max(tmax,t_data(i,j,1))
 240  tmin=min(tmin,t_data(i,j,1))
!      write(*,245) m,tmin,tmax
! 245  format('month',i3,' tmin,tmax=',2f6.2)

! ===============================================
! interpolate from levitus levels to model levels
! ===============================================
      z2=100.*nzdat(1)
      k=1
      zc=z(2)
      n=1
 320  continue
      if(n.lt.33) then
      n=n+1
      z1=z2
      z2=100.*nzdat(n)
      endif

      if (z2.lt.zc.and.n.lt.33) go to 320
 322  cc1=(z2-zc)/(z2-z1)
      cc2=(zc-z1)/(z2-z1)
!      write(14,323) k,n
! 323  format('for diecast model level',i3,', lower input level=',i3)
      do 325 j=1,j0
      do 325 i=1,i0
      t1_data(i,j,k)=cc1*t_data(i,j,n-1)+cc2*t_data(i,j,n)
 325  s1_data(i,j,k)=cc1*s_data(i,j,n-1)+cc2*s_data(i,j,n)

      k=k+1
      if (k.eq.k0) go to 390
      zc=z(2*k)
      if (zc.gt.z2) go to 320
      go to 322

! write all seasons or months so we can diagnose annual cycle
! 390  write(70) s1_data,t1_data
      
 390  clidata_t1(:,:,:,m) = t1_data
      clidata_s1(:,:,:,m) = s1_data


! wenien
!      open(777,file='/home/opt/s1t1'//grd,form='unformatted')
!      write(777) s1,t1,s,t,swoa,twoa
!      close(777)
!      call mpi_finalize(ierr)
!      stop 390
      enddo
! increment westdeg by one grid interval, and start at second latitude
! (no ghost zones in winds)
      do m=1,12
      if (myid .eq. 0) then
      open(15,file=txc(m))
      open(16,file=tyc(m))
      read(15,400) i,j,x_tmp,x_tmp,x_tmp,x_tmp,fmt
      read(16,400) i,j,x_tmp,x_tmp,x_tmp,x_tmp,fmt
 400  format(2i4/4(3x,f10.7),a10)

! ================
! read hellerman winds
! ================
      read(15,415) tx_tmp
      read(16,415) ty_tmp
      tx_data = tx_tmp
      ty_data = ty_tmp
 415  format(5e15.7)
      write(*,416)
 416  format('read wind stress')

      endif
 
      call mpi_bcast(tx_data,nxhel*nyhel,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(ty_data,nxhel*nyhel,mpi_real8,0,mpi_comm_world,ierr)

! xlon is model data point, measured in degrees east of grenwich
! ylat is model data point, measured in degrees north of south pole
! ydeg is measured from -90 at south pole to 90 at north pole
! tlev(1,1) is at 0.5 deg e, -89.5 deg lat
! tlev(360,180,) is at 359.5 deg e, 89.5 deg lat
! allowed model latitudes: -89.5.lt.ydeg.lt.89.5
! allowed model longitudes: 0.le.xlon.le.720

      slon=.5
      slat=-89.5
      ilon=1.
      ilat=1.

      do jj=2,j1
!      ylat=90.+ydeg(jj)
!      j=ylat+0.5
      ylat=(ydeg(jj)-slat)*ilat+1.
      j=ylat
      jp=j+1
      do ii=2,i1
!      xlon=xdeg(ii+1)
      xlon=(xdeg(ii)-slon)*ilon+1.

      if (xlon.gt.dble(nxhel)) xlon=xlon-dble(nxhel)

!      xplus=xlon+0.5
      xplus=xlon
      i=xplus
      ip=i+1
      tmp1=xplus-dble(i)
      tmp2=1.-tmp1
      if (ip.eq.(nxhel+1)) ip=1
!      if (i.eq.0) i=nxhel-1
      if (i.eq.0) i=nxhel
 
      dnw=tx_data(i,jp)
      dne=tx_data(ip,jp)
      dsw=tx_data(i,j)
      dse=tx_data(ip,j)
      dn=tmp1*dne+tmp2*dnw
      ds=tmp1*dse+tmp2*dsw
!      taux(ii-1,jj-1)=(.5+ylat-j)*dn+(.5-ylat+j)*ds
      taux_data(ii-1,jj-1)=(ylat-j)*dn+(jp-ylat)*ds


      dnw=ty_data(i,jp)
      dne=ty_data(ip,jp)
      dsw=ty_data(i,j)
      dse=ty_data(ip,j)
      dn=tmp1*dne+tmp2*dnw
      ds=tmp1*dse+tmp2*dsw
      tauy_data(ii-1,jj-1)=(ylat-j)*dn+(jp-ylat)*ds

      enddo
      enddo


!     
 
      if(lscowf)then
      if (myid .eq. 0) then
! hellerman
      ierr=nf_get_vara_real(nctx,m+2,istart,icount,sctx_tmp)
      ierr=nf_get_vara_real(ncty,m+2,istart,icount,scty_tmp)
      scowtx_data = sctx_tmp
      scowty_data = scty_tmp 
!dog
!      print*,nf_strerror(ierr)
!      write(1233)tx,ty

      do 420 j=1,nyscow
      do 420 i=1,nxscow
      if(abs(scowtx_data(i,j)).lt.10..and.abs(scowty_data(i,j)).lt.10.)go to 420
      scowtx_data(i,j)=0.
      scowty_data(i,j)=0.
 420  continue
      endif

!      stop

      call mpi_bcast(scowtx_data,nxscow*nyscow,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(scowty_data,nxscow*nyscow,mpi_real8,0,mpi_comm_world,ierr)

! ===================================================
! interpolate scow data to model horizontal grid
! ===================================================
      slon=.125
      slat=-69.875
      ilon=4.
      ilat=4.

      do jj=2,j1
! ============================================================
! tx_data(i,1) is at latitude -69.875; tx_data(i,560) is at latitude 69.875 
! tx_data(1,j) is at longitude 0.125; tx_data(1440,j) is at longitude 359.875
! xlon is measured in degrees east from 0. (xlon=0)
! ydeg is measured in degrees north from -69.875 (ydeg=0)
! ============================================================
! xlon is measured in degrees east from 0 at grenwich
! ydeg is measured from -90 at south pole to 90 at north pole
          
      ylat=(min(ydeg(jj),69.875))*ilat+1.
!      ylat=(ydeg(jj)-slat)*ilat+1.
!      ylat=(70.125+ydeg(jj))*4.
      j=ylat
      if(j.le.0) j=1
      jp=j+1
      do ii=2,i1
      xlon=(xdeg(ii)-slon)*ilon+1.
!      xlon=xdeg(ii)*4.+0.5
!      xlon=xdeg(ii)*4.
      if (xlon.gt.dble(nxscow)) xlon=xlon-dble(nxscow)
! xplus is model data point distance east of 0.5 deg west(!)
! note that xplus is always positive due to xlon restrictions
      xplus=xlon
      i=xplus
      ip=i+1
! tmp1 is weighting of east point
      tmp1=xplus-i
      tmp2=1.-tmp1
 
! special treatments when 359.5 < xlon < 360 and 0 < xlon < 0.5
      if (ip.eq.(nxscow+1)) ip=1
!      if (i.eq.0) i=nxscow-1
      if (i.eq.0) i=nxscow
 
      dnw=scowtx_data(i,jp)
      dne=scowtx_data(ip,jp)
      dsw=scowtx_data(i,j)
      dse=scowtx_data(ip,j)
      dn=tmp1*dne+tmp2*dnw
      ds=tmp1*dse+tmp2*dsw
      tauxs(ii-1,jj-1)=(ylat-j)*dn+(jp-ylat)*ds

      dnw=scowty_data(i,jp)
      dne=scowty_data(ip,jp)
      dsw=scowty_data(i,j)
      dse=scowty_data(ip,j)
      dn=tmp1*dne+tmp2*dnw
      ds=tmp1*dse+tmp2*dsw
      tauys(ii-1,jj-1)=(ylat-j)*dn+(jp-ylat)*ds

      enddo
      enddo
      endif
      
    
      if(lscowf)then

      do i=1,i2
      do j=1,j2
      maxj(i,j)=0
      if(abs(tauxs(i,j)*tauys(i,j)).le.1.e-4)then
      iss=i-1
      iff=i+1
      if(iss.eq.0) iss=1
      if(iff.eq.i1) iff=i2
      jss=j-1
      jff=j+1
      if(jss.eq.0) jss=1
      if(jff.eq.j1) jff=j2
      do ii=iss,iff
      do jj=jss,jff
      maxj(ii,jj)=1
      enddo
      enddo
      tauxs(i,j)=taux_data(i,j)
      tauys(i,j)=tauy_data(i,j)
      endif      
      enddo
      enddo

! lf smoothing !dog ???
      call mpi_type_vector(2*j0,1,i0,mpi_real8,m_vlon3d,ierr)
      call mpi_type_vector(2,i0,i0*j0,mpi_real8,m_vlat3d,ierr)
      call mpi_type_commit(m_vlon3d,ierr)
      call mpi_type_commit(m_vlat3d,ierr)

      js=2
      jf=j1
      if(mylat.eq.0) js=3
      if(mylat.eq.ngy1) jf=j2

      mxpas=10
      do n=1,mxpas

      do i=2,i1
      do j=2,j1
      scr_data(i,j,1)=tauxs(i-1,j-1)
      scr_data(i,j,2)=tauys(i-1,j-1)
      enddo
      enddo
      call mpi_sendrecv(scr_data(2,1,1),1,m_vlon3d,m_w,1,scr_data(i0,1,1),1,m_vlon3d,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(scr_data(i1,1,1),1,m_vlon3d,m_e,1,scr_data(1,1,1),1,m_vlon3d,m_w,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(scr_data(1,2,1),1,m_vlat3d,m_s,1,scr_data(1,j0,1),1,m_vlat3d,m_n,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(scr_data(1,j1,1),1,m_vlat3d,m_n,1,scr_data(1,1,1),1,m_vlat3d,m_s,1,mpi_comm_world,istat,ierr)

      do i=2,i1
      do j=js,jf
      if(maxj(i-1,j-1).eq.1) then
      tauxs(i-1,j-1)=.2*(scr_data(i,j,1)+scr_data(i,j-1,1)+scr_data(i,j+1,1)+scr_data(i-1,j,1)+scr_data(i+1,j,1))
      tauys(i-1,j-1)=.2*(scr_data(i,j,2)+scr_data(i,j-1,2)+scr_data(i,j+1,2)+scr_data(i-1,j,2)+scr_data(i+1,j,2))
      endif
      enddo
      enddo
 
      enddo


!      write(72) tauxs,tauys
      if(myid.eq.0)print*,'scow windstress month:',m


      else
!      write(72) taux_data,tauy_data
      if(myid.eq.0)print*,'hr windstress month:',m
      endif

      cli_tauys(:,:,m) = tauys
      cli_tauxs(:,:,m) = tauxs
      clidata_tauy(:,:,m) = tauy_data
      clidata_taux(:,:,m) = taux_data

      enddo

      call wr8nc_winds(cli_tauxs,cli_tauys,clidata_taux,clidata_tauy) !dog
      call wr8nc_levitus(clidata_s1,clidata_t1) !dog

end subroutine

       
function r(t,s,p)
      p02=1747.4508988+t*(11.51588-0.046331033*t)-s*(3.85429655+0.01353985*t)
      p01=p/10d0+5884.81703666+t*(39.803732+t*(-0.3191477+t*0.0004291133))+2.6126277*s
      r=p01/(p02+0.7028423*p01)
      return
end function


function theta(sal,ts,ps,pr)
    !     ...sal=salinity
    !     ...ts=in situ temp (deg cel)
    !     ...ps=in situ pressure (dbars)
    !     ...pr=reference pressure (dbars)
    delp=pr-ps
    hafp=ps+.5*delp
    delt1=delp*gamma0(sal,ts,ps)
    t1_tmp=ts+.5*delt1
    delt2=delp*gamma0(sal,t1_tmp,hafp)
    t2_tmp=t1_tmp+.2928932*(delt2-delt1)
    delt3=delp*gamma0(sal,t2_tmp,hafp)
    t3_tmp=t2_tmp+1.707107*(delt3-0.5857864*delt2-0.1213203*delt1)
    delt4=delp*gamma0(sal,t3_tmp,pr)
    t4_tmp=t3_tmp+0.16666667*(delt4-6.828427*delt3+4.828427*delt2+delt1)
    theta=t4_tmp
    return
end function


function gamma0(ss,tt,p)
    !     ...adiabati! temperature gradient(deg c/dbar)
    !     ...according to bryden (1973),dsr,20,401-408
    !     ...copied from bio auxilary library
    xx=ss-35
    gamma0=0.35803e-4+0.18932e-5*xx+p*(0.18741e-7- &
         0.11351e-9*xx-0.46206e-12*p) &
         +tt*(0.85258e-5-0.42393e-7*xx+p*(-0.67795e-9+ &
         0.27759e-11*xx+0.18676e-13*p) &
         +tt*(-0.68360e-7+p*(0.87330e-11-0.21687e-15*p) &
         +tt*(0.66228e-9-0.54481e-13*p)))
    return
end function

subroutine wr8nc_depth(rmets)
use const, only: i0,j0,ngx0,ngy0
use comm, only: m_cart,mylat,mylon
implicit none
!declare pnetcdf var
character :: fn*13
integer :: ierr,ncid,nd(2)
integer :: i0id,j0id
integer :: irmt
integer(kind=mpi_offset_kind) :: dim_nn(2)
integer(kind=mpi_offset_kind) :: ns_ij(2),nc_ij(2)
real*8,dimension(i0,j0) :: rmets

      !output file name
      fn = 'data/depth.nc'
      !create ncfile and id
      ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid)

      !define dimension and assign dimID
      dim_nn(1) = i0*ngx0
      dim_nn(2) = j0*ngy0
      ierr=nfmpi_def_dim(ncid,"i0_global",dim_nn(1),i0id)
      ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn(2),j0id)

      !define variables and assign varID
      nd(1) = i0id
      nd(2) = j0id
      ierr=nfmpi_def_var(ncid,"rmets",nf_double,2,nd,irmt)
      ierr=nfmpi_enddef(ncid)

      !put variables
      ns_ij(1) = mylon*i0+1
      ns_ij(2) = mylat*j0+1
      nc_ij(1) = i0
      nc_ij(2) = j0
      ierr=nfmpi_put_vara_double_all(ncid,irmt,ns_ij,nc_ij,rmets)

      ierr=nfmpi_close(ncid)

end subroutine wr8nc_depth

subroutine wr8nc_levitus(clidata_s1,clidata_t1)
use const, only: i0,j0,k1,ngx0,ngy0
use comm, only: m_cart,mylat,mylon,myid
implicit none
!declare pnetcdf var
character :: fn*20
integer :: ierr,ncid,nq(4)
integer :: i0id,j0id,k1id,moid
integer :: is1d,it1d,itmp
integer(kind=mpi_offset_kind) :: dim_nn(4)
integer(kind=mpi_offset_kind) :: ns_ijkt(4),nc_ijkt(4)
real*8,dimension(i0,j0,k1,12),intent(in) :: clidata_s1,clidata_t1
      !output file name
      fn = 'data/annualevitus.nc'
      !create ncfile and id
      ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid)

      !define dimension and assign dimID
      dim_nn(1) = i0*ngx0
      dim_nn(2) = j0*ngy0
      dim_nn(3) = k1
      dim_nn(4) = 12
      ierr=nfmpi_def_dim(ncid,"i0_global",dim_nn(1),i0id)
      ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn(2),j0id)
      ierr=nfmpi_def_dim(ncid,"k1",dim_nn(3),k1id)
      ierr=nfmpi_def_dim(ncid,"month",dim_nn(4),moid)
      !define variables and assign varID
      nq(1) = i0id
      nq(2) = j0id
      nq(3) = k1id
      nq(4) = moid
      ierr=nfmpi_def_var(ncid,"s1_data",nf_double,4,nq,is1d)
      ierr=nfmpi_def_var(ncid,"t1_data",nf_double,4,nq,it1d)
      ierr=nfmpi_enddef(ncid)

      !put variables
      ns_ijkt(1) = mylon*i0+1
      ns_ijkt(2) = mylat*j0+1
      ns_ijkt(3) = 1
      ns_ijkt(4) = 1
      nc_ijkt(1) = i0
      nc_ijkt(2) = j0
      nc_ijkt(3) = k1
      nc_ijkt(4) = 12
      ierr=nfmpi_put_vara_double_all(ncid,is1d,ns_ijkt,nc_ijkt,clidata_s1)
      ierr=nfmpi_put_vara_double_all(ncid,it1d,ns_ijkt,nc_ijkt,clidata_t1)
      ierr=nfmpi_close(ncid)

end subroutine wr8nc_levitus

subroutine wr8nc_winds(cli_tauxs,cli_tauys,clidata_taux,clidata_tauy)
use const, only: i2,j2,ngy0,ngx0
use comm, only: m_cart
implicit none
!declare pnetcdf var
character :: fn*13
integer :: ierr,ncid,nt(3)
integer :: i2id,j2id,moid
integer :: ixs,iys,ixd,iyd
integer(kind=mpi_offset_kind) :: dim_nn(3)
integer(kind=mpi_offset_kind) :: ns_ijk(3),nc_ijk(3)
real*8,dimension(i2,j2,12),intent(in) :: cli_tauxs,cli_tauys
real*8,dimension(i2,j2,12),intent(in) :: clidata_taux,clidata_tauy

      !output file name
      fn = 'data/winds.nc'
      !create ncfile and id
      ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid)

      !define dimension and assign dimID
      dim_nn(1) = i2*ngx0
      dim_nn(2) = j2*ngy0
      dim_nn(3) = 12
      ierr=nfmpi_def_dim(ncid,"i2_global",dim_nn(1),i2id)
      ierr=nfmpi_def_dim(ncid,"j2_global",dim_nn(2),j2id)
      ierr=nfmpi_def_dim(ncid,"month",dim_nn(3),moid)

      !define variables and assign varID
      nt(1) = i2id
      nt(2) = j2id
      nt(3) = moid
      ierr=nfmpi_def_var(ncid,"tauxs",nf_double,3,nt,ixs)
      ierr=nfmpi_def_var(ncid,"tauys",nf_double,3,nt,iys)
      ierr=nfmpi_def_var(ncid,"taux_data",nf_double,3,nt,ixd)
      ierr=nfmpi_def_var(ncid,"tauy_data",nf_double,3,nt,iyd)
      ierr=nfmpi_enddef(ncid)

      !put variables
      ns_ijk(1) = mylon*i2+1
      ns_ijk(2) = mylat*j2+1
      ns_ijk(3) = 1
      nc_ijk(1) = i2
      nc_ijk(2) = j2
      nc_ijk(3) = 12
      ierr=nfmpi_put_vara_double_all(ncid,ixs,ns_ijk,nc_ijk,cli_tauxs)
      ierr=nfmpi_put_vara_double_all(ncid,iys,ns_ijk,nc_ijk,cli_tauys)
      ierr=nfmpi_put_vara_double_all(ncid,ixd,ns_ijk,nc_ijk,clidata_taux)
      ierr=nfmpi_put_vara_double_all(ncid,iyd,ns_ijk,nc_ijk,clidata_tauy)

      ierr=nfmpi_close(ncid)
end subroutine wr8nc_winds

subroutine rdnc_westdeg
use comm, only: m_cart,myid
implicit none
character :: fn*15
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_i,nc_i
real :: test
      fn = 'data/westdeg.nc'
      ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)
      ierr=nfmpi_inq_varid(ncid,'westdeg',itmp)
      ns_i = myid+1
      nc_i = 1
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,westdeg)
      ierr=nfmpi_close(ncid)
end subroutine rdnc_westdeg

subroutine rdnc_yzgrid(j,k,dxmnut)
use comm, only: m_cart,myid,mylat,mylon
use const, only: j0,j1
implicit none
character :: fn*14
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_i,nc_i
integer,intent(inout) :: j,k
real*8,intent(inout) :: dxmnut
      fn = 'data/yzgrid.nc'
      ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)

      ierr=nfmpi_inq_varid(ncid,'j0',itmp)
      ierr=nfmpi_get_var_int_all(ncid,itmp,j)
      ierr=nfmpi_inq_varid(ncid,'k0',itmp)
      ierr=nfmpi_get_var_int_all(ncid,itmp,k)
      ierr=nfmpi_inq_varid(ncid,'dxmnut',itmp)
      ierr=nfmpi_get_var_double_all(ncid,itmp,dxmnut)
      ierr=nfmpi_inq_varid(ncid,'tlz_depth',itmp)
      ierr=nfmpi_get_var_double_all(ncid,itmp,tlz)
      ierr=nfmpi_inq_varid(ncid,'z',itmp)
      ierr=nfmpi_get_var_double_all(ncid,itmp,z)
      ierr=nfmpi_inq_varid(ncid,'odz',itmp)
      ierr=nfmpi_get_var_double_all(ncid,itmp,odz)
      ierr=nfmpi_inq_varid(ncid,'odzw',itmp)
      ierr=nfmpi_get_var_double_all(ncid,itmp,odzw)

      ns_i=mylon*i0+1
      nc_i=i0
      ierr=nfmpi_inq_varid(ncid,'xdeg',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,xdeg)

      ns_i=mylat*j0+1
      nc_i=j0
      ierr=nfmpi_inq_varid(ncid,'y',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,y)
      ierr=nfmpi_inq_varid(ncid,'yv',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,yv)
      ierr=nfmpi_inq_varid(ncid,'yvdeg',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,yvdeg)
      ierr=nfmpi_inq_varid(ncid,'ydeg',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,ydeg)
      ierr=nfmpi_inq_varid(ncid,'cs',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,cs)
      ierr=nfmpi_inq_varid(ncid,'ocs',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,ocs)
      ierr=nfmpi_inq_varid(ncid,'dx',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,dx)
      ierr=nfmpi_inq_varid(ncid,'odx',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,odx)
      ierr=nfmpi_inq_varid(ncid,'dy',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,dy)
      ierr=nfmpi_inq_varid(ncid,'ody',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,ody)

      ns_i=mylat*j1+1
      nc_i=j1
      ierr=nfmpi_inq_varid(ncid,'csv',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,csv)
      ierr=nfmpi_inq_varid(ncid,'ocsv',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,ocsv)
      ierr=nfmpi_inq_varid(ncid,'dxv',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,dxv)
      ierr=nfmpi_inq_varid(ncid,'odxv',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,odxv)
      ierr=nfmpi_inq_varid(ncid,'dyv',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,dyv)
      ierr=nfmpi_inq_varid(ncid,'odyv',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,odyv)
      ierr=nfmpi_close(ncid)
end subroutine rdnc_yzgrid



end module
