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

    integer(kind=mpi_offset_kind) nn,ns_kji,nc_kji,ns_tji,nc_tji,ns_ji,nc_ji,ns_1,nc_1,nf_kji
    common/pnetcdf/nn(6),ns_kji(3),nc_kji(3),ns_tji(3),nc_tji(3),ns_ji(2),nc_ji(2),ns_1,nc_1,nf_kji(3)
    integer(kind=mpi_offset_kind) ns_kji2,nc_kji2,ns_tji2,nc_tji2,ns_ji2,nc_ji2,nf_kji2
    common/pnetcdf2/ns_kji2(3),nc_kji2(3),ns_tji2(3),nc_tji2(3),ns_ji2(2),nc_ji2(2),nf_kji2(3)


contains


subroutine init_io
    open(14,file='data/glo'//grd,form='formatted')
end subroutine


! ----------------------------------------------------------------------
subroutine write_hist
    parameter(i01=i0+1,j01=j0+1,nbg0=nb0*ngy0,nbg1=nbg0-1,ibir=i2t/ngy0)

    character fn*22
    integer*2 nbf

    real*8 rhoavg

    common/histavg/pavg(i0,j0,k1),uavg(i0,j0,k1), &
    vavg(i0,j0,k1),tavg(i0,j0,k1),cavg(i0,j0,k1), &
    savg(i0,j0,k1),wcavg(i0,j0,k0), &
    evavg(i0,j0,k0),hvavg(i0,j0,k0), &
    dmxavg(i0,j0,k1),dmyavg(i0,j0,k1),rhoavg(i0,j0,k1)

    dimension nbf(i2*j2*k1),nda(3),ndc(3),rbf(i2*j2*k1),rbfw(i2*j2*k0),rmbf(i2*j2*12)
    real*8 drbf(i2*j2*k1)
    call timeavg(p,pavg,k1,wopt)
    call timeavg(u2,uavg,k1,wopt)
    call timeavg(v2,vavg,k1,wopt)
    call timeavg(t2,tavg,k1,wopt)
    call timeavg(s2,savg,k1,wopt)
    call timeavg(dmx,dmxavg,k1,wopt)
    call timeavg(dmy,dmyavg,k1,wopt)
    call timeavgd(rho,rhoavg,k1,wopt)
!yc
    call timeavg(c2,cavg,k1,wopt)
!yc
    do 100 k=1,k0
    do 100 j=2,j1
    do 100 i=2,i1
    100  wcavg(i,j,k)=wcavg(i,j,k)+wopt*w(i,j,k)

    do 101 k=2,k1
    do 101 j=2,j1
    do 101 i=2,i1
    evavg(i,j,k)=evavg(i,j,k)+wopt*ev(i-1,j-1,k-1)
    101  hvavg(i,j,k)=hvavg(i,j,k)+wopt*hv(i-1,j-1,k-1)

    if (mod(itf,nopt*itfday) .ne. 0) return
    if(myid.eq.0) print*,'write output'

    ! wenien
    !      open(1111,file='data/his_check'//grd,form='unformatted')
    !      write(1111) uavg,vavg,tavg,savg,t2,s2
    !      close(1111)
    !      call mpi_finalize(ierr)
    !      stop 666

    write(fn,'(a11,i3.3,a1,i3.3,a3)'),'../output/data',nyr,'_',n360,'.nc'

    ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid)
    ierr=nfmpi_def_dim(ncid, "level-center", nn(1), kcid)
    ierr=nfmpi_def_dim(ncid, "level-face", nn(5), kfid)
    ierr=nfmpi_def_dim(ncid, "latitude", nn(2), jid)
    ierr=nfmpi_def_dim(ncid, "longitude", nn(3), iid)
    ierr=nfmpi_def_dim(ncid, "month", nn(4), mid)

    nda(3) = kcid
    nda(2) = jid
    nda(1) = iid
    ierr=nfmpi_def_var(ncid, "pressure", nf_double, 3, nda, idp)
    ierr=nfmpi_def_var(ncid, "u-velocity", nf_double, 3, nda, idu)
    ierr=nfmpi_def_var(ncid, "v-velocity", nf_double, 3, nda, idv)
    ierr=nfmpi_def_var(ncid, "temperature", nf_double, 3, nda, idt)
    ierr=nfmpi_def_var(ncid, "salinity", nf_double, 3, nda, ids)
    ierr=nfmpi_def_var(ncid, "dmx", nf_double, 3, nda, idmx)
    ierr=nfmpi_def_var(ncid, "dmy", nf_double, 3, nda, idmy)
    ierr=nfmpi_def_var(ncid, "rho", nf_double, 3, nda, irho)
!yc
    ierr=nfmpi_def_var(ncid, "cfc11", nf_double, 3, nda, idc)
!yc
    ndc(3) = kfid
    ndc(2) = jid
    ndc(1) = iid
    ierr=nfmpi_def_var(ncid, "w-velocity", nf_double, 3, ndc, idw)
    ierr=nfmpi_def_var(ncid, "hv", nf_double, 3, ndc, ihv)
    ierr=nfmpi_def_var(ncid, "ev", nf_double, 3, ndc, iev)
    ndc(3) = mid
    ndc(2) = jid
    ndc(1) = iid
    ierr=nfmpi_def_var(ncid, "qavg", nf_double, 3, ndc, iqavg)

    ierr=nfmpi_enddef(ncid)

    call mtx2vec_c(pavg,in,rbf)
    ierr=nfmpi_put_vara_double_all(ncid, idp, ns_kji2, nc_kji2, rbf)

    call mtx2vec_c(uavg,in,rbf)
    ierr=nfmpi_put_vara_double_all(ncid, idu, ns_kji2, nc_kji2, rbf)

    call mtx2vec_c(vavg,in,rbf)
    ierr=nfmpi_put_vara_double_all(ncid, idv, ns_kji2, nc_kji2, rbf)

    call mtx2vec_c(tavg,in,rbf)
    ierr=nfmpi_put_vara_double_all(ncid, idt, ns_kji2, nc_kji2, rbf)

    call mtx2vec_c(savg,in,rbf)
    ierr=nfmpi_put_vara_double_all(ncid, ids, ns_kji2, nc_kji2, rbf)
!yc
    call mtx2vec_c(cavg,in,rbf)
    ierr=nfmpi_put_vara_double_all(ncid, idc, ns_kji2, nc_kji2, rbf)
!yc
    call mtx2vec_c(dmxavg,in,rbf)
    ierr=nfmpi_put_vara_double_all(ncid, idmx, ns_kji2, nc_kji2, rbf)

    call mtx2vec_c(dmyavg,in,rbf)
    ierr=nfmpi_put_vara_double_all(ncid, idmy, ns_kji2, nc_kji2, rbf)

    call mtx2dec_c(rhoavg,in,drbf)
    ierr=nfmpi_put_vara_double_all(ncid,irho,ns_kji2,nc_kji2,drbf)

    call mtx2vec_f(wcavg,iw,rbfw)
    ierr=nfmpi_put_vara_double_all(ncid, idw, ns_kji2, nf_kji2, rbfw)

    call mtx2vec_f(hvavg,iw,rbfw)
    ierr=nfmpi_put_vara_double_all(ncid, ihv, ns_kji2, nf_kji2, rbfw)

    call mtx2vec_f(evavg,iw,rbfw)
    ierr=nfmpi_put_vara_double_all(ncid, iev, ns_kji2, nf_kji2, rbfw)

    call mtx2rvec2(qavg,rmbf,i2,j2,i2,j2,12,0,0)
    ierr=nfmpi_put_vara_double_all(ncid, iqavg, ns_tji2, nc_tji2, rmbf)

    call pnetcdferrck(ierr)

    ierr=nfmpi_close(ncid)

    pavg  = 0.
    uavg  = 0.
    vavg  = 0.
    tavg  = 0.
    savg  = 0.
!yc
    cavg  = 0.
!yc
    wcavg = 0.
    hvavg = 0.
    evavg = 0.
    dmxavg= 0.
    dmyavg= 0.
    rhoavg= 0.
end subroutine

! ----------------------------------------------------------------------
subroutine write_topo
    parameter(i01=i0+1,j01=j0+1,nbg0=nb0*ngy0,nbg1=nbg0-1,ibir=i2t/ngy0)

    character fn*17
    integer*2 nbf

    common/windsglo/taux(i2,j2,12),tauy(i2,j2,12)

    dimension z_c(k1),z_f(k0),nda(3),ndb(3),ndc(3),ndd(4),ndw(3)

    integer(kind=mpi_offset_kind) st_z,ct_z,nss(4),nee(4),nws(3),nwe(3),nks(2),nke(2)

    dimension xlon(i2),ylat(j2)

    dimension wbf(i2*j2),sbf(i2*j2*k1),rbf(i2*j2*12),vbf(i2*j2*k1*12)

    dimension inf(i2*j2*k1),iwf(i2*j2*k0),kbf(i2*j2)


    fn='../output/topo.nc'

    do 999 j=1,j2
    999 ylat(j)=ydeg(j+1)
    do 1000 i=1,i2
    1000 xlon(i)=xdeg(i+1)
    do 1001 k=1,k1
    1001 z_c(k)=z(2*k)
    do 1002 k=1,k0
    1002 z_f(k)=z(1+2*(k-1))


    ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid)
    ierr=nfmpi_def_dim(ncid, "level-center", nn(1), kcid)
    ierr=nfmpi_def_dim(ncid, "level-face", nn(5), kfid)
    ierr=nfmpi_def_dim(ncid, "latitude", nn(2), jid)
    ierr=nfmpi_def_dim(ncid, "longitude", nn(3), iid)
    ierr=nfmpi_def_dim(ncid, "month", nn(4), mid)


    ierr=nfmpi_def_var(ncid, "lev_c", nf_double, 1, kcid, izc)
    ierr=nfmpi_def_var(ncid, "lev_f", nf_double, 1, kfid, izf)
    ierr=nfmpi_def_var(ncid, "lat", nf_double, 1, jid, iyd)
    ierr=nfmpi_def_var(ncid, "lon", nf_double, 1, iid, ixd)

    nda(1)=iid
    nda(2)=jid
    ierr=nfmpi_def_var(ncid, "s_nudge", nf_double, 2, nda, isn)
    ierr=nfmpi_def_var(ncid, "t_nudge", nf_double, 2, nda, itn)

    ndb(1)=iid
    ndb(2)=jid
    ndb(3)=kcid
    ierr=nfmpi_def_var(ncid, "hmean_nudge", nf_double, 3, ndb, ihn)

    ndc(1)=iid
    ndc(2)=jid
    ndc(3)=mid
    ierr=nfmpi_def_var(ncid, "taux", nf_double, 3, ndc, itx)
    ierr=nfmpi_def_var(ncid, "tauy", nf_double, 3, ndc, ity)

    ndd(1)=iid
    ndd(2)=jid
    ndd(3)=kcid
    ndd(4)=mid
    ierr=nfmpi_def_var(ncid, "sclis", nf_double, 4, ndd, isc)
    ierr=nfmpi_def_var(ncid, "tclis", nf_double, 4, ndd, itc)

    ierr=nfmpi_def_var(ncid, "kb", nf_int, 2, ndb, iik)
    ierr=nfmpi_def_var(ncid, "mask", nf_int,3, ndb, iin)

    ndw(1)=iid
    ndw(2)=jid
    ndw(3)=kfid
    ierr=nfmpi_def_var(ncid, "mask_w", nf_int,3, ndw, iiw)

    ierr=nfmpi_enddef(ncid)

    ! z! zf ydeg xdeg
    st_z=1
    ct_z = k1
    ierr=nfmpi_put_vara_double_all(ncid,izc,st_z,ct_z,z_c)
    ct_z = k0
    ierr=nfmpi_put_vara_double_all(ncid,izf,st_z,ct_z,z_f)
    ierr=nfmpi_put_vara_double_all(ncid,iyd,ns_kji(2),nc_kji(2),ylat)
    ierr=nfmpi_put_vara_double_all(ncid,ixd,ns_kji(3),nc_kji(3),xlon)
    ! s_nudge t_nudge hmean_nudge
    call mtx2rvec2(s_nudge,wbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_put_vara_double_all(ncid, isn, ns_ji2,nc_ji2,wbf)
    call mtx2rvec2(t_nudge,wbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_put_vara_double_all(ncid, itn, ns_ji2, nc_ji2,wbf)
    call mtx2rvec2(hmean_nudge,sbf,i0,j0,i2,j2,k1,1,1)
    ierr=nfmpi_put_vara_double_all(ncid, ihn, ns_kji2, nc_kji2,sbf)
    ! taux tauy
    call mtx2rvec2(taux,rbf,i2,j2,i2,j2,12,0,0)
    ierr=nfmpi_put_vara_double_all(ncid, itx, ns_tji2, nc_tji2, rbf) 
    call mtx2rvec2(tauy,rbf,i2,j2,i2,j2,12,0,0)
    ierr=nfmpi_put_vara_double_all(ncid, ity, ns_tji2, nc_tji2, rbf)

    ! sclis tclis
    nss(4)=ns_tji(1)
    nss(3)=ns_kji(1)
    nss(2)=ns_tji(2)
    nss(1)=ns_tji(3)
    nee(4)=nc_tji(1)
    nee(3)=nc_kji(1)
    nee(2)=nc_tji(2)
    nee(1)=nc_tji(3)
    call mtx2rvec_m(sclis,vbf,i0,j0,k1,i2,j2,k1,12,1,1)
    ierr=nfmpi_put_vara_double_all(ncid,isc,nss,nee,vbf)    
    call mtx2rvec_m(tclis,vbf,i0,j0,k1,i2,j2,k1,12,1,1)  
    ierr=nfmpi_put_vara_double_all(ncid,itc,nss,nee,vbf)


    ! kb in iw
    nks(1)=ns_kji2(1)
    nks(2)=ns_kji2(2)
    nke(1)=nc_kji2(1)
    nke(2)=nc_kji2(2)
    call mtx2ivec(kb,kbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_put_vara_int_all(ncid,iik,ns_kji2,nc_kji2,kbf)

    call mtx2ivec(in,inf,i0,j0,i2,j2,k1,1,1)
    ierr=nfmpi_put_vara_int_all(ncid,iin,ns_kji2,nc_kji2,inf)

    nws(3)=ns_kji(1)
    nws(2)=ns_kji(2)
    nws(1)=ns_kji(3)
    nwe(3)=nn(5)
    nwe(2)=nc_kji(2)
    nwe(1)=nc_kji(3)
    call mtx2ivec(iw,iwf,i0,j0,i2,j2,k0,1,1)
    ierr=nfmpi_put_vara_int_all(ncid,iiw,nws,nwe,iwf)


    ierr=nfmpi_close(ncid)
end subroutine


! ----------------------------------------------------------------------
subroutine write_rerun
    parameter(nbg0=nb0*ngy0,nbg1=nbg0-1,i01=i0+1,ibir=i2t/ngy0,j01=j0+1)

    character fn*17
    integer*2 nbf
    real*8 dbf

    dimension dbf(i2*j2)
    dimension wbf(i2*j2*k0)

    dimension nd(4)


    !      write(*,'(a20,12i3)') 'kilmer:nsombo=', nsombo
    ierr=nfmpi_create(m_cart,'data/rerun.nc',nf_64bit_offset,mpi_info_null,ncid)
    call pnetcdferrck(ierr)
    ierr=nfmpi_def_dim(ncid, "latitude", nn(2), idj)
    ierr=nfmpi_def_dim(ncid, "longitude", nn(3), idi)
    ierr=nfmpi_def_dim(ncid, "month", nn(4), idt)
    ierr=nfmpi_def_dim(ncid, "global" , nn(6), id1)
    call pnetcdferrck(ierr)
    nd(1)=id1
    ierr=nfmpi_def_var(ncid, "itf", nf_int, 1, nd, i_itf)
    ierr=nfmpi_def_var(ncid, "days", nf_double, 1, nd, i_days)
    ierr=nfmpi_def_var(ncid, "evap", nf_double, 1, nd, i_evap)
    ierr=nfmpi_def_var(ncid, "qsum", nf_double, 1, nd, i_qsum)
    ierr=nfmpi_def_var(ncid, "olde", nf_double, 1, nd, i_olde)
    ierr=nfmpi_def_var(ncid, "oldq", nf_double, 1, nd, i_oldq)
    nd(1)=idt
    ierr=nfmpi_def_var(ncid, "nsombo", nf_int, 1, nd, i_nsombo)
    nd(1)=idj
    nd(2)=idi
    ierr=nfmpi_def_var(ncid, "x", nf_double, 2, nd, ix)
    !       ierr=nfmpi_def_var(ncid, "fsf", nf_double, 2, nd, ifsf)
    ierr=nfmpi_def_var(ncid, "p0", nf_double, 2, nd, ip0)
    nd(1)=idt
    nd(2)=idj
    nd(3)=idi
    ierr=nfmpi_def_var(ncid, "qvag", nf_double, 3, nd, iqavg)
    nd(1)=idj
    nd(2)=idi
    ierr=nfmpi_def_var(ncid, "tnudge", nf_double, 2, nd, itnudge)
    ierr=nfmpi_def_var(ncid, "pbar", nf_double, 2, nd, ipbar)
    ierr=nfmpi_def_var(ncid, "pvar", nf_double, 2, nd, ipvar)
    ierr=nfmpi_def_var(ncid, "xbar", nf_double, 2, nd, ixbar)
    ierr=nfmpi_def_var(ncid, "ucli", nf_double, 2, nd, iucli)
    ierr=nfmpi_def_var(ncid, "vcli", nf_double, 2, nd, ivcli)
    ierr=nfmpi_def_var(ncid, "rmsv", nf_double, 2, nd, irmsv)
    ierr=nfmpi_enddef(ncid)
    ierr=nfmpi_put_vara_int_all(ncid, i_itf, ns_1, nc_1, itf)
    ierr=nfmpi_put_vara_double_all(ncid, i_days, ns_1, nc_1, days)
    ierr=nfmpi_put_vara_double_all(ncid, i_evap, ns_1, nc_1, evap)
    ierr=nfmpi_put_vara_double_all(ncid, i_qsum, ns_1, nc_1, qsum)
    ierr=nfmpi_put_vara_double_all(ncid, i_olde, ns_1, nc_1, olde)
    ierr=nfmpi_put_vara_double_all(ncid, i_oldq, ns_1, nc_1, oldq)
    ierr=nfmpi_put_vara_int_all(ncid, i_nsombo, ns_tji(1), nc_tji(1), nsombo)
    call mtx2dvec(x,dbf)
    ierr=nfmpi_put_vara_double_all(ncid,ix,ns_ji,nc_ji,dbf)
    !       call mtx2rvec(fsf,wbf,i0,j0,i2,j2,1,1,1)
    !       ierr=nfmpi_put_vara_double_all(ncid, ifsf, ns_ji, nc_ji, wbf)
    call mtx2rvec(p0,wbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_put_vara_double_all(ncid, ip0, ns_ji, nc_ji, wbf)
    call mtx2rvec(qavg,wbf,i2,j2,i2,j2,12,0,0)
    ierr=nfmpi_put_vara_double_all(ncid, iqavg, ns_tji, nc_tji, wbf)
    call mtx2rvec(tnudge,wbf,i2,j2,i2,j2,1,0,0)
    ierr=nfmpi_put_vara_double_all(ncid, itnudge, ns_ji, nc_ji, wbf)
    call mtx2rvec(pbar,wbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_put_vara_double_all(ncid, ipbar, ns_ji, nc_ji, wbf)
    call mtx2rvec(pvar,wbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_put_vara_double_all(ncid, ipvar, ns_ji, nc_ji, wbf)
    call mtx2rvec(xbar,wbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_put_vara_double_all(ncid, ixbar, ns_ji, nc_ji, wbf)
    call mtx2rvec(ucli,wbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_put_vara_double_all(ncid, iucli, ns_ji, nc_ji, wbf)
    call mtx2rvec(vcli,wbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_put_vara_double_all(ncid, ivcli, ns_ji, nc_ji, wbf)
    call mtx2rvec(rmsv,wbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_put_vara_double_all(ncid, irmsv, ns_ji, nc_ji, wbf)
    ierr=nfmpi_close(ncid)

    call write2nc('u1','u1',u1,k1,1)
    call write2nc('u2','u2',u2,k1,1)
    call write2nc('ul','ulf',ulf,k1,1)
    call write2nc('v1','v1',v1,k1,1)
    call write2nc('v2','v2',v2,k1,1)
    call write2nc('vl','vlf',vlf,k1,1)
    call write2nc('t1','t1',t1,k1,1)
    call write2nc('t2','t2',t2,k1,1)
    call write2nc('tl','tlf',tlf,k1,1)
    call write2nc('s1','s1',s1,k1,1)
    call write2nc('s2','s2',s2,k1,1)
    call write2nc('sl','slf',slf,k1,1)
    call write2nc('pp','p',p,k1,1)
    call write2nc('uc','u',u,k1,1)
    call write2nc('wc','w',w,k0,5)
    call write2nc('vc','v',v,k1,1)
!yc
    call write2nc('c1','c1',c1,k1,1)
    call write2nc('c2','c2',c2,k1,1)
    call write2nc('cl','clf',clf,k1,1)
!yc
end subroutine

! ----------------------------------------------------------------------
subroutine read_gas(fice,xkw,patm)
!
    dimension fice(i2,j2,12),xkw(i2,j2,12),patm(i2,j2,12)
    parameter(nbg0=nb0*ngy0,nbg1=nbg0-1,i01=i0+1,ibir=i2t/ngy0,j01=j0+1)
    dimension rrbf(i2*j2*12)
    real*8 :: rrbf
!
    ierr=nfmpi_open(m_cart,'data/timcom_gasx.nc',nf_nowrite,mpi_info_null,ncid)
    call pnetcdferrck(ierr)
!
    ierr=nfmpi_inq_varid(ncid, "fice", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_tji,nc_tji,rrbf)

    call pnetcdferrck(ierr)
    call rvec2mtx(fice,rrbf,i2,j2,i2,j2,12,0,0)
!
    ierr=nfmpi_inq_varid(ncid, "xkw", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_tji,nc_tji,rrbf)
    call pnetcdferrck(ierr)
    call rvec2mtx(xkw,rrbf,i2,j2,i2,j2,12,0,0)
!
    ierr=nfmpi_inq_varid(ncid, "patm", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_tji,nc_tji,rrbf)
    call pnetcdferrck(ierr)
    call rvec2mtx(patm,rrbf,i2,j2,i2,j2,12,0,0)
!
end subroutine

! ----------------------------------------------------------------------


! ----------------------------------------------------------------------
subroutine read_rerun
    parameter(nbg0=nb0*ngy0,nbg1=nbg0-1,i01=i0+1,ibir=i2t/ngy0,j01=j0+1)
    character fn*17
    integer*2 nbf

    real*8 dbf

    dimension dbf(i2*j2)
    dimension rbf(i2*j2*k0),rrbf(i2*j2*12)
    dimension nd(4)

    ierr=nfmpi_open(m_cart,'data/rerun.nc',nf_nowrite,mpi_info_null,ncid)
    call pnetcdferrck(ierr)

    ierr=nfmpi_inq_varid(ncid, "itf", itmp)
    ierr=nfmpi_get_var_int_all(ncid, itmp, itf)
    call pnetcdferrck(ierr)
    ierr=nfmpi_inq_varid(ncid, "days", itmp)
    ierr=nfmpi_get_var_double_all(ncid, itmp, days)
    call pnetcdferrck(ierr)
    ierr=nfmpi_inq_varid(ncid, "evap", itmp)
    ierr=nfmpi_get_var_double_all(ncid, itmp, evap)
    call pnetcdferrck(ierr)
    ierr=nfmpi_inq_varid(ncid, "qsum", itmp)
    ierr=nfmpi_get_var_double_all(ncid, itmp, qsum)
    call pnetcdferrck(ierr)
    ierr=nfmpi_inq_varid(ncid, "olde", itmp)
    ierr=nfmpi_get_var_double_all(ncid, itmp, olde)
    call pnetcdferrck(ierr)
    ierr=nfmpi_inq_varid(ncid, "oldq", itmp)
    ierr=nfmpi_get_var_double_all(ncid, itmp, oldq)
    call pnetcdferrck(ierr)
    ierr=nfmpi_inq_varid(ncid, "nsombo", itmp)
    ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_tji(1),nc_tji(1),nsombo)
    call pnetcdferrck(ierr)
    ierr=nfmpi_inq_varid(ncid, "x", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ji,nc_ji,dbf)
    call pnetcdferrck(ierr)
    call dvec2mtx(x,dbf)
    !       ierr=nfmpi_inq_varid(ncid, "fsf", itmp)
    !       ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ji,nc_ji,rbf)
    !       call rvec2mtx(fsf,rbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_inq_varid(ncid, "p0", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ji,nc_ji,rbf)
    call pnetcdferrck(ierr)
    call rvec2mtx(p0,rbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_inq_varid(ncid, "qvag", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_tji,nc_tji,rrbf)
    call pnetcdferrck(ierr)
    call rvec2mtx(qavg,rrbf,i2,j2,i2,j2,12,0,0)
    ierr=nfmpi_inq_varid(ncid, "tnudge", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ji,nc_ji,rbf)
    call pnetcdferrck(ierr)
    call rvec2mtx(tnudge,rbf,i2,j2,i2,j2,1,0,0)
    ierr=nfmpi_inq_varid(ncid, "pbar", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ji,nc_ji,rbf)
    call pnetcdferrck(ierr)
    call rvec2mtx(pbar,rbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_inq_varid(ncid, "pvar", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ji,nc_ji,rbf)
    call pnetcdferrck(ierr)
    call rvec2mtx(pvar,rbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_inq_varid(ncid, "xbar", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ji,nc_ji,rbf)
    call pnetcdferrck(ierr)
    call rvec2mtx(xbar,rbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_inq_varid(ncid, "ucli", itmp)

    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ji,nc_ji,rbf)
    call pnetcdferrck(ierr)
    call rvec2mtx(ucli,rbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_inq_varid(ncid, "vcli", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ji,nc_ji,rbf)
    call pnetcdferrck(ierr)
    call rvec2mtx(vcli,rbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_inq_varid(ncid, "rmsv", itmp)
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ji,nc_ji,rbf)
    call pnetcdferrck(ierr)
    call rvec2mtx(rmsv,rbf,i0,j0,i2,j2,1,1,1)
    ierr=nfmpi_close(ncid)
    call pnetcdferrck(ierr)

    call readnc('u1','u1',u1,k1)
    call readnc('u2','u2',u2,k1)
    call readnc('ul','ulf',ulf,k1)
    call readnc('v1','v1',v1,k1)
    call readnc('v2','v2',v2,k1)
    call readnc('vl','vlf',vlf,k1)
    call readnc('t1','t1',t1,k1)
    call readnc('t2','t2',t2,k1)
    call readnc('tl','tlf',tlf,k1)
    call readnc('s1','s1',s1,k1)
    call readnc('s2','s2',s2,k1)
    call readnc('sl','slf',slf,k1)
    call readnc('pp','p',p,k1)
    call readnc('uc','u',u,k1)
    call readnc('wc','w',w,k0)
    call readnc('vc','v',v,k1)
!yc
    call readnc('c1','c1',c1,k1)
    call readnc('c2','c2',c2,k1)
    call readnc('cl','clf',clf,k1)   
!yc
    call mpi_sendrecv(u(i1,1,1),1,m_vlon,m_e,1,u(1,1,1),1,m_vlon,m_w,1,m_cart,istat,ierr)
    call mpi_sendrecv(v(1,j1,1),1,m_vlat,m_n,1,v(1,1,1),1,m_vlat,m_s,1,m_cart,istat,ierr)
end subroutine


! -------------------------------------------------------------
subroutine write_init2nc(stm,ttm)
    dimension stm(i0,j0,k1),ttm(i0,j0,k1),wbf(i2*j2*k1),nd(3)

    ierr=nfmpi_create(m_cart,'data/initial.nc',nf_64bit_offset,mpi_info_null,ncid)

    call pnetcdferrck(ierr)
    ierr=nfmpi_def_dim(ncid, "level", nn(1), idk)
    ierr=nfmpi_def_dim(ncid, "latitude", nn(2), idj)
    ierr=nfmpi_def_dim(ncid, "longitude", nn(3), idi)
    call pnetcdferrck(ierr)
    nd(1)=idk
    nd(2)=idj
    nd(3)=idi
    ierr=nfmpi_def_var(ncid, 'ti', nf_double, 3, nd, itmp)
    ierr=nfmpi_def_var(ncid, 'si', nf_double, 3, nd, ismp)
    ierr=nfmpi_enddef(ncid)

    ! for temperature
    l=1
    do 100 i=2,i1
    do 100 j=2,j1
    do 100 k=1,k1
    wbf(l)=ttm(i,j,k)
    100  l=l+1
    nc_kji(1)=k1
    ierr=nfmpi_put_vara_double_all(ncid, itmp, ns_kji, nc_kji, wbf)
    call pnetcdferrck(ierr)
    ! for salinity
    l=1
    do 101 i=2,i1
    do 101 j=2,j1
    do 101 k=1,k1
    wbf(l)=stm(i,j,k)
    101  l=l+1
    ierr=nfmpi_put_vara_double_all(ncid, ismp, ns_kji, nc_kji, wbf)
    ierr=nfmpi_close(ncid)
end subroutine


! ----------------------------------------------------------------------
subroutine initbred
    parameter(nbg0=nb0*ngy0,nbg1=nbg0-1,i01=i0+1,ibir=i2t/ngy0,j01=j0+1)

    character fn*17,filename*30

    real*8 dbf

    dimension dbf(i2*j2)
    dimension rbf(i2*j2*k0),rrbf(i2*j2*12)
    dimension nd(4)


    dimension ub(i0,j0,k1),vb(i0,j0,k1)

    jt=mylat*j2
    it=mylon*i2
    do 50 j=1,j2
    do 50 i=1,i2
    ibred(i,j)=0
    50   if (((it+i.ge.540).and.(it+i).le.720).and.((jt+j.ge.524).and.(jt+j).le.572))ibred(i,j)=1




    do 100 k=1,30
    do 100 j=2,j1
    do 100 i=2,i1
    ub(i,j,k)=u1(i,j,k)
    100  vb(i,j,k)=v1(i,j,k)

    call readnc('u1','u1',u1,k1)
    call readnc('v1','v1',v1,k1)
    call readnc('t1','t1',t1,k1)
    call readnc('s1','s1',s1,k1)

    bf=0.
    eke=0.
    tmp=0.
    do 150 k=1,30
    do 150 j=2,j1
    do 150 i=2,i1
    ub(i,j,k)=ibred(i-1,j-1)*(ub(i,j,k)-u1(i,j,k))
    vb(i,j,k)=ibred(i-1,j-1)*(vb(i,j,k)-v1(i,j,k))
    temp=in(i,j,k)*(ub(i,j,k)*ub(i,j,k)+vb(i,j,k)*vb(i,j,k))
    150  tmp=tmp+temp
    call mpi_allreduce(tmp,temp,1,mpi_real8,mpi_sum,m_cart,ierr)
    bf1=sqrt(temp)

    if (myid .eq. 0) print*,'bred ke1 factor =',bf1,'bred ke1 = ',temp

    do 200 k=1,30
    do 200 j=2,j1
    do 200 i=2,i1
    ub(i,j,k)=u2(i,j,k)
    200  vb(i,j,k)=v2(i,j,k)

    call readnc('u2','u2',u2,k1)
    call readnc('v2','v2',v2,k1)
    call readnc('t2','t2',t2,k1)
    call readnc('s2','s2',s2,k1)

    bf=0.
    eke=0.
    tmp=0.
    do 250 k=1,30
    do 250 j=2,j1
    do 250 i=2,i1
    ub(i,j,k)=ibred(i-1,j-1)*(ub(i,j,k)-u2(i,j,k))
    vb(i,j,k)=ibred(i-1,j-1)*(vb(i,j,k)-v2(i,j,k))
    temp=in(i,j,k)*(ub(i,j,k)*ub(i,j,k)+vb(i,j,k)*vb(i,j,k))
    250  tmp=tmp+temp
    call mpi_allreduce(tmp,temp,1,mpi_real8,mpi_sum,m_cart,ierr)
    bf2=sqrt(temp)
    if (myid .eq. 0) print*,'bred ke2 factor =',bf2,'bred ke2 = ',temp

    do 300 k=1,30
    do 300 j=2,j1
    do 300 i=2,i1
    ub(i,j,k)=ulf(i,j,k)
    300  vb(i,j,k)=vlf(i,j,k)

    call readnc('ul','ulf',ulf,k1)
    call readnc('vl','vlf',vlf,k1)
    call readnc('tl','tlf',tlf,k1)
    call readnc('sl','slf',slf,k1)

    bf=0.
    eke=0.
    tmp=0.
    do 350 k=1,30
    do 350 j=2,j1
    do 350 i=2,i1
    ub(i,j,k)=ibred(i-1,j-1)*(ub(i,j,k)-ulf(i,j,k))
    vb(i,j,k)=ibred(i-1,j-1)*(vb(i,j,k)-vlf(i,j,k))
    temp=in(i,j,k)*(ub(i,j,k)*ub(i,j,k)+vb(i,j,k)*vb(i,j,k))
    350  tmp=tmp+temp
    call mpi_allreduce(tmp,temp,1,mpi_real8,mpi_sum,m_cart,ierr)
    bflf=sqrt(temp)

    if (myid .eq. 0)print*,'bred kelf factor =',bflf,'bred kelf = ',temp
end subroutine


! ----------------------------------------------------------------------
subroutine rescbred
    parameter(nbg0=nb0*ngy0,nbg1=nbg0-1,i01=i0+1,ibir=i2t/ngy0,j01=j0+1)

    character fn*17,filename*30

    real*8 dbf
    dimension dbf(i2*j2)
    dimension rbf(i2*j2*k0),rrbf(i2*j2*12)
    dimension nd(4)

    dimension ub(i0,j0,k1),vb(i0,j0,k1),tb(i0,j0,k1),sb(i0,j0,k1)

    do 100 k=1,30
    do 100 j=2,j1
    do 100 i=2,i1
    ub(i,j,k)=u1(i,j,k)
    vb(i,j,k)=v1(i,j,k)
    tb(i,j,k)=t1(i,j,k)
    100  sb(i,j,k)=s1(i,j,k)

    call readnc('u1','u1',u1,k1)
    call readnc('v1','v1',v1,k1)
    call readnc('t1','t1',t1,k1)
    call readnc('s1','s1',s1,k1)

    tmp=0.

    do 130 k=1,30
    do 130 j=2,j1
    do 130 i=2,i1
    ub(i,j,k)=ibred(i-1,j-1)*(ub(i,j,k)-u1(i,j,k))
    vb(i,j,k)=ibred(i-1,j-1)*(vb(i,j,k)-v1(i,j,k))
    tb(i,j,k)=ibred(i-1,j-1)*(tb(i,j,k)-t1(i,j,k))
    sb(i,j,k)=ibred(i-1,j-1)*(sb(i,j,k)-s1(i,j,k))
    temp=in(i,j,k)*(ub(i,j,k)*ub(i,j,k)+vb(i,j,k)*vb(i,j,k))
    130  tmp=tmp+temp


    call mpi_allreduce(tmp,temp,1,mpi_real8,mpi_sum,m_cart,ierr)

    tp=bf1/sqrt(temp)
    temp=0.
    do 160 k=1,30
    do 160 j=2,j1
    do 160 i=2,i1
    temp=temp+in(i,j,k)*ibred(i-1,j-1)*ub(i,j,k)*ub(i,j,k)
    temp=temp+in(i,j,k)*ibred(i-1,j-1)*vb(i,j,k)*vb(i,j,k)
    u1(i,j,k)=u1(i,j,k)+in(i,j,k)*ub(i,j,k)*tp
    v1(i,j,k)=v1(i,j,k)+in(i,j,k)*vb(i,j,k)*tp
    t1(i,j,k)=t1(i,j,k)+in(i,j,k)*tb(i,j,k)*tp
    160  s1(i,j,k)=s1(i,j,k)+in(i,j,k)*sb(i,j,k)*tp

    call mpi_allreduce(temp,tmp,1,mpi_real8,mpi_sum,m_cart,ierr)
    if (myid .eq. 0) print*,'bred ke1 factor =',tp,'bred ke1 = ',tmp


    do 200 k=1,30
    do 200 j=2,j1
    do 200 i=2,i1
    ub(i,j,k)=u2(i,j,k)
    vb(i,j,k)=v2(i,j,k)
    tb(i,j,k)=t2(i,j,k)
    200  sb(i,j,k)=s2(i,j,k)

    call readnc('u2','u2',u2,k1)
    call readnc('v2','v2',v2,k1)
    call readnc('t2','t2',t2,k1)
    call readnc('s2','s2',s2,k1)

    tmp=0.

    do 230 k=1,30
    do 230 j=2,j1
    do 230 i=2,i1
    ub(i,j,k)=ibred(i-1,j-1)*(ub(i,j,k)-u2(i,j,k))
    vb(i,j,k)=ibred(i-1,j-1)*(vb(i,j,k)-v2(i,j,k))
    tb(i,j,k)=ibred(i-1,j-1)*(tb(i,j,k)-t2(i,j,k))
    sb(i,j,k)=ibred(i-1,j-1)*(sb(i,j,k)-s2(i,j,k))
    temp=in(i,j,k)*(ub(i,j,k)*ub(i,j,k)+vb(i,j,k)*vb(i,j,k))
    230  tmp=tmp+temp

    call mpi_allreduce(tmp,temp,1,mpi_real8,mpi_sum,m_cart,ierr)

    tp=bf2/sqrt(temp)
    temp=0.
    do 260 k=1,30
    do 260 j=2,j1
    do 260 i=2,i1
    temp=temp+in(i,j,k)*ibred(i-1,j-1)*ub(i,j,k)*ub(i,j,k)
    temp=temp+in(i,j,k)*ibred(i-1,j-1)*vb(i,j,k)*vb(i,j,k)
    u2(i,j,k)=u2(i,j,k)+in(i,j,k)*ub(i,j,k)*tp
    v2(i,j,k)=v2(i,j,k)+in(i,j,k)*vb(i,j,k)*tp
    t2(i,j,k)=t2(i,j,k)+in(i,j,k)*tb(i,j,k)*tp
    260  s2(i,j,k)=s2(i,j,k)+in(i,j,k)*sb(i,j,k)*tp

    call mpi_allreduce(temp,tmp,1,mpi_real8,mpi_sum,m_cart,ierr)
    if (myid .eq. 0) print*,'bred ke2 factor =',tp,'bred ke2 = ',tmp

    do 300 k=1,30
    do 300 j=2,j1
    do 300 i=2,i1
    ub(i,j,k)=ulf(i,j,k)
    vb(i,j,k)=vlf(i,j,k)
    tb(i,j,k)=tlf(i,j,k)
    300  sb(i,j,k)=slf(i,j,k)

    call readnc('ul','ulf',ulf,k1)
    call readnc('vl','vlf',vlf,k1)
    call readnc('tl','tlf',tlf,k1)
    call readnc('sl','slf',slf,k1)

    tmp=0.

    do 330 k=1,30
    do 330 j=2,j1
    do 330 i=2,i1
    ub(i,j,k)=ibred(i-1,j-1)*(ub(i,j,k)-ulf(i,j,k))
    vb(i,j,k)=ibred(i-1,j-1)*(vb(i,j,k)-vlf(i,j,k))
    tb(i,j,k)=ibred(i-1,j-1)*(tb(i,j,k)-tlf(i,j,k))
    sb(i,j,k)=ibred(i-1,j-1)*(sb(i,j,k)-slf(i,j,k))
    temp=in(i,j,k)*(ub(i,j,k)*ub(i,j,k)+vb(i,j,k)*vb(i,j,k))
    330  tmp=tmp+temp

    call mpi_allreduce(tmp,temp,1,mpi_real8,mpi_sum,m_cart,ierr)

    tp=bflf/sqrt(temp)
    temp=0.
    do 360 k=1,30
    do 360 j=2,j1
    do 360 i=2,i1
    temp=temp+in(i,j,k)*ibred(i-1,j-1)*ub(i,j,k)*ub(i,j,k)
    temp=temp+in(i,j,k)*ibred(i-1,j-1)*vb(i,j,k)*vb(i,j,k)
    ulf(i,j,k)=ulf(i,j,k)+in(i,j,k)*ub(i,j,k)*tp
    vlf(i,j,k)=vlf(i,j,k)+in(i,j,k)*vb(i,j,k)*tp
    tlf(i,j,k)=tlf(i,j,k)+in(i,j,k)*tb(i,j,k)*tp
    360  slf(i,j,k)=slf(i,j,k)+in(i,j,k)*sb(i,j,k)*tp

    call mpi_allreduce(temp,tmp,1,mpi_real8,mpi_sum,m_cart,ierr)
    if (myid .eq. 0) print*,'bred kel factor =',tp,'bred kel = ',tmp
end subroutine

! ----------------------------------------------------------------------
subroutine write2nc(fn,vn,fld,l0,nl)
    parameter(nbg0=nb0*ngy0,nbg1=nbg0-1,i01=i0+1,ibir=i2t/ngy0,j01=j0+1)

    character fn*(*),vn*(*),filename*10

    dimension fld(i0,j0,l0),wbf(i2*j2*l0),nd(3)

    write(filename,'(a5,a2,a3)'),'data/',fn,'.nc'

    ierr=nfmpi_create(m_cart,filename,nf_64bit_offset,mpi_info_null,ncid)
    if (ierr .ne. 0) print*,fn,nfmpi_strerror(ierr)

    if (nl .eq. 1) then
	ierr=nfmpi_def_dim(ncid, "level", nn(1), idk)
    else
	ierr=nfmpi_def_dim(ncid, "level", nn(nl), idk)
    endif
    ierr=nfmpi_def_dim(ncid, "latitude", nn(2), idj)
    ierr=nfmpi_def_dim(ncid, "longitude", nn(3), idi)
    nd(1)=idk
    nd(2)=idj
    nd(3)=idi
    ierr=nfmpi_def_var(ncid, vn, nf_double, 3, nd, itmp)
    if (ierr .ne. 0) print*,fn,nfmpi_strerror(ierr)

    ierr=nfmpi_enddef(ncid)
    if (ierr .ne. 0) print*,fn,nfmpi_strerror(ierr)

    l=1
    do 100 i=2,i1
    do 100 j=2,j1
    do 100 k=1,l0
    wbf(l)=fld(i,j,k)
    100  l=l+1

    nc_kji(1)=l0
    ierr=nfmpi_put_vara_double_all(ncid, itmp, ns_kji, nc_kji, wbf)

    if (ierr .ne. 0) print*,fn,nfmpi_strerror(ierr)

    ierr=nfmpi_close(ncid)
    if (ierr .ne. 0) print*,fn,nfmpi_strerror(ierr)
end subroutine


! -------------------------------------------------------------
subroutine read_initnc(stm,ttm)
    dimension stm(i0,j0,k1),ttm(i0,j0,k1),rbf(i2*j2*k1),nd(3)

    ierr=nfmpi_open(m_cart,'data/initial.nc',nf_nowrite,mpi_info_null,ncid)
    if(ierr.ne.0)then
	call pnetcdferrck(ierr)
	call mpi_finalize(ierr)
	stop
    endif
    ierr=nfmpi_inq_varid(ncid,'ti',itmp)
    nc_kji(1)=k1
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_kji,nc_kji,rbf)
    call pnetcdferrck(ierr)
    l=1
    do 11 i=2,i1
    do 11 j=2,j1
    do 11 k=1,k1
    ttm(i,j,k)=rbf(l)
    11   l=l+1

    ierr=nfmpi_inq_varid(ncid,'si',ismp)
    ierr=nfmpi_get_vara_double_all(ncid,ismp,ns_kji,nc_kji,rbf)
    call pnetcdferrck(ierr)
    l=1
    do 12 i=2,i1
    do 12 j=2,j1
    do 12 k=1,k1
    stm(i,j,k)=rbf(l)
    12   l=l+1
    ierr=nfmpi_close(ncid)
end subroutine


! ----------------------------------------------------------------------
subroutine readnc(fn,vn,fld,l0)
    parameter(nbg0=nb0*ngy0,nbg1=nbg0-1,i01=i0+1,ibir=i2t/ngy0,j01=j0+1)

    character fn*(*),vn*(*),filename*10

    dimension fld(i0,j0,l0),rbf(i2*j2*l0),nd(3)

    write(filename,'(a4,a2,a3)'),'data/',fn,'.nc'
    ierr=nfmpi_open(m_cart,filename,nf_nowrite,mpi_info_null,ncid)
    if (ierr .ne. 0) print*,fn,nfmpi_strerror(ierr)
    ierr=nfmpi_inq_varid(ncid,vn,itmp)
    if (ierr .ne. 0) print*,'nfmpi_inq_var ',fn,nfmpi_strerror(ierr)
    nc_kji(1)=l0
    ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_kji,nc_kji,rbf)
    if (ierr .ne. 0) print*,'nfmpi_get_vara',fn,nfmpi_strerror(ierr)

    ierr=nfmpi_close(ncid)
    if (ierr .ne. 0) print*,fn,nfmpi_strerror(ierr)

    l=1
    do 11 i=2,i1
    do 11 j=2,j1
    do 11 k=1,l0
    fld(i,j,k)=rbf(l)
    11   l=l+1
end subroutine


! ----------------------------------------------------------------------
subroutine timeavg(fd,fld,kl,ct)
    dimension fld(i0,j0,k1),fd(i0,j0,kl)

    do 100 k=1,kl
    do 100 j=2,j1
    do 100 i=2,i1
    100  fld(i,j,k)=fld(i,j,k)+ct*fd(i,j,k)
end subroutine


! ----------------------------------------------------------------------
subroutine timeavgd(fd,fld,kl,ct)
    real*8 fld,fd
    dimension fld(i0,j0,k1),fd(i0,j0,kl)
    do 100 k=1,kl
    do 100 j=2,j1
    do 100 i=2,i1
    100  fld(i,j,k)=fld(i,j,k)+ct*fd(i,j,k)
end subroutine


! ----------------------------------------------------------------------
subroutine mtx2rvec(fd,rbf,in0,jn0,nx,ny,nl,ip,jp)
    dimension fd(in0,jn0,nl),rbf(nx*ny*nl)
    l=1
    ibegin=1+ip
    iend=nx+ibegin-1
    jb=1+jp
    je=ny+jb-1
    do 11 i=ibegin,iend
    do 11 j=jb,je
    do 11 k=1,nl
    rbf(l)=fd(i,j,k)
    11   l=l+1
end subroutine


! ----------------------------------------------------------------------
subroutine mtx2rvec2(fd,rbf,in0,jn0,nx,ny,nl,ip,jp)
    dimension fd(in0,jn0,nl),rbf(nx*ny*nl)
    l=1
    ibegin=1+ip
    iend=nx+ibegin-1
    jb=1+jp
    je=ny+jb-1
    do 11 k=1,nl
    do 11 j=jb,je
    do 11 i=ibegin,iend
    rbf(l)=fd(i,j,k)
    11   l=l+1
end subroutine


! ----------------------------------------------------------------------
subroutine mtx2rvec_m(fd,rbf,in0,jn0,kn1,nx,ny,nk,nl,ip,jp)
    dimension fd(in0,jn0,kn1,nl),rbf(nx*ny*nk*nl)
    l=1
    ibegin=1+ip
    iend=nx+ibegin-1
    jb=1+jp
    je=ny+jb-1
    do 12 m=1,nl
    do 12 k=1,nk
    do 12 j=jb,je
    do 12 i=ibegin,iend
    rbf(l)=fd(i,j,k,m)
    12   l=l+1
end subroutine


! ----------------------------------------------------------------------
subroutine mtx2ivec(fd,rbf,in0,jn0,nx,ny,nl,ip,jp)
    integer*2 fd
    integer rbf
    dimension fd(in0,jn0,nl),rbf(nx*ny*nl)
    l=1
    ibegin=1+ip
    iend=nx+ibegin-1
    jb=1+jp
    je=ny+jb-1
    do 11 k=1,nl
    do 11 j=jb,je
    do 11 i=ibegin,iend
    rbf(l)=int(fd(i,j,k))
    11   l=l+1
end subroutine


! ----------------------------------------------------------------------
subroutine mtx2dvec(fd,dbf)
    parameter(i2j2=i2*j2)
    real*8 fd,dbf
    dimension fd(i0,j0),dbf(i2j2)
    l=1
    do 11 i=2,i1
    do 11 j=2,j1
    dbf(l)=fd(i,j)
    11   l=l+1
end subroutine


! ----------------------------------------------------------------------
subroutine rvec2mtx(fd,rbf,in0,jn0,nx,ny,nl,ip,jp)
    dimension fd(in0,jn0,nl),rbf(nx*ny*nl)
    ibegin=1+ip
    iend=nx+ibegin-1
    jb=1+jp
    je=ny+jb-1
    l=1
    do 11 i=ibegin,iend
    do 11 j=jb,je
    do 11 k=1,nl
    fd(i,j,k)=rbf(l)
    11   l=l+1
end subroutine


! ----------------------------------------------------------------------
subroutine rvec2mtx2(fd,rbf,in0,jn0,nx,ny,nl,ip,jp)
    dimension fd(nx,ny,nl),rbf(nx*ny*nl)
    ibegin=1+ip
    iend=nx+ibegin-1
    jb=1+jp
    je=ny+jb-1
    l=1
    do 11 i=ibegin,iend
    do 11 j=jb,je
    do 11 k=1,nl
    fd(i-1,j-1,k)=rbf(l)
    11   l=l+1
end subroutine


! ----------------------------------------------------------------------
subroutine mtx2vec_c(fd,in,rbf)
    parameter(i2j2k1=i2*j2*k1)
    integer*2 in,nbf
    dimension fd(i0,j0,k1)
    dimension in(i0,j0,k1),rbf(i2j2k1)
    !      rm=rmin+1.e-5*rf
    !      rn=9999./rf
    l=1
    do 11 k=1,k1
    do 11 j=2,j1
    do 11 i=2,i1
    !      fd(i,j,k)=max(min(fd(i,j,k),rmin+rf),rmin)
    !      nbf(l)=-999*(1-in(i,j,k))+in(i,j,k)*rn*(fd(i,j,k)-rm)
    rbf(l)=fd(i,j,k)*in(i,j,k)
    ! nbf(l)=in(i,j,k)*rn*(fd(i,j,k)-rmin)
    11   l=l+1
end subroutine


! ----------------------------------------------------------------------
subroutine mtx2vec_f(fd,iw,rbfw)
    parameter(i2j2k0=i2*j2*k0)
    integer*2 iw
    dimension fd(i0,j0,k0)
    dimension iw(i0,j0,k0),rbfw(i2j2k0)
    !      rm=rmin+1.e-5*rf
    !      rn=9999./rf
    l=1
    do 11 k=1,k0
    do 11 j=2,j1
    do 11 i=2,i1
    !      fd(i,j,k)=max(min(fd(i,j,k),rmin+rf),rmin)
    !      nbf(l)=-999*(1-in(i,j,k))+in(i,j,k)*rn*(fd(i,j,k)-rm)
    rbfw(l)=fd(i,j,k)*iw(i,j,k)
    ! nbf(l)=in(i,j,k)*rn*(fd(i,j,k)-rmin)
    11   l=l+1
end subroutine


! ----------------------------------------------------------------------
subroutine mtx2dec_c(fd,in,rbf)
    parameter(i2j2k1=i2*j2*k1)
    integer*2 in,nbf
    real*8 fd
    dimension fd(i0,j0,k1)
    dimension in(i0,j0,k1)
    real*8 rbf(i2j2k1)
    l=1
    do 11 k=1,k1
    do 11 j=2,j1
    do 11 i=2,i1
    rbf(l)=in(i,j,k)*fd(i,j,k)
    11   l=l+1
end subroutine



! ----------------------------------------------------------------------
subroutine dvec2mtx(fd,dbf)
    parameter(i2j2=i2*j2)
    real*8 fd,dbf
    dimension fd(i0,j0),dbf(i2j2)
    l=1
    do 11 i=2,i1
    do 11 j=2,j1
    fd(i,j)=dbf(l)
    11   l=l+1
end subroutine


! ----------------------------------------------------------------------
subroutine pnetcdferrck(ierr)
    if(ierr.ne.nf_noerr)write(*,'(/2a)')'pnetcdf error: ',nfmpi_strerror(ierr)
end subroutine


subroutine restart
    if (mod(itf,nopt*itfday).eq.0.or.itf.eq.mxit)then
	if (myid.eq.0) print*,'write rerun file'
	if (lbred .eq. 1) call rescbred
	if (lbred .eq. 0) call write_rerun
    endif

    !  nopt: years of sving rerun files
    if(mod(itf,ncopt*360*itfday).eq.0) then
	ifl_sys=1
	nys=itf/(360*itfday)
	if(myid.eq.0) then
	    write(sys1,'(a11,i3.3,a13,i3.3)'),'test -e opt',nys,' || mkdir opt',nys
	    write(sys2,'(a19,i3.3)'),'cp ./data/*.n! ./opt',nys
	    call system(sys1)
	    call system(sys2)
	    print*,'saving rerun files in opt',nys
	    ifl_sys=0
	endif
	call mpi_bcast(ifl_sys,1,mpi_integer,0,m_cart,ierr)
    endif

    if (itf .ge. mxit) then
	call mpi_finalize(ierr)
	stop
    endif

    write(14,405) n360,nyr
    405  format('save at day',i4,', year',i3)
end subroutine


subroutine finalize_io

end subroutine

subroutine rdnc_evp
use comm, only: m_cart,myid,mylon,mylat
use const, only: i2,j2,ngy0,ngx0
implicit none
character :: fn*19
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_ij(2),nc_ij(2)
      fn = 'prepglo/data/evp.nc'
      ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)
      ns_ij(1) = mylon*i2+1
      ns_ij(2) = mylat*j2+1
      nc_ij(1) = i2
      nc_ij(2) = j2
      ierr=nfmpi_inq_varid(ncid,'al',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,al)
      ierr=nfmpi_inq_varid(ncid,'ar',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,ar)
      ierr=nfmpi_inq_varid(ncid,'ab',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,ab)
      ierr=nfmpi_inq_varid(ncid,'ac',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,ac)
      ierr=nfmpi_inq_varid(ncid,'cl',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,cl)
      ierr=nfmpi_inq_varid(ncid,'cb',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,cb)
      ierr=nfmpi_inq_varid(ncid,'cc',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,cc)
      ierr=nfmpi_inq_varid(ncid,'cr',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,cr)
      ierr=nfmpi_inq_varid(ncid,'ct',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,ct)

      ns_ij(1) = mylon*i2+1
      ns_ij(2) = mylat*(j2+1)+1
      nc_ij(1) = i2
      nc_ij(2) = j2+1
      ierr=nfmpi_inq_varid(ncid,'at',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,at)

      ierr=nfmpi_close(ncid)
end subroutine rdnc_evp

subroutine rdnc_windmix
use comm, only: m_cart
use const, only: k2
implicit none
character :: fn*23
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_i,nc_i
      fn = 'prepglo/data/windmix.nc'
      ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)
      ns_i = 1
      nc_i = k2
      ierr=nfmpi_inq_varid(ncid,'vbk',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,vbk)
      ierr=nfmpi_inq_varid(ncid,'hbk',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,hbk)
      ierr=nfmpi_close(ncid)
end subroutine rdnc_windmix

subroutine rdnc_winds
use comm, only: m_cart,mylon,mylat
use const, only: i2,j2,ngx0,ngy0
implicit none
character :: fn*21
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_ijt(3),nc_ijt(3)
      fn = 'prepglo/data/winds.nc'
      ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)
      ns_ijt(1) = mylon*i2+1
      ns_ijt(2) = mylat*j2+1
      ns_ijt(3) = 1
      nc_ijt(1) = i2
      nc_ijt(2) = j2
      nc_ijt(3) = 12

!      if (lscowf) then
      ierr=nfmpi_inq_varid(ncid,'tauxs',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ijt,nc_ijt,taux)
      ierr=nfmpi_inq_varid(ncid,'tauys',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ijt,nc_ijt,tauy)
!      else
!      ierr=nfmpi_inq_varid(ncid,'taux_data',itmp)
!      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ijt,nc_ijt,taux)
!      ierr=nfmpi_inq_varid(ncid,'tauy_data',itmp)
!      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ijt,nc_ijt,tauy)
!      endif
      ierr=nfmpi_close(ncid)
end subroutine rdnc_winds

subroutine rdnc_levitus(nm)
use comm, only: m_cart,mylon,mylat
use const, only: i0,j0,ngx0,ngy0,k1
implicit none
character :: fn*28
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_ijkt(4),nc_ijkt(4)
integer,intent(in) :: nm
      fn = 'prepglo/data/annualevitus.nc'
      ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)
      ns_ijkt(1) = mylon*i0+1
      ns_ijkt(2) = mylat*j0+1
      ns_ijkt(3) = 1
      ns_ijkt(4) = nm
      nc_ijkt(1) = i0
      nc_ijkt(2) = j0
      nc_ijkt(3) = k1
      nc_ijkt(4) = 1
      ierr=nfmpi_inq_varid(ncid,'s1_data',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ijkt,nc_ijkt,sclim)
      ierr=nfmpi_inq_varid(ncid,'t1_data',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ijkt,nc_ijkt,tclim)
      ierr=nfmpi_close(ncid)
end subroutine rdnc_levitus

subroutine rdnc_zkb
use comm, only: m_cart,mylon,mylat
use const, only: i0,j0,ngx0,ngy0,k1,k0,i01
implicit none
character :: fn*19
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_i,nc_i,ns_ij(2),nc_ij(2)
integer(kind=mpi_offset_kind) :: ns_ijk(3),nc_ijk(3)
      fn = 'prepglo/data/zkb.nc'
      ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)
      ns_i = 1
      nc_i = 1
      ierr=nfmpi_inq_varid(ncid,'tlz',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,tlz)
      ns_i = 1
      nc_i = k0+k1
      ierr=nfmpi_inq_varid(ncid,'z',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,z)
      ns_i = 1
      nc_i = k1
      ierr=nfmpi_inq_varid(ncid,'odz',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,odz)
      ns_i = 1
      nc_i = k0
      ierr=nfmpi_inq_varid(ncid,'odzw',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,odzw)
      ns_ij(1) = mylon*i0+1
      ns_ij(2) = mylat*j0+1
      nc_ij(1) = i0
      nc_ij(2) = j0
      ierr=nfmpi_inq_varid(ncid,'kb',itmp)
      ierr=nfmpi_get_vara_int2_all(ncid,itmp,ns_ij,nc_ij,kb)
      ns_ij(1) = mylon*i0+1
      ns_ij(2) = mylat*(j0+1)+1
      nc_ij(1) = i0
      nc_ij(2) = j0+1
      ierr=nfmpi_inq_varid(ncid,'iv0',itmp)
      ierr=nfmpi_get_vara_int2_all(ncid,itmp,ns_ij,nc_ij,iv0)
      ns_ij(1) = mylon*(i01+1)+1
      ns_ij(2) = mylat*j0+1
      nc_ij(1) = i01+1
      nc_ij(2) = j0
      ierr=nfmpi_inq_varid(ncid,'iu0',itmp)
      ierr=nfmpi_get_vara_int2_all(ncid,itmp,ns_ij,nc_ij,iu0)
      ns_ijk(1) = mylon*i0+1
      ns_ijk(2) = mylat*j0+1
      ns_ijk(3) = 1
      nc_ijk(1) = i0
      nc_ijk(2) = j0
      nc_ijk(3) = k0
      ierr=nfmpi_inq_varid(ncid,'iw',itmp)
      ierr=nfmpi_get_vara_int2_all(ncid,itmp,ns_ijk,nc_ijk,iw)
      nc_ijk(3) = k1
      ierr=nfmpi_inq_varid(ncid,'in',itmp)
      ierr=nfmpi_get_vara_int2_all(ncid,itmp,ns_ijk,nc_ijk,in)
      ns_ijk(2) = mylat*(j0+1)+1
      nc_ijk(2) = j0+1
      ierr=nfmpi_inq_varid(ncid,'iv',itmp)
      ierr=nfmpi_get_vara_int2_all(ncid,itmp,ns_ijk,nc_ijk,iv)
      ns_ijk(1) = mylon*(i01+1)+1
      ns_ijk(2) = mylat*j0+1
      nc_ijk(1) = i01+1
      nc_ijk(2) = j0
      ierr=nfmpi_inq_varid(ncid,'iu',itmp)
      ierr=nfmpi_get_vara_int2_all(ncid,itmp,ns_ijk,nc_ijk,iu)

      ierr=nfmpi_close(ncid)
end subroutine rdnc_zkb

subroutine rdnc_zkb2
use comm, only: m_cart,mylon,mylat
use const, only: i0,j0
implicit none
character :: fn*19
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_ij(2),nc_ij(2)
      fn = 'prepglo/data/zkb.nc'
      ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)
      ns_ij(1) = mylon*i0+1
      ns_ij(2) = mylat*j0+1
      nc_ij(1) = i0
      nc_ij(2) = j0
      ierr=nfmpi_inq_varid(ncid,'rmets',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,scr(:,:,1))
      ierr=nfmpi_close(ncid)
end subroutine rdnc_zkb2

subroutine rdnc_rundata
use comm, only: m_cart,mylon,mylat
use const, only: i0,j0,j1,j2,ngx0,ngy0
implicit none
character :: fn*23
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_i,nc_i
      fn = 'prepglo/data/rundata.nc'
      ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)
      ns_i = 1
      nc_i = 5
      ierr=nfmpi_inq_varid(ncid,'case',itmp)
      ierr=nfmpi_get_vara_text_all(ncid,itmp,ns_i,nc_i,case)
      nc_i = 63
      ierr=nfmpi_inq_varid(ncid,'describ',itmp)
      ierr=nfmpi_get_vara_text_all(ncid,itmp,ns_i,nc_i,dscrib)
      nc_i = 1
      ierr=nfmpi_inq_varid(ncid,'dt',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,dt)
      ierr=nfmpi_inq_varid(ncid,'tlz',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,tlz)
      ierr=nfmpi_inq_varid(ncid,'b',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,b)
      ierr=nfmpi_inq_varid(ncid,'g',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,g)
      ierr=nfmpi_inq_varid(ncid,'dm0',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,dm0)
      ierr=nfmpi_inq_varid(ncid,'de0',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,de0)
      ierr=nfmpi_inq_varid(ncid,'dmz0',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,dmz0)
      ierr=nfmpi_inq_varid(ncid,'rzmx',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,rzmx)
      ierr=nfmpi_inq_varid(ncid,'drag',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,drag)
      ierr=nfmpi_inq_varid(ncid,'fltw',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,fltw)
      ierr=nfmpi_inq_varid(ncid,'daodt',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,daodt)
      ierr=nfmpi_inq_varid(ncid,'odt',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,odt)
      ierr=nfmpi_inq_varid(ncid,'prn',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,prn)
      ierr=nfmpi_inq_varid(ncid,'orzmx',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,orzmx)
      ierr=nfmpi_inq_varid(ncid,'ofltw',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,ofltw)
      ierr=nfmpi_inq_varid(ncid,'rawf',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,rawf)

      ierr=nfmpi_inq_varid(ncid,'ktrm',itmp)
      ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_i,nc_i,ktrm)
      ierr=nfmpi_inq_varid(ncid,'lrstrt',itmp)
      ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_i,nc_i,lrstrt)
      ierr=nfmpi_inq_varid(ncid,'mxit',itmp)
      ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_i,nc_i,mxit)
      ierr=nfmpi_inq_varid(ncid,'isav',itmp)
      ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_i,nc_i,isav)
      ierr=nfmpi_inq_varid(ncid,'mxsav',itmp)
      ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_i,nc_i,mxsav)
      ierr=nfmpi_inq_varid(ncid,'lopen',itmp)
      ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_i,nc_i,lopen)
      ierr=nfmpi_inq_varid(ncid,'lmovi',itmp)
      ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_i,nc_i,lmovi)
      ierr=nfmpi_inq_varid(ncid,'mvi',itmp)
      ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_i,nc_i,mvi)
      ierr=nfmpi_inq_varid(ncid,'lfsrf',itmp)
      ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_i,nc_i,lfsrf)
      ierr=nfmpi_inq_varid(ncid,'lsolver',itmp)
      ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_i,nc_i,lsolver)
      ierr=nfmpi_inq_varid(ncid,'lprecon',itmp)
      ierr=nfmpi_get_vara_int_all(ncid,itmp,ns_i,nc_i,lprecon)
      ns_i = mylat*j2+1
      nc_i = j2
      ierr=nfmpi_inq_varid(ncid,'f',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,f)
      ierr=nfmpi_inq_varid(ncid,'tanphi',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,tanphi)

      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)

      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,'ydeg',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,ydeg)
      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,'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 = 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)


      ierr=nfmpi_close(ncid)
end subroutine rdnc_rundata

end module
