module presolver  
    use const
    use comm
    use var
implicit real*8 (a-h,o-z)
#   include "pnetcdf.inc"
    
    ! preprocessor for 1/2 deg res global model                   july, 2007
    !!cyclic longitudinal b.c.'s (i.e., global or periodic)
    ! ----------------------------------------------------------------------
    !      block data nishal
    ! ----------------------------------------------------------------------
    real*8 hh,d1
    common/contro/r0
    common/lpi/mvfreq
    common/boyanc/d1(i0,j0,k1)
    common/bathy/hh(i0,j0),rmets(i0,j0)
    common/bcs/tsurf(i2,j2,12),ssurf(i2,j2,12),snsp(i2,10,k2,12),sssp(i2,10,k2,12),tnsp(i2,10,k2,12),tssp(i2,10,k2,12)

! lsolver solver control
! lsolver:0 no solver setting
! lsolver:1 evp
! lsolver:2 bi-cgstab
! ie(nb0)=j1
! we initialize ie in do loop 250 of subroutine initfs
!     data ie/6,11,16,...,j1/

! ----------------------------
! explanation of parameters
! ----------------------------
! i0,j0: horizontal dimension parameters (includes ghost zones around
! entire model domain; climatological boundary values are stored in
! ghost zones)
! k0: vertical dimension parameter (number of layer interfaces, equal
! the number of layers or pressure levels plus 1, does not include ghost
! zone)
! nb0: number of evp ellipti! solver restart vectors, used to control
! solver roundoff error (decreases with increasing nb0)
! case=5-character case name
! dscrib=63-character case descriptor
! r0=earth radius
! daodt=number of time steps per day
! b=thermal expansion coefficient when linear equation of state is used
! g=gravity
! dm0,de0 are horizontal heat and momentum diffusivities (cm-cm/sec)
! rzmx=vertical cell reynolds number limit
! drag=drag coefficient (standard value is .002)
! fltw=time filter coefficient on time filtered leapfrog method
! ktrm=thermocline level (used only for animation graphics)
! lrstrt=flag to initialize from restart file; if lrstrt=0, initial vel
!         is zero (but not necessarily boundary values)
! mxit=last time step number of current run (last time step of previous
!      run is saved on restart file)
! lopen=flag for open lateral boundaries (derive geostrophi! inflows)
! lmovi=flag to save data for animation of the results
! lsolver=flag for this program to initialize ellipti! solver data for the
! lprecon=flag for precodition of bi-cgstab(lsolver=2)
!  given geometry and also to generate initial temperature and salinity
! ie=array for evp ellipti! solver (controls roundoff error); ie(nb0)=j1

contains
! ----------------------------------------------------------------------
subroutine solvergen
! ----------------------------------------------------------------------
! input data: etopo5, levitus (t and s)
! the model now reads input wind data directly
! output data:
! input data translated to model variable latitudinal resolution grid

    integer*2 intmp(i2,j2),ing(i0t,j0t),inbuf(i2t*j2t)
    integer*2 isw(j0*k1),ise(j0*k1),isn(i0*k1),iss(i0*k1),irw(j0*k1),ire(j0*k1),irn(i0*k1),irs(i0*k1)
    real*8 scr_tmp(i0,j0) 

    daodt=576
    r0=6.4e8
    b=2.05e-4
    g=980.
    dm0=5.e6
    de0=25.e5
    dmz0=.01
    rzmx=50
    drag=.001
    fltw=.1
    rawf=.53

    ! this may be unstable
    !data ktrm,lrstrt,mxit,isav,mxsav/5,0,0,5400,5400/
   
    ktrm=5
    lrstrt=0
    mxit=0
    isav=5400
    mxsav=5400
    
    !data lopen,lmovi,mvfreq,lfsrf,lsolver,lprecon/0,1,5,1,2,1/
    lopen=0
    lmovi=1
    mvfreq=5
    lfsrf=1
    lsolver=2
    lprecon=1
    case="GLOBO"
    dscrib="global model (.5 deg longitudinal resolution)                  "
    
    if (myid .eq. 0) then
	open(79,file='namelist.txt',form='formatted')
	read(79,*)lrstrt
	read(79,*)daodt
	read(79,*)irun
	close(79)
	mxit=irun*daodt*360
    endif
    call mpi_bcast(lrstrt,1,mpi_integer,0,mpi_comm_world,ierr)
    call mpi_bcast(daodt,1,mpi_real8,0,mpi_comm_world,ierr)
    call mpi_bcast(mxit,1,mpi_integer,0,mpi_comm_world,ierr)

    ! output data
!    open(10,file='data/yzgrid'//grd,form='unformatted')
!    open(11,file='data/depth'//grd,form='unformatted')
!    open(70,file='data/annualevitus'//grd,form='unformatted')
    ! output data
!    open(7,file='data/vbk_ascii'//grd)
!    open(8,file='data/windmix'//grd,form='unformatted')
!    open(9,file='data/rundata'//grd,form='unformatted')
!    open(17,file='data/kbview'//grd)
!    open(60,file='data/zkb'//grd,form='unformatted')
!    open(62,file='data/boundaries'//grd,form='unformatted')
    ! check solver 
    if(lsolver.eq.1 .or. lsolver.eq.2) then
!	open(99,file='data/evp'//grd,form='unformatted')
!	rewind 99
    elseif(lsolver.eq.0) then
	if(myid.eq.0) print*,'no solver setting'
    else
	if(myid.eq.0) print*,'lsolver=',lsolver,' err'
	stop 0
    endif

!    rewind 10
!    rewind 11
!    rewind 70
!    rewind 7 
!    rewind 8 
!    rewind 9 
!    rewind 17 
!    rewind 60 
!    rewind 62

    call rdnc_yzgrid2(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 i0,k0 to values from metgen.f.')
    stop 0
    endif
!   continue   

!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

    !k01=k0+k1
    !write(21,10) (.01*z(k),k=1,k01,2)
    !10   format('z-levels(m)'/(10f8.2))
!coriolis parameter array
      omega2=3.141592654/21600.
! degrees to radians
      pi_180=3.141592654/180.
      do 12 j=2,j1
! spherical curvature parameter
      tanphi(j-1)=tan(pi_180*ydeg(j))/r0
 12   f(j-1)=omega2*sin(pi_180*ydeg(j))

     ! write(21,14) i0,j0
 !14   format('longitude(x) grid dimension:',i4/'latitude(y) grid dimension:',i4/)

! rmets data
! rmets(i,1) corresponds to southern most internal grid point
!      stop
      call rdnc_depth(rmets)
!      read(11) rmets


      write(*,*) '==============',maxval(rmets),minval(rmets)
      write(*,*) 'rmets filtered'
      write(*,*) '=============='
      call mpi_type_vector(j0,1,i0,mpi_real8,m_vlon_2d,ierr) !dog
      call mpi_type_vector(1,i0,j0,mpi_real8,m_vlat_2d,ierr)
      call mpi_type_commit(m_vlon_2d,ierr)
      call mpi_type_commit(m_vlat_2d,ierr)

      call mpi_sendrecv(rmets(2,1),1,m_vlon_2d,m_w,1,rmets(i0,1),1,m_vlon_2d,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(rmets(i1,1),1,m_vlon_2d,m_e,1,rmets(1,1),1,m_vlon_2d,m_w,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(rmets(1,2),1,m_vlat_2d,m_s,1,rmets(1,j0),1,m_vlat_2d,m_n,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(rmets(1,j1),1,m_vlat_2d,m_n,1,rmets(1,1),1,m_vlat_2d,m_s,1,mpi_comm_world,istat,ierr)


!cts change from mxpas=2 to 4
      mxpas=4
      do n=1,mxpas
      do 15 j=2,j1
      do 15 i=2,i1

! don't filter izu ridge and se honshu coastal region
      if(xdeg(i).gt.134.5.and.xdeg(i).lt.147..and.ydeg(j).gt.-62.5.and.ydeg(j).lt.-51.5) goto 15
! don't filter gs window region      
      if(xdeg(i).gt.284.5.and.xdeg(i).lt.320..and.ydeg(j).gt.-54.5.and.ydeg(j).lt.-29.) goto 15

! don't filter regions near land (they are adequately known)
! filter only points farther than nwet grid intervals from land

      if (rmets(i,j).lt. 300.) goto 15
      scr_tmp(i,j)=.25*(rmets(i,j-1)+rmets(i,j+1)+rmets(i-1,j)+rmets(i+1,j))
 15   continue

      do 16 j=2,j1
      do 16 i=2,i1

! don't filter izu ridge and se honshu coastal region
      if(xdeg(i).gt.134.5.and.xdeg(i).lt.147..and.ydeg(j).gt.-62.5.and.ydeg(j).lt.-51.5) goto 16
! don't filter gs window region
      if(xdeg(i).gt.284.5.and.xdeg(i).lt.320..and.ydeg(j).gt.-54.5.and.ydeg(j).lt.-29.) goto 16

      if (rmets(i,j).lt. 300.) goto 16
      rmets(i,j)=scr_tmp(i,j)
 16   continue
     
      call mpi_sendrecv(rmets(2,1),1,m_vlon_2d,m_w,1,rmets(i0,1),1,m_vlon_2d,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(rmets(i1,1),1,m_vlon_2d,m_e,1,rmets(1,1),1,m_vlon_2d,m_w,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(rmets(1,2),1,m_vlat_2d,m_s,1,rmets(1,j0),1,m_vlat_2d,m_n,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(rmets(1,j1),1,m_vlat_2d,m_n,1,rmets(1,1),1,m_vlat_2d,m_s,1,mpi_comm_world,istat,ierr)
      enddo


!cts20150921
      if(.false.)then
! 555 510  64  53
! make izu ridge into a narrow ridge at sill level
!cw      do 18 j=547,527,-1
!cw      jloc=j-mylat*j2
      do 18 j=2,j1
      do 18 i=3,i2
! rmets is in meters
! 15   write(*,16) j,(int(.1*rmets(i,j)),i=555,565)
! 16   format(i3,1x,11i4)
! find shallowest sill level in longitudinal range i=555-565
      kmn=10000
!cw      do 17 i=556,564
!cw      iloc=i-mylon*i2
      
!      if(iloc.le.2.or.iloc.ge.i1.or.jloc.le.2.and.jloc.ge.j1)
!cw      if(iloc.lt.2.or.iloc.gt.i1.or.jloc.lt.2.or.jloc.gt.j1)
!cw     1  goto 18
       if(xdeg(i).gt.138..and.xdeg(i).lt.142..and.ydeg(j).gt.30..and.ydeg(j).lt.35.) then
!      k=kmn
!      kmn=min(kmn,int(.1*rmets(i,j)))
! 17   if (kmn.ne.k) imn=i
! make narrow ridge
!      tmp=max(rmets(imn-1,j),rmets(imn-2,j))
!      rmets(imn-2,j)=tmp
!      rmets(imn-1,j)=tmp
!      tmp=max(rmets(imn+1,j),rmets(imn+2,j))
!      rmets(imn+2,j)=tmp
!      rmets(imn+1,j)=tmp
!     write(*,*) imn,kmn

      k=kmn
!cw      kmn=min(kmn,int(.1*rmets(iloc,jloc)))
      kmn=min(kmn,int(.1*rmets(i,j)))
!cw 17   if (kmn.ne.k) imn=iloc
 17   if (kmn.ne.k) imn=i
! make narrow ridge
      tmp=max(rmets(imn-1,j),rmets(imn-2,j))
      rmets(imn-2,j)=tmp
      rmets(imn-1,j)=tmp
      tmp=max(rmets(imn+1,j),rmets(imn+2,j))
      rmets(imn+2,j)=tmp
      rmets(imn+1,j)=tmp
!cw      tmp=max(rmets(imn-1,jloc),rmets(imn-2,jloc))
!cw      rmets(imn-2,jloc)=tmp
!cw      rmets(imn-1,jloc)=tmp
!cw      tmp=max(rmets(imn+1,jloc),rmets(imn+2,jloc))
!cw      rmets(imn+2,jloc)=tmp
!cw      rmets(imn+1,jloc)=tmp
!     write(*,*) imn,kmn
      endif
 18   continue


      do i=2,i1
      do j=2,j1
! smoother coastline near cozumel
       if(xdeg(i).gt.273.2.and.xdeg(i).lt.273.5.and.ydeg(j).gt.21.1.and.ydeg(j).lt.21.5) then
       rmets(i,j)=0.
       endif

! eliminate shallow water northwest of easternmost florida point
! rmets(1118:1120,516:530)=0.
       if(xdeg(i).gt.279.1.and.xdeg(i).lt.279.9.and.ydeg(j).gt.28.1.and.ydeg(j).lt.31.6) then
       rmets(i,j)=0.
       endif
! rmets(1121,516:530)=28.
       if(xdeg(i).gt.280.0.and.xdeg(i).lt.280.3.and.ydeg(j).gt.28.1.and.ydeg(j).lt.31.6) then
       rmets(i,j)=28.
       endif
! rmets(1121,521:525)=45.
       if(xdeg(i).gt.280.0.and.xdeg(i).lt.280.3.and.ydeg(j).gt.29.2.and.ydeg(j).lt.30.5) then
       rmets(i,j)=45.
       endif
! rmets(1120,522:523)=28.
       if(xdeg(i).gt.279.5.and.xdeg(i).lt.280.0.and.ydeg(j).gt.29.5.and.ydeg(j).lt.30.0) then
       rmets(i,j)=28.
       endif         
       if(ydeg(j).gt.85.0)then
       rmets(i,j)=0.
       endif
       enddo
       enddo

!      stop
! artificial high latitude shelves (shortcircuits arcti! ocean,
! truncates antarcti! seas and bays, and avoids artificial excessive bt
! mode vortex stretching at latitudinal boundaries)

!cw new shallowing scheme
      jtt=min(j2t/60,k2)
!cw      jtt=min(j2t/30,k2)
      kinv=int(dble(k2)/dble(jtt)+.5)
      kzl=int(70.*dble(k1)/60.)
      if(myid.eq.0) print*,'kzl=',kzl

      do 23 j=2,j1

      if(mylat*j2+j-1.le.jtt) then
      jin=mylat*j2+j-2      
      kin=1+kinv*jin
!cw      kk=min(k0+k1,2*kin+3)
      kk=min(k0+k1,2*kin+kzl+1)
      tmp=.01*z(kk)
      do 21 i=1,i0
 21   rmets(i,j)=min(rmets(i,j),tmp)  
      endif

      if(mylat*j2+j.ge.j0t-jtt) then
      jin=j0t-mylat*j2-j-1
      kin=1+kinv*jin
!cw      kk=min(k0+k1,2*kin+3)
      kk=min(k0+k1,2*kin+kzl+1)

      tmp=.01*z(kk)
      do 22 i=1,i0
 22   rmets(i,j)=min(rmets(i,j),tmp)
      endif
 23   continue   


!cts20150921 if(.false.)then
      endif

!      if (mylat .eq. 0) then
!      n=0
!      do 22 k=1,k1
!      n=n+1
!      kk=min(k0+k1,2*k+3)
!      tmp=.01*z(kk)
!      do 22 i=1,i0
! 22   rmets(i,n+1)=min(rmets(i,n+1),tmp)
!      endif
!      
!      if (mylat .eq. ngy1) then
!      n=0
!      do 23 k=1,k1
!      n=n+1
!      kk=min(k0+k1,2*k+3)
!      tmp=.01*z(kk)
!      do 23 i=1,i0
! 23   rmets(i,j0-n)=min(rmets(i,j0-n),tmp)
!      endif

      

! bad point is (39,401) starting at step 91100
!      write(*,24) (j,(int(rmets(i,j)),i=31,47),j=j1,j1-50,-1)
! 24   format(i3,4x,17i4)

!     write(*,24) (j,(int(rmets(i,j)),i=261,300),j=206,167,-1)
!c24   format(/(i3,1x,40i1))

! make all water least 2 layers deep
! to allow baroclini! ventillation all the way to coast
      tmp=.01*z(5)
!convert 1-layer deep water to 2-layer-deep water
! alternatively, !convert water less than 2 layers deep to land

!cw20150910
!cw      tmp1=.01*(z(3)+z(5))/2.
      tmp1=.01*(z(5)+z(7))/2. 

      do 25 j=2,j1
      do 25 i=2,i1
!convert land to water
! 25   if (rmets(i,j).gt.tmp1) rmets(i,j)=max(rmets(i,j),tmp)
!convert water to land
 25   if (rmets(i,j).lt.tmp1) rmets(i,j)=0.

! periodic
      call mpi_sendrecv(rmets(2,1),1,m_vlon_2d,m_w,1,rmets(i0,1),1,m_vlon_2d,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(rmets(i1,1),1,m_vlon_2d,m_e,1,rmets(1,1),1,m_vlon_2d,m_w,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(rmets(1,2),1,m_vlat_2d,m_s,1,rmets(1,j0),1,m_vlat_2d,m_n,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(rmets(1,j1),1,m_vlat_2d,m_n,1,rmets(1,1),1,m_vlat_2d,m_s,1,mpi_comm_world,istat,ierr)

      do 33 j=1,j0
      do 33 i=1,i0
 33   rmets(i,j)=min(tlz,100.*rmets(i,j))
  
! determine "logical grid" rmets array kb(i,j)
! kb(i,j)=number of layers at pressure point (i,j)
      do 34 k=1,k1
      zz=.5*(z(2*k-1)+z(2*k+1))
      do 34 j=1,j0
      do 34 i=1,i0
! if rmets is deeper than mid-layer (zz), increment kb by 1
      l=rmets(i,j)/zz
 34   kb(i,j)=kb(i,j)+min(1,l)

! widen northern izu gap to allow more flow through to entrain more
! oyashio water. this was added at day 90 of year 6. before that the
! kuroshio extension was good, but a small southward loop downstream from
! the gap is decreased by this change.
! deepen widened gap at start of year 14
!     n=1
! further deepen gap at start of year 16
      n=2
      do i=2,i1
      do j=2,j1

      if(.false.)then
!cts 20150911 remove new york bight singularity
       if(xdeg(i).gt.289.5.and.xdeg(i).lt.293.5.and.ydeg(j).gt.42.2.and.ydeg(j).lt.45.2) then
       kb(i,j)=min(kb(i,j),3)
       endif
       if(xdeg(i).gt.289.5.and.xdeg(i).lt.293.5.and.ydeg(j).gt.41.5.and.ydeg(j).lt.45.2) then
       kb(i,j)=min(kb(i,j),7)
       endif
      endif
!cts seamounts around 44.7n, 321e

      if(xdeg(i).gt.319.and.xdeg(i).lt.322.and.ydeg(j).gt.44.8.and.ydeg(j).lt.45.2) then
       kb(i,j)=min(kb(i,j),32)
      endif

      if(.false.)then
!cts kuroshio seperation-japan coast should be sahllower
      if(xdeg(i).gt.138.8.and.xdeg(i).lt.141.4.and.ydeg(j).gt.34.5.and.ydeg(j).lt.35.) then
       kb(i,j)=min(kb(i,j),18)
      endif
      if(xdeg(i).gt.140.9.and.xdeg(i).lt.141.8.and.ydeg(j).gt.35..and.ydeg(j).lt.36.) then
       kb(i,j)=min(kb(i,j),15)
      endif
      if(xdeg(i).gt.140.9.and.xdeg(i).lt.141.8.and.ydeg(j).gt.36.and.ydeg(j).lt.37.) then
       kb(i,j)=min(kb(i,j),11)
      endif
      if(xdeg(i).gt.141.8.and.xdeg(i).lt.142.8.and.ydeg(j).gt.35.1.and.ydeg(j).lt.39.4) then
       kb(i,j)=kb(i,j)-5
      endif
      endif
!cw 20150910
      if(.false.) then

!cw      do j=539,541  i=561,563
!cw      kb(iloc,jloc)=13+n
       if(xdeg(i).gt.140.0.and.xdeg(i).lt.140.8.and.ydeg(j).gt.33.2.and.ydeg(j).lt.34.0) then
       kb(i,j)=13+n
       endif
!cw      iloc=561-mylon*i2    jloc=539-mylat*j2
!cw      kb(iloc,jloc)=13+n
       if(xdeg(i).gt.140.0.and.xdeg(i).lt.140.3.and.ydeg(j).gt.33.2.and.ydeg(j).lt.33.6) then
       kb(i,j)=13+n
       endif
!cw      iloc=562-mylon*i2    jloc=539-mylat*j2
!cw      kb(iloc,jloc)=13+n
       if(xdeg(i).gt.140.3.and.xdeg(i).lt.140.5.and.ydeg(j).gt.33.2.and.ydeg(j).lt.33.6) then
       kb(i,j)=13+n
       endif
!cw      iloc=561-mylon*i2    jloc=540-mylat*j2
!cw      kb(iloc,jloc)=13+n
       if(xdeg(i).gt.140.0.and.xdeg(i).lt.140.3.and.ydeg(j).gt.33.5.and.ydeg(j).lt.33.7) then
       kb(i,j)=13+n
       endif
!cw      iloc=562-mylon*i2    jloc=540-mylat*j2
!cw      kb(iloc,jloc)=13+n
       if(xdeg(i).gt.140.3.and.xdeg(i).lt.140.5.and.ydeg(j).gt.33.5.and.ydeg(j).lt.33.7) then
       kb(i,j)=13+n
       endif
!cw      iloc=562-mylon*i2    jloc=541-mylat*j2
!cw      kb(iloc,jloc)=13+n
       if(xdeg(i).gt.140.3.and.xdeg(i).lt.140.5.and.ydeg(j).gt.33.7.and.ydeg(j).lt.33.9) then
       kb(i,j)=13+n
       endif

! eliminate anomalous ny bight island
!cw      iloc=1142-mylon*i2   jloc=557-mylat*j2
!cw      kb(iloc,jloc)=2
       if(xdeg(i).gt.285.3.and.xdeg(i).lt.285.5.and.ydeg(j).gt.36.9.and.ydeg(j).lt.37.1) then
       kb(i,j)=2
       endif
!cw      iloc=1142-mylon*i2   jloc=558-mylat*j2
!cw      kb(iloc,jloc)=2
       if(xdeg(i).gt.285.3.and.xdeg(i).lt.285.5.and.ydeg(j).gt.37.1.and.ydeg(j).lt.37.3) then
       kb(i,j)=2
       endif
!cw      iloc=1142-mylon*i2   jloc=559-mylat*j2
!cw      kb(iloc,jloc)=2
       if(xdeg(i).gt.285.3.and.xdeg(i).lt.285.5.and.ydeg(j).gt.37.3.and.ydeg(j).lt.37.5) then
       kb(i,j)=2
       endif
!cw      iloc=1144-mylon*i2   jloc=561-mylat*j2
!cw      kb(iloc,jloc)=2
       if(xdeg(i).gt.285.8.and.xdeg(i).lt.290.1.and.ydeg(j).gt.37.7.and.ydeg(j).lt.37.9) then
       kb(i,j)=2
       endif
!cts no above modification
       endif



!cts new check and topographical modification
       if(xdeg(i).gt.140.5.and.xdeg(i).lt.141.5.and.ydeg(j).gt.41.0.and.ydeg(j).lt.42.) then
       kb(i,j)=10
       endif
!cw20150902
       if(xdeg(i).gt.120.5.and.xdeg(i).lt.121.5.and.ydeg(j).gt.24.5.and.ydeg(j).lt.25.0) then
       kb(i,j)=9
       endif
  
       if(xdeg(i).gt.121.5.and.xdeg(i).lt.122.5.and.ydeg(j).gt.25.5.and.ydeg(j).lt.26.0) then
       kb(i,j)=13
       endif

!cw indonesia throughflow
       if(xdeg(i).gt.119.5.and.xdeg(i).lt.120.5.and.ydeg(j).gt.-6..and.ydeg(j).lt.0.25) then
       kb(i,j)=23
       endif

       if(xdeg(i).gt.119.5.and.xdeg(i).lt.121.5.and.ydeg(j).gt.0.5.and.ydeg(j).lt.0.75) then
       kb(i,j)=20
       endif

       if(xdeg(i).gt.119.5.and.xdeg(i).lt.121.5.and.ydeg(j).gt.0.75.and.ydeg(j).lt.1.) then
       kb(i,j)=24
       endif

       if(xdeg(i).gt.119.5.and.xdeg(i).lt.121.5.and.ydeg(j).gt.1..and.ydeg(j).lt.1.5) then
       kb(i,j)=26
       endif

!cw allow lombok
       if(xdeg(i).gt.118.5.and.xdeg(i).lt.119.5.and.ydeg(j).gt.-12.5.and.ydeg(j).lt.-11.75) then
       kb(i,j)=26
       endif

       if(xdeg(i).gt.118.5.and.xdeg(i).lt.119.5.and.ydeg(j).gt.-10.75.and.ydeg(j).lt.-10.25) then
       kb(i,j)=28
       endif

!cw allow ombai
       if(xdeg(i).gt.122.5.and.xdeg(i).lt.123.5.and.ydeg(j).gt.-8.5.and.ydeg(j).lt.-8.) then
       kb(i,j)=26
       endif

       if(xdeg(i).gt.123.5.and.xdeg(i).lt.124.5.and.ydeg(j).gt.-8.5.and.ydeg(j).lt.-8.) then
       kb(i,j)=25
       endif

       if(xdeg(i).gt.126.5.and.xdeg(i).lt.127.5.and.ydeg(j).gt.-9.and.ydeg(j).lt.-8.5) then
       kb(i,j)=24
       endif

       if(xdeg(i).gt.129.5.and.xdeg(i).lt.131.5.and.ydeg(j).gt.-3.75.and.ydeg(j).lt.0.75) then
       kb(i,j)=24
       endif


       enddo
       enddo

       do j=2,j1
       do i=2,i1
       if(ydeg(j).ge.34..and.ydeg(j).le.36.5.and.xdeg(i).ge.352.5.and.xdeg(i).le.356)then
        kb(i,j)=max(kb(i,j),19)
       endif
       enddo
       enddo

       do j=2,j1
       do i=2,i1
       if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
        write(*,*)'topo1',i,j

        do jj=j+6,j-8,-1
        write(*,'(a,11i4)')'topo1',kb(i-5:i+5,jj)
        enddo
       endif
       enddo
       enddo


!cw check
!       do i=2,i1
!       do j=2,j1
!       it=i+mylon*i2-1
!       jt=j+mylat*j2-1
!       if(it.eq.107.and.jt.eq.144) print*,'1 kb 107 144=',kb(i,j)
!     1 ,xdeg(i),ydeg(j)
!       if(it.eq.116.and.jt.eq.144) print*,'1 kb 116 144=',kb(i,j)
!     1 ,xdeg(i),ydeg(j)
!       if(it.eq.117.and.jt.eq.144) print*,'1 kb 117 144=',kb(i,j)
!       enddo
!       enddo

! derive land/sea mask array
      do 52 j=1,j0
      do 52 i=1,i0
      it=kb(i,j)
      if (it.ne.0) go to 50
      kb(i,j)=1
      go to 52
 50   do 51 k=1,it
 51   in(i,j,k)=1
 52   continue

      do 53 j=2,j1
      do 53 i=2,i1
 53   intmp(i-1,j-1)=in(i,j,1)

       do j=2,j1
       do i=2,i1
       if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
        write(*,*)'topo1',i,j
        do jj=j+6,j-8,-1
        write(*,'(a,11i4)')'topo1',in(i-5:i+5,jj,1)
        enddo
       endif
       enddo
       enddo

      nsend=i2*j2

      call mpi_gather(intmp,nsend,mpi_integer2,inbuf,nsend,mpi_integer2,0,mpi_comm_world,ierr)

      if (myid .eq. 0) then
      do 57 nx=1,ngx0
      ngxt=(nx-1)*ngy0
      do 57 ny=1,ngy0
      ngt=((ny-1)+ngxt)*i2*j2
      do 57 j=1,j2
      nt=ngt+(j-1)*i2
      jt=j+(ny-1)*j2
      do 57 i=1,i2
!      it=mod(i+(nx-1)*i2+i2t/2,i2t)
      it=i+(nx-1)*i2
 57   ing(it+1,jt+1)=inbuf(nt+i)
      
! in and ing is the same with in in global_serial1


!remark that boundary condition is computed in value of ghost cell.
! ------------------------------------------------------
! start with single point in central equatorial pacific
! check pt
      n=0
!      ing(i2t/6,j0t/2)=2


      ing(2*i0t/3,j0t/2)=2


! ------------------------------------------------------
! due to periodi! logical grid boundary at 0 deg longitude
! we need a second point in the arcti! ocean
! rmk: i=2 alway lies in processor 0.

!cw20140604      ing(2,j1t)=2

 153  n=n+1
      nsum=0
      do 155 j=2,j1t
      do 154 i=2,i1t
      nn=ing(i+1,j)
      if (ing(i+1,j).eq.1.and.ing(i,j).eq.2) ing(i+1,j)=2
      nsum=nsum+ing(i+1,j)-nn
      nn=ing(i-1,j)
      if (ing(i-1,j).eq.1.and.ing(i,j).eq.2) ing(i-1,j)=2
      nsum=nsum+ing(i-1,j)-nn
      nn=ing(i,j+1)
      if (ing(i,j+1).eq.1.and.ing(i,j).eq.2) ing(i,j+1)=2
      nsum=nsum+ing(i,j+1)-nn
      nn=ing(i,j-1)
      if (ing(i,j-1).eq.1.and.ing(i,j).eq.2) ing(i,j-1)=2
 154  nsum=nsum+ing(i,j-1)-nn
      do 155 i=i1t,2,-1
      nn=ing(i+1,j)
      if (ing(i+1,j).eq.1.and.ing(i,j).eq.2) ing(i+1,j)=2
      nsum=nsum+ing(i+1,j)-nn
      nn=ing(i-1,j)
      if (ing(i-1,j).eq.1.and.ing(i,j).eq.2) ing(i-1,j)=2
      nsum=nsum+ing(i-1,j)-nn
      nn=ing(i,j+1)
      if (ing(i,j+1).eq.1.and.ing(i,j).eq.2) ing(i,j+1)=2
      nsum=nsum+ing(i,j+1)-nn
      nn=ing(i,j-1)
      if (ing(i,j-1).eq.1.and.ing(i,j).eq.2) ing(i,j-1)=2
 155  nsum=nsum+ing(i,j-1)-nn

      write(*,156) nsum,n
 156  format(i8,'changes during iteration',i4)
      if (nsum.ne.0) go to 153


      do 157 j=1,j0t
      ing(i0t,j)=ing(2,j)
 157  ing(1,j)=ing(i1t,j)
      endif
      call mpi_bcast(ing,i0t*j0t,mpi_integer2,0,mpi_comm_world,ierr)

! widen northern izu gap to allow more flow through to entrain more
! oyashio water. this was added at day 90 of year 6. before that the
! kuroshio extension was good, but a small southward loop downstream from
! the gap is decreased by this change.

!cw      iloc=561-mylon*i2
!cw      jloc=539-mylat*j2
!cw      if(iloc.ge.2.and.iloc.le.i1.and.jloc.ge.2.and.jloc.le.j1)
!cw     1 kb(iloc,jloc)=13
!cw      iloc=562-mylon*i2
!cw      jloc=539-mylat*j2
!cw      if(iloc.ge.2.and.iloc.le.i1.and.jloc.ge.2.and.jloc.le.j1)
!cw     1 kb(iloc,jloc)=13
!cw      iloc=561-mylon*i2
!cw      jloc=540-mylat*j2
!cw      if(iloc.ge.2.and.iloc.le.i1.and.jloc.ge.2.and.jloc.le.j1)
!cw     1 kb(iloc,jloc)=13
!cw      iloc=562-mylon*i2
!cw      jloc=540-mylat*j2
!cw      if(iloc.ge.2.and.iloc.le.i1.and.jloc.ge.2.and.jloc.le.j1)
!cw     1 kb(iloc,jloc)=13
!cw      iloc=562-mylon*i2
!cw      jloc=541-mylat*j2
!cw      if(iloc.ge.2.and.iloc.le.i1.and.jloc.ge.2.and.jloc.le.j1)
!cw     1 kb(iloc,jloc)=13


      it=mylon*i2
      jt=mylat*j2
      do 158 j=1,j0
      do 158 i=1,i0
      in(i,j,1)=ing(it+i,jt+j)
      n=in(i,j,1)
 158  if (in(i,j,1).ne.2) kb(i,j)=0
       do j=2,j1
       do i=2,i1
       if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
        write(*,*)'topo2',i,j
        do jj=j+6,j-8,-1
         write(*,'(a,11i4)')'toin2',in(i-5:i+5,jj,1)
        enddo
       endif
       enddo
       enddo

      do 159 k=1,k1
      do 159 j=1,j0
      do 159 i=1,i0
 159  in(i,j,k)=0

       do j=2,j1
       do i=2,i1
       if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
        write(*,*)'topo2',i,j
        do jj=j+6,j-8,-1
         write(*,'(a,11i4)')'topo2',kb(i-5:i+5,jj)
        enddo
       endif
       enddo
       enddo

!     do 8886 j=j1,2,-1
!8886 write(*,8887) j,(kb(i,j),i=1,122)
!8887 format(i3,1x,122i1)
!     do 8888 j=j1,2,-1
!8888 write(*,8887) j,(kb(i,j),i=121,242)
!     do 8889 j=j1,2,-1
!8889 write(*,8887) j,(kb(i,j),i=241,362)

! redefine in with exterraneous pools eliminated

      do 162 j=2,j1
      do 162 i=1,i0
      it=kb(i,j)
      if (it.ne.0) go to 160
      kb(i,j)=1
      go to 162
 160  do 161 k=1,it
 161  in(i,j,k)=1
 162  continue


       do j=2,j1
       do i=2,i1
       if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
        write(*,*)'topo3',i,j
        do jj=j+6,j-8,-1
        write(*,'(a,11i4)')'topo3',kb(i-5:i+5,jj)
        enddo
       endif
       enddo
       enddo

! periodic
      call mpi_type_vector(1,i0,j0,mpi_integer2,m_vlat_2di,ierr)
      call mpi_type_commit(m_vlat_2di,ierr)
      call mpi_sendrecv(kb(1,2),1,m_vlat_2di,m_s,1,kb(1,j0),1,m_vlat_2di,m_n,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(kb(1,j1),1,m_vlat_2di,m_n,1,kb(1,1),1,m_vlat_2di,m_s,1,mpi_comm_world,istat,ierr)

      call mpi_type_vector(k1,i0,i0*j0,mpi_integer2, m_vlat_2di,ierr)
      call mpi_type_commit(m_vlat_2di,ierr)
      call mpi_sendrecv(in(1,2,1),1,m_vlat_2di,m_s,1,in(1,j0,1),1,m_vlat_2di,m_n,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(in(1,j1,1),1,m_vlat_2di,m_n,1,in(1,1,1),1,m_vlat_2di,m_s,1,mpi_comm_world,istat,ierr)


!cw check
!       do i=2,i1
!       do j=2,j1
!       it=i+(mylon)*i2-1
!       jt=j+(mylat)*j2-1
!       if(it.eq.107.and.jt.eq.144) print*,'2 kb 107 144=',kb(i,j)
!       if(it.eq.116.and.jt.eq.144) print*,'2 kb 116 144=',kb(i,j)
!       if(it.eq.117.and.jt.eq.144) print*,'2 kb 117 144=',kb(i,j)
!       enddo
!       enddo






! mask out points having no lateral ventillation
 69   n=0
      do 70 k=2,k1
      do 70 j=2,j1
      do 70 i=2,i1
      if (in(i-1,j,k)+in(i+1,j,k)+in(i,j-1,k)+in(i,j+1,k).ne.0) go to 70
      if (in(i,j,k).ne.0) n=n+1
      in(i,j,k)=0
 70   continue
      write(*,71) n
 71   format('mask out',i4,' points having no lateral ventillation')
      if (n.ne.0) go to 69
! redefine kb based on modified in
      do 72 k=2,k1
      do 72 j=2,j1
      do 72 i=2,i1
      if (in(i,j,k).eq.1) go to 72
      if (in(i,j,k-1).eq.1) kb(i,j)=k-1
! this may be redundant
      kb(i,j)=max(1,kb(i,j))
 72   continue

       do j=2,j1
       do i=2,i1
       if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
        write(*,*)'topo4',i,j
        do jj=j+6,j-8,-1
        write(*,'(a,11i4)')'topo4',kb(i-5:i+5,jj)
        enddo
       endif
       enddo
       enddo

! periodic
      call mpi_barrier(mpi_comm_world,ierr)
      call mpi_type_vector(1,i0,j0,mpi_integer2,m_vlat_2di,ierr)
      call mpi_type_commit(m_vlat_2di,ierr)
      call mpi_type_vector(j0,1,i0,mpi_integer2,m_vlon_2di,ierr)
      call mpi_type_commit(m_vlon_2di,ierr)
      call mpi_sendrecv(kb(1,2),1,m_vlat_2di,m_s,1,kb(1,j0),1,m_vlat_2di,m_n,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(kb(1,j1),1,m_vlat_2di,m_n,1,kb(1,1),1,m_vlat_2di,m_s,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(kb(2,1),1,m_vlon_2di,m_w,1,kb(i0,1),1,m_vlon_2di,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(kb(i1,1),1,m_vlon_2di,m_e,1,kb(1,1),1,m_vlon_2di,m_w,1,mpi_comm_world,istat,ierr)


      call mpi_barrier(mpi_comm_world,ierr)
      call mpi_type_vector(k1,i0,i0*j0,mpi_integer2,m_vlat_2di,ierr)
      call mpi_type_commit(m_vlat_2di,ierr)
      call mpi_type_vector(k1*j0,1,i0,mpi_integer2,m_vlon_2di,ierr)
      call mpi_type_commit(m_vlon_2di,ierr)
      call mpi_sendrecv(in(1,2,1),1,m_vlat_2di,m_s,1,in(1,j0,1),1,m_vlat_2di,m_n,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(in(1,j1,1),1,m_vlat_2di,m_n,1,in(1,1,1),1,m_vlat_2di,m_s,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(in(2,1,1),1,m_vlon_2di,m_w,1,in(i0,1,1),1,m_vlon_2di,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(in(i1,1,1),1,m_vlon_2di,m_e,1,in(1,1,1),1,m_vlon_2di,m_w,1,mpi_comm_world,istat,ierr)
      call mpi_barrier(mpi_comm_world,ierr)



!cw check
!       do i=2,i1
!       do j=2,j1
!       it=i+(mylon)*i2-1
!       jt=j+(mylat)*j2-1
!       if(it.eq.107.and.jt.eq.144) print*,'3 kb 107 144=',kb(i,j)
!       if(it.eq.116.and.jt.eq.144) print*,'3 kb 116 144=',kb(i,j)
!       if(it.eq.117.and.jt.eq.144) print*,'3 kb 117 144=',kb(i,j)
!       enddo
!       enddo


! other mask arrays

      do 64 k=1,k1
      do 64 j=1,j0
      do 64 i=1,i1
 64   iu(i,j,k)=in(i,j,k)*in(i+1,j,k)

      do 65 k=1,k1
      do 65 j=1,j1
      do 65 i=1,i0
 65   iv(i,j,k)=in(i,j,k)*in(i,j+1,k)

      do 66 k=1,k2
      do 66 j=1,j0
      do 66 i=1,i0
 66   iw(i,j,k+1)=in(i,j,k)*in(i,j,k+1)


      do k=1,k1
       jt=(k-1)*j0
       do j=1,j0
        ise(jt+j)=iu(i1,j,k)
        isw(jt+j)=iu(2,j,k)
       enddo
      enddo
      do k=1,k1
       it=(k-1)*i0
       do i=1,i0
        isn(it+i)=iv(i,j2,k)
        iss(it+i)=iv(i,2,k)
       enddo
      enddo

      nlon=j0*k1
      nlat=i0*k1

      call mpi_barrier(mpi_comm_world,ierr)
      call mpi_sendrecv(isw,nlon,mpi_integer2,m_w,1,ire,nlon,mpi_integer2,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(ise,nlon,mpi_integer2,m_e,1,irw,nlon,mpi_integer2,m_w,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(iss,nlat,mpi_integer2,m_s,1,irn,nlat,mpi_integer2,m_n,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(isn,nlat,mpi_integer2,m_n,1,irs,nlat,mpi_integer2,m_s,1,mpi_comm_world,istat,ierr)

      do k=1,k1
       jt=(k-1)*j0
        do j=1,j0
         ise(jt+j)=iu(i2,j,k)
         isw(jt+j)=iu(3,j,k)
         iu(i0,j,k)=ire(jt+j)
         iu(1,j,k)=irw(jt+j)
       enddo
      enddo
      if (mylat .ne.ngy1) then
      do k=1,k1
       it=(k-1)*i0
       do i=1,i0
        iv(i,j0,k)=irn(it+i)
       enddo
      enddo
      endif

      if (mylat .ne. 0) then
      do k=1,k1
       it=(k-1)*i0
       do i=1,i0
        iv(i,0,k)=irs(it+i)
       enddo
      enddo
      endif

      call mpi_barrier(mpi_comm_world,ierr)
      call mpi_sendrecv(isw,nlon,mpi_integer2,m_w,1,ire,nlon,mpi_integer2,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(ise,nlon,mpi_integer2,m_e,1,irw,nlon,mpi_integer2,m_w,1,mpi_comm_world,istat,ierr)

      do k=1,k1
       jt=(k-1)*j0
       do j=1,j0
        iu(i0+1,j,k)=ire(jt+j)
        iu(0,j,k)=irw(jt+j)
       enddo
      enddo
      do j=2,j1
       do i=0,i0+1
        iu0(i,j)=iu(i,j,1)
        iu(i,j,1)=1
       enddo
      enddo
      do j=0,j0
      do i=2,i1
        iv0(i,j)=iv(i,j,1)
        iv(i,j,1)=1
       enddo
      enddo

      if(mylat.eq.0)then
        do i=1,i0
         iv(i,0,1)=0
         iv(i,1,1)=0
        enddo    
      endif
      if(mylat.eq.ngy1)then
        do i=1,i0
         iv(i,j0,1)=0
         iv(i,j1,1)=0
        enddo
      endif                    

!      write(60) tlz,z,odz,odzw,kb,in,iu,iv,iw,iu0,iv0
!      write(60) rmets
      call wr8nc_zkb  !dog
! ===================================================================
! coupling to 1/8 deg resolution tai and gom grids; one may still run
! entire medina model with coupling to global model, because it uses
! no equatorial shelf.
! this global model, the medina model and norpa! tai grid all fit in
! the memory of a single "fast" processor, but may run less than a
! day per hour of cpu. with one way coupling and separate processors
! for each model, the speed will be almost as fast as the full norpac
! duo grid model.
! ===================================================================
! write rmets data for tai grid (at 150 deg e)
! this covers coupling region with tai grid (j0=488 in tai grid)
! from equator northward
!      open(80,file='tai_kbe')
!      write(80,80) (kb(601,j),kb(601,j),j=399,642)
! 80   format(40i3)
!      open(81,file='gom_kbe')
! write rmets data for gom grid (at 300 deg e)
! this covers coupling region with gom grid (j0=336 in gom grid)
! from 10 deg latitude to 46 deg n
!      write(81,80) (kb(1201,j),kb(1201,j),j=439,606)
!      write(*,99) (j,(in(i,j,1),i=1,140),j=j0,1,-1)
! 99   format(i3,1x,140i1)

      call initfs

end subroutine
! ****************
! initialization *
! ****************
! ----------------------------------------------------------------------
subroutine pre(ax,ay,bb,cx,cy,rinv,rinv1,hh,ie)
! ----------------------------------------------------------------------
      real*8 rinv,rinv1,hh,run,rtp,rtemp,rp
      dimension ax(i2,j2),ay(i2,j2),bb(i2,j2),cx(i2,j2),cy(i2,0:j2),rinv(ibir,i0,nbg0),rinv1(ibir,i0,nbg1),hh(i0,j0),ie(nb0),rp(i0t*i0t),rtemp(i2t,i0t),rtp(i2t,i2t),ipvt(i2t)
      dimension axt(i2,j2),ayt(i2,j2),bbt(i2,j2),cxt(i2,j2),cyt(i2,j2),ct(i2)

      do 100 ng=1,ngy0
      jl=1
      if (mylat .eq. ng-1) then
      do 99 j=1,j2
      do 99 i=1,i2
      axt(i,j)=ax(i,j)
      ayt(i,j)=ay(i,j)
      bbt(i,j)=bb(i,j)
      cxt(i,j)=cx(i,j)
 99   cyt(i,j)=cy(i,j)
      endif
      if (mylat .eq. ng-2) then
      do 98 i=1,i2
 98   ct(i)=cy(i,j2)
      endif
      nsd=i2*j2
      call mpi_bcast(axt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
      call mpi_bcast(ayt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
      call mpi_bcast(bbt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
      call mpi_bcast(cxt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
      call mpi_bcast(cyt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
      if (ng .ne. 1)call mpi_bcast(ct,i2,mpi_real8,ng-2,m_comm_lat,ierr)

      if (mylat .eq. ng-1) then
      do 97 i=1,i2
 97   cy(i,0)=ct(i)
      endif
      do 100 nb=1,nb0
      nbg=(ng-1)*nb0+nb
      write(*,111) nbg,nbg0
 111  format ('processing block #',i3,' out of ',i3,' total evp solver')
      jh=ie(nb)
      jhp=jh+1
      jhm=jh-2
      jg=jl+1
      do 250 ii=1,ibir
      ig=mylat*ibir+ii
      idxy=(ig-1)/i2
      do 210 j=jl,jhp
      do 210 i=1,i0
 210  hh(i,j)=0.d0
      if (idxy .eq. mylon) then
      if (mod(ig,i2) .eq. 0) then
      hh(i1,jg)=1.d0
      else
      hh(mod(ig,i2)+1,jg)=1.d0
      endif
      endif
      call mpi_sendrecv(hh(2,jg),1,mpi_real8,m_w,1,hh(i0,jg),1,mpi_real8,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(hh(i1,jg),1,mpi_real8,m_e,1,hh(1,jg),1,mpi_real8,m_w,1,mpi_comm_world,istat,ierr)
      if (nbg .eq. 1) go to 220
      if (idxy .eq. mylon) then
      if (nb .ne. 1) then
       if (mod(ig,i2) .eq. 0) then
       ctemp=cyt(i2,jg-2)
       else
       ctemp=cyt(mod(ig,i2),jg-2)
       endif
      else
       if (mod(ig,i2) .eq. 0) then
       ctemp=ct(i2)
       else
       ctemp=ct(mod(ig,i2))
       endif
      endif
      endif
      call mpi_bcast(ctemp,1,mpi_real8,idxy,m_comm_lon,ierr)
      do 218 n=1,i0
 218  hh(n,jl)=rinv1(ii,n,nbg-1)*ctemp
 220  continue
      do 226 j=jl,jhm
      do 225 i=1,i2
 225  hh(i+1,j+2)=-(axt(i,j)*hh(i,j+1)+ayt(i,j)*hh(i+1,j)+bbt(i,j)*hh(i+1,j+1)+cxt(i,j)*hh(i+2,j+1))/cyt(i,j)
      call mpi_sendrecv(hh(2,j+2),1,mpi_real8,m_w,1,hh(i0,j+2),1,mpi_real8,m_e,1,mpi_comm_world,istat,ierr)
 226  call mpi_sendrecv(hh(i1,j+2),1,mpi_real8,m_e,1,hh(1,j+2),1,mpi_real8,m_w,1,mpi_comm_world,istat,ierr)
      j=jh-1
      do 230 i=1,i2
 230  rinv(ii,i+1,nbg)=axt(i,j)*hh(i,j+1)+ayt(i,j)*hh(i+1,j)+bbt(i,j)*hh(i+1,j+1)+cxt(i,j)*hh(i+2,j+1)

      if (nbg.eq.nbg0) go to 250
      j=ie(nb)
      do 240 n=1,i0
 240  rinv(ii,n,nbg0)=hh(n,j)
      if (idxy .eq. mylon) then
      if (mod(ig,i2) .eq. 0) then
      hh(i1,jg)=0.d0
      else
      hh(mod(ig,i2)+1,jg)=0.d0
      endif
      endif
 250  continue

      nsend=i2*ibir

      call mpi_gather(rinv(1,2,nbg),nsend,mpi_real8,rp(1),nsend,mpi_real8,0,mpi_comm_world,ierr)
      if (myid .eq. 0) then
      do 497 nx=1,ngx0
      ngxt=(nx-1)*ngy0
      do 497 ny=1,ngy0
      ngt=((ny-1)+ngxt)*ibir*i2
      do 497 j=1,i2
      nt=ngt+(j-1)*ibir
      jt=j+(nx-1)*i2
      do 497 i=1,ibir
      it=i+(ny-1)*ibir
 497  rtemp(it,jt)=rp(nt+i)

      do 498 i=1,i2t
      do 498 j=1,i2t
 498  rtp(j,i)=0.d0
      do 499 i=1,i2t
 499  rtp(i,i)=1.d0
      call dgesv(i2t,i2t,rtemp,i2t,ipvt,rtp,i2t,info)
      do 500 i=1,i2t
      do 500 j=1,i2t
 500  rtemp(j,i)=rtp(i,j)
      endif

      call mpi_bcast(info,1,mpi_integer,0,mpi_comm_world,ierr)
      if (info .ne. 0) then
      print*,'preporcessor failed because rinv is singular at nbg =',nbg
      stop
      endif
      nsend=i2t*i2t
      call mpi_bcast(rtp,nsend,mpi_real8,0,mpi_comm_world,ierr)
      it=mylon*i2
      jt=mylat*ibir

      jss=1
      jff=i0
      if (mylon .eq.0) then
      do 501 i=1,ibir
 501  rinv(i,1,nbg)=rtp(jt+i,i2t)
      jss=2
      endif
      if (mylon .eq. ngx1) then
      do 502 i=1,ibir
 502  rinv(i,i0,nbg)=rtp(jt+i,1)
      jff=i1
      endif
      do 503 j=jss,jff
      do 503 i=1,ibir
 503  rinv(i,j,nbg)=rtp(jt+i,it+j-1)

      if (nbg.eq.nbg0) return
      nsend=i0*ibir
      call mpi_allgather(rinv(1,1,nbg0),nsend,mpi_real8,rp,nsend,mpi_real8,m_comm_lat,ierr)
      do 504 ny=1,ngy0
      ngt=(ny-1)*ibir*i0
      do 504 j=1,i0
      nt=ngt+(j-1)*ibir
      do 504 i=1,ibir
      it=i+(ny-1)*ibir
 504  rtemp(it,j)=rp(nt+i)

      it=ibir*mylat
      do 300 i=1,ibir
      do 300 j=1,i0
      rinv1(i,j,nbg)=0.d0
      do 300 k=1,i2t
 300  rinv1(i,j,nbg)=rinv1(i,j,nbg)-rtp(it+i,k)*rtemp(k,j)
      jl=jh

 100  call mpi_barrier(mpi_comm_world,ierr)

end subroutine

! ----------------------------------------------------------------------
subroutine pre2(ax,ay,bb,cx,cy,rinv,rinv1,hh,ie)
! ----------------------------------------------------------------------
      real*8 rinv,rinv1,h,run,rtp,rtemp,rp
      dimension ax(i2,j2),ay(i2,j2),bb(i2,j2),cx(i2,j2),cy(i2,0:j2),rinv(ibir,i0,nbg0),rinv1(ibir,i0,nbg1),hh(i0,j0),ie(nb0),rp(ig0*ig0),rtemp(ig2,ig0),rtp(ig2,ig2),ipvt(ig2)
      dimension axt(i2,j2),ayt(i2,j2),bbt(i2,j2),cxt(i2,j2),cyt(i2,j2),ct(i2)

!      if (myid .eq. 0)
      do 100 ng=1,ngy0
      jl=1
      if (mylat .eq. ng-1) then
      do 99 j=1,j2
      do 99 i=1,i2
      axt(i,j)=ax(i,j)
      ayt(i,j)=ay(i,j)
      bbt(i,j)=bb(i,j)
      cxt(i,j)=cx(i,j)
 99   cyt(i,j)=cy(i,j)
      endif
      if (mylat .eq. ng-2) then
      do 98 i=1,i2
 98   ct(i)=cy(i,j2)
      endif
      nsd=i2*j2
      call mpi_bcast(axt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
      call mpi_bcast(ayt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
      call mpi_bcast(bbt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
      call mpi_bcast(cxt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
      call mpi_bcast(cyt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
      if (ng .ne. 1) call mpi_bcast(ct,i2,mpi_real8,ng-2,m_comm_lat,ierr)

      if (mylat .eq. ng-1) then
      do 97 i=1,i2
 97   cy(i,0)=ct(i)
      endif
      do 100 nb=1,nb0
      nbg=(ng-1)*nb0+nb
      write(*,111) nbg,nbg0
 111  format ('processing block #',i3,' out of ',i3,' total evp solver')
      jh=ie(nb)
      jhp=jh+1
      jhm=jh-2
      jg=jl+1
      do 250 ii=1,ibir
      ig=mylat*ibir+ii
      idxy=(ig-1)/i2
      do 210 j=jl,jhp
      do 210 i=1,i0
 210  hh(i,j)=0.d0
      if (idxy .eq. mylon) then
      if (mod(ig,i2) .eq. 0) then
      hh(i1,jg)=1.d0
      else
      hh(mod(ig,i2)+1,jg)=1.d0
      endif
      endif
      call mpi_sendrecv(hh(2,jg),1,mpi_real8,m_w,1,hh(i0,jg),1,mpi_real8,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(hh(i1,jg),1,mpi_real8,m_e,1,hh(1,jg),1,mpi_real8,m_w,1,mpi_comm_world,istat,ierr)
      if (nbg .eq. 1) go to 220
      if (idxy .eq. mylon) then
      if (nb .ne. 1) then
       if (mod(ig,i2) .eq. 0) then
       ctemp=cyt(i2,jg-2)
       else
       ctemp=cyt(mod(ig,i2),jg-2)
       endif
      else
       if (mod(ig,i2) .eq. 0) then
       ctemp=ct(i2)
       else
       ctemp=ct(mod(ig,i2))
       endif
      endif
      endif
      call mpi_bcast(ctemp,1,mpi_real8,idxy,m_comm_lon,ierr)
      do 218 n=1,i0
 218  hh(n,jl)=rinv1(ii,n,nbg-1)*ctemp
 220  continue
      do 226 j=jl,jhm
      do 225 i=1,i2
 225  hh(i+1,j+2)=-(axt(i,j)*hh(i,j+1)+ayt(i,j)*hh(i+1,j)+bbt(i,j)*hh(i+1,j+1)+cxt(i,j)*hh(i+2,j+1))/cyt(i,j)
      call mpi_sendrecv(hh(2,j+2),1,mpi_real8,m_w,1,hh(i0,j+2),1,mpi_real8,m_e,1,mpi_comm_world,istat,ierr)
 226  call mpi_sendrecv(hh(i1,j+2),1,mpi_real8,m_e,1,hh(1,j+2),1,mpi_real8,m_w,1,mpi_comm_world,istat,ierr)
      j=jh-1
      do 230 i=1,i2
 230  rinv(ii,i+1,nbg)=axt(i,j)*hh(i,j+1)+ayt(i,j)*hh(i+1,j)+bbt(i,j)*hh(i+1,j+1)+cxt(i,j)*hh(i+2,j+1)

      if (nbg.eq.nbg0) go to 250
      j=ie(nb)
      do 240 n=1,i0
 240  rinv(ii,n,nbg0)=hh(n,j)
      if (idxy .eq. mylon) then
      if (mod(ig,i2) .eq. 0) then
      hh(i1,jg)=0.d0
      else
      hh(mod(ig,i2)+1,jg)=0.d0
      endif
      endif
 250  continue

      nsend=i2*ibir

      call mpi_gather(rinv(1,2,nbg),nsend,mpi_real8,rp(1),nsend,mpi_real8,0,mpi_comm_world,ierr)
      if (myid .eq. 0) then
      do 497 nx=1,ngx0
      ngxt=(nx-1)*ngy0
      do 497 ny=1,ngy0
      ngt=((ny-1)+ngxt)*ibir*i2
      do 497 j=1,i2
      nt=ngt+(j-1)*ibir
      jt=j+(nx-1)*i2
      do 497 i=1,ibir
      it=i+(ny-1)*ibir
 497  rtemp(it,jt)=rp(nt+i)

      do 498 i=1,ig2
      do 498 j=1,ig2
 498  rtp(j,i)=0.d0
      do 499 i=1,ig2
 499  rtp(i,i)=1.d0
!      write(1233)rtemp
      call dgesv(ig2,ig2,rtemp,ig2,ipvt,rtp,ig2,info)
      do 500 i=1,ig2
      do 500 j=1,ig2
 500  rtemp(j,i)=rtp(i,j)
      endif

      call mpi_bcast(info,1,mpi_integer,0,mpi_comm_world,ierr)
      if (info .ne. 0) then
      print*,'preporcessor failed because rinv is singular at nbg =',nbg
      stop
      endif
      nsend=ig2*ig2
      call mpi_bcast(rtp,nsend,mpi_real8,0,mpi_comm_world,ierr)
      it=mylon*i2
      jt=mylat*ibir

      jss=1
      jff=i0
      if (mylon .eq.0) then
      do 501 i=1,ibir
 501  rinv(i,1,nbg)=rtp(jt+i,ig2)
      jss=2
      endif
      if (mylon .eq. ngx1) then
      do 502 i=1,ibir
 502  rinv(i,i0,nbg)=rtp(jt+i,1)
      jff=i1
      endif
      do 503 j=jss,jff
      do 503 i=1,ibir
 503  rinv(i,j,nbg)=rtp(jt+i,it+j-1)

      if (nbg.eq.nbg0) return
      nsend=i0*ibir
      call mpi_allgather(rinv(1,1,nbg0),nsend,mpi_real8,rp,nsend,mpi_real8,m_comm_lat,ierr)
      do 504 ny=1,ngy0
      ngt=(ny-1)*ibir*i0
      do 504 j=1,i0
      nt=ngt+(j-1)*ibir
      do 504 i=1,ibir
      it=i+(ny-1)*ibir
 504  rtemp(it,j)=rp(nt+i)

      it=ibir*mylat
      do 300 i=1,ibir
      do 300 j=1,i0
      rinv1(i,j,nbg)=0.d0
      do 300 k=1,ig2
 300  rinv1(i,j,nbg)=rinv1(i,j,nbg)-rtp(it+i,k)*rtemp(k,j)
      jl=jh

 100  call mpi_barrier(mpi_comm_world,ierr)

end subroutine

! ----------------------------------------------------------------------
subroutine sipinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct,alpha)
! ----------------------------------------------------------------------
      dimension al(i2,j2),ar(i2,j2),ab(i2,j2),ac(i2,j2),at(i2,0:j2),cl(i2,j2),cr(i2,j2),cb(i2,j2),cc(i2,j2),ct(i2,j2),act(i2,j2)

      do 50 j=1,j2
      do 50 i=1,i2
 50   act(i,j)=ac(i,j)

      cl(1,1)=al(1,1)
      cb(1,1)=ab(1,1)
      cc(1,1)=1/ac(1,1)
      ct(1,1)=at(1,1)*cc(1,1)
      cr(1,1)=ar(1,1)*cc(1,1)
      do 100 i=2,i2
      cl(i,1)=al(i,1)
      cb(i,1)=ab(i,1)/(1+alpha*cr(i-1,1))
      cc(i,1)=1/(act(i,1)+alpha*(cb(i,1)*cr(i-1,1))-cb(i,1)*ct(i-1,1))
      ct(i,1)=at(i,1)*cc(i,1)
 100  cr(i,1)=(ar(i,1)-alpha*cb(i,1)*ar(i-1,1))*cc(i,1)

      do 200 j=2,j2
      cl(1,j)=al(1,j)/(1+alpha*ct(1,j-1))
      cb(1,j)=ab(1,j)
      cc(1,j)=1/(act(1,j)+alpha*(cl(1,j)*ct(1,j-1))-cl(1,j)*cr(1,j-1))
      ct(1,j)=(at(1,j)-alpha*cl(1,j)*ct(1,j-1))*cc(1,j)
      cr(1,j)=ar(1,j)*cc(1,j)
      do 200 i=2,i2
      cl(i,j)=al(i,j)/(1+alpha*ct(i,j-1))
      cb(i,j)=ab(i,j)/(1+alpha*cr(i-1,j))
      cc(i,j)=1/(act(i,j)+alpha*(cl(i,j)*ct(i,j-1)+cb(i,j)*cr(i-1,j))-cl(i,j)*cr(i,j-1)-cb(i,j)*ct(i-1,j))
      ct(i,j)=(at(i,j)-alpha*cl(i,j)*ct(i,j-1))*cc(i,j)
 200  cr(i,j)=(ar(i,j)-alpha*cb(i,j)*cr(i-1,j))*cc(i,j)
end subroutine

! ----------------------------------------------------------------------
subroutine cholinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct)
! ----------------------------------------------------------------------
      dimension al(i2,j2),ar(i2,j2),ab(i2,j2),ac(i2,j2),at(i2,0:j2),cl(i2,j2),cr(i2,j2),cb(i2,j2),cc(i2,j2),ct(i2,j2),act(i2,j2)

      do 50 j=1,j2
      do 50 i=1,i2
 50   act(i,j)=ac(i,j)

      do 100 j=1,j0-3
      do 150 i=1,i0-3
      cc(i,j)=1/sqrt(act(i,j))
      cr(i,j)=ar(i,j)*cc(i,j)
      cl(i,j)=cr(i,j)
      ct(i,j)=at(i,j)*cc(i,j)
      cb(i,j)=ct(i,j)
      act(i+1,j)=act(i+1,j)-ar(i,j)*ar(i,j)*cc(i,j)*cc(i,j)
 150  act(i,j+1)=act(i,j+1)-at(i,j)*at(i,j)*cc(i,j)*cc(i,j)
      cc(i2,j)=1/sqrt(act(i2,j))
      ct(i2,j)=at(i2,j)*cc(i2,j)
      cb(i2,j)=ct(i2,j)
 100  act(i2,j+1)=act(i2,j+1)-at(i2,j)*at(i2,j)*cc(i2,j)*cc(i2,j)

      j=j2
      do 200 i=1,i3
      cc(i,j)=1/sqrt(act(i,j))
      cr(i,j)=ar(i,j)*cc(i,j)
      cl(i,j)=cr(i,j)
 200  act(i+1,j)=act(i+1,j)-ar(i,j)*ar(i,j)*cc(i,j)*cc(i,j)

      cc(i2,j2)=1/sqrt(ac(i2,j2))


end subroutine

! ----------------------------------------------------------------------
subroutine luinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct)
! ----------------------------------------------------------------------
      dimension al(i2,j2),ar(i2,j2),ab(i2,j2),ac(i2,j2),at(i2,0:j2),cl(i2,j2),cr(i2,j2),cb(i2,j2),cc(i2,j2),ct(i2,j2),act(i2,j2)

      do 50 j=1,j2
      do 50 i=1,i2
 50   act(i,j)=ac(i,j)

      do 100 j=1,j0-3
      do 150 i=1,i0-3
      cc(i,j)=1/act(i,j)
      cl(i,j)=al(i,j)*cc(i,j)
      cr(i,j)=ar(i,j)
      cb(i,j)=ab(i,j)*cc(i,j)
      ct(i,j)=at(i,j)
      act(i+1,j)=act(i+1,j)-ar(i,j)*ar(i,j)*cc(i,j)
 150  act(i,j+1)=act(i,j+1)-at(i,j)*at(i,j)*cc(i,j)
      cc(i2,j)=1/act(i2,j)
      cb(i2,j)=ab(i2,j)*cc(i2,j)
      ct(i2,j)=at(i2,j)
 100  act(i2,j+1)=act(i2,j+1)-at(i2,j)*at(i2,j)*cc(i2,j)
      j=j2
      do 200 i=1,i3
      cc(i,j)=1/act(i,j)
      cl(i,j)=al(i,j)*cc(i,j)
      cr(i,j)=ar(i,j)
 200  act(i+1,j)=act(i+1,j)-ar(i,j)*ar(i,j)*cc(i,j)

      cc(i2,j2)=1/ac(i2,j2)
end subroutine

! ---------------------------------------------------------------------
subroutine bounds
! ----------------------------------------------------------------------
      real*8 mpd,sltd,pisd,p01,p02,sb8,rb8,dg
      dimension ttemp(i0,j0,k1),stemp(i0,j0,k1)

      call mpi_type_vector(k1*j0,1,i0,mpi_real8,m_vlon,ierr) !dog
      call mpi_type_vector(k1,i0,i0*j0,mpi_real8,m_vlat,ierr)
      call mpi_type_commit(m_vlon,ierr)
      call mpi_type_commit(m_vlat,ierr)
      
      call rdnc_levitus(1)   !dog
!      read(70) s1,t1


      do 32 n=1,50

      write(*,20) n
 20   format('iteration #',i3)


      do 31 k=1,k1
      smin=100.
      smax=0.
      tmin=100.
      tmax=0.
      do 30 j=2,j1
      do 30 i=2,i1
      m=in(i-1,j,k)+in(i+1,j,k)+in(i,j-1,k)+in(i,j+1,k)
      if (m*in(i,j,k).eq.0) go to 30
      tmp=1./m
      ttemp(i,j,k)=tmp*(in(i+1,j,k)*t1(i+1,j,k)+in(i-1,j,k)*t1(i-1,j,k)+in(i,j+1,k)*t1(i,j+1,k)+in(i,j-1,k)*t1(i,j-1,k))
      stemp(i,j,k)=tmp*(in(i+1,j,k)*s1(i+1,j,k)+in(i-1,j,k)*s1(i-1,j,k)+in(i,j+1,k)*s1(i,j+1,k)+in(i,j-1,k)*s1(i,j-1,k))
      smin=min(smin,stemp(i,j,k))
      smax=max(smax,stemp(i,j,k))
      tmin=min(tmin,ttemp(i,j,k))
      tmax=max(tmax,ttemp(i,j,k))
 30   continue
 31   continue
      do 33 k=1,k1
      do 33 j=2,j1
      do 33 i=2,i1
      t1(i,j,k)=ttemp(i,j,k)
      s1(i,j,k)=stemp(i,j,k)
 33   continue      
      call mpi_sendrecv(s1(2,1,1),1,m_vlon,m_w,1,s1(i0,1,1),1,m_vlon,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(s1(i1,1,1),1,m_vlon,m_e,1,s1(1,1,1),1,m_vlon,m_w,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(s1(1,2,1),1,m_vlat,m_s,1,s1(1,j0,1),1,m_vlat,m_n,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(s1(1,j1,1),1,m_vlat,m_n,1,s1(1,1,1),1,m_vlat,m_s,1,mpi_comm_world,istat,ierr)

      call mpi_sendrecv(t1(2,1,1),1,m_vlon,m_w,1,t1(i0,1,1),1,m_vlon,m_e,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(t1(i1,1,1),1,m_vlon,m_e,1,t1(1,1,1),1,m_vlon,m_w,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(t1(1,2,1),1,m_vlat,m_s,1,t1(1,j0,1),1,m_vlat,m_n,1,mpi_comm_world,istat,ierr)
      call mpi_sendrecv(t1(1,j1,1),1,m_vlat,m_n,1,t1(1,1,1),1,m_vlat,m_s,1,mpi_comm_world,istat,ierr)
      call mpi_barrier(mpi_comm_world,ierr)
 32   continue
! levitus initial conditions have been modified
! now, we check for statically unstable points that may exist
 50   do 55 k=1,k1
      do 55 j=2,j1
      do 55 i=2,i1
      tmpd=t1(i,j,k)
      sltd=s1(i,j,k)
      pisd=.001*z(2*k)
! dan wright's full e.o.s.
!     rho(i,j,k)=r(tmp,slt,pis)
!     function r(t,s,p)
      p02=1747.4508988+tmpd*(11.51588-0.046331033*tmpd)-sltd*(3.85429655+0.01353985*tmpd)
      p01=pisd+5884.81703666+tmpd*(39.803732+tmpd*(-0.3191477+tmpd*0.0004291133))+2.6126277*sltd
 55   d1(i,j,k)=p01/(p02+0.7028423*p01)

      do 57 k=2,k1
      n=0
      do 56 j=2,j1
      do 56 i=2,i1
      if (in(i,j,k).eq.0.or.d1(i,j,k).ge.d1(i,j,k-1)) go to 56
      n=n+1
      t1(i,j,k)=t1(i,j,k-1)
      s1(i,j,k)=s1(i,j,k-1)
      tmpd=t1(i,j,k)
      sltd=s1(i,j,k)
      pisd=.001*z(2*k)
      p02=1747.4508988+tmpd*(11.51588-0.046331033*tmpd)-sltd*(3.85429655+0.01353985*tmpd)
      p01=pisd+5884.81703666+tmpd*(39.803732+tmpd*(-0.3191477+tmpd*0.0004291133))+2.6126277*sltd
      d1(i,j,k)=p01/(p02+0.7028423*p01)
 56   continue
 57   write(*,58) n,k
 58   format(i6,' unstable points at level',i3)

! periodic
! remark that there doesn't not exchange the b.c. the b! will be exchange in fsglo

! write climatology based initial conditions
! climatology is modified in prep program

! wenien
!      open(69,file='data/initial'//grd,form='unformatted')
!      write(69) s1,t1
!      close(69)
      call wr8nc_preini  !dog

! 199  rewind 70

!      if (.false.) then
! seasonal boundary conditions


      do 7000 nm=1,12
      write(*,238) nm
 238  format('month',i3,' in bounds')

! read in monthly t,s boundary climatologies
      call rdnc_levitus(nm) !dog
!      read(70) s1,t1

      smin=100.
      smax=0.
      tmin=100.
      tmax=0.

      do 2400 k=1,k1
      do 2400 j=2,j1
      do 2400 i=2,i1
      s1(i,j,k)=min(max(s1(i,j,k),smax),smin)
 2400 t1(i,j,k)=min(max(t1(i,j,k),tmax),tmin)

      do 240 j=2,j1
      do 240 i=2,i1
      tsurf(i-1,j-1,nm)=t1(i,j,1)
 240  ssurf(i-1,j-1,nm)=s1(i,j,1)


      do 260 k=1,k1
      avg1=0.
      avg2=0.
      areaxy=0.
      smax=0.
      smin=100.
      do 250 j=2,j1
      tmp=dx(j)*dy(j)

      do 250 i=2,i1
      temp=s1(i,j,k)
      s1(i,j,k)=max(s1(i,j,k),25.)
      if (in(i,j,k).eq.0) go to 250
      smax=max(smax,s1(i,j,k))
      smin=min(smin,s1(i,j,k))

      if (s1(i,j,k).eq.25.) write(*,29) i,j,k,temp,s1(i,j,k)
 29   format('fixed bad point i,j,k,s1bad,s1fix=',3i5,2f5.1)

      areaxy=areaxy+tmp
      avg1=avg1+tmp*t1(i,j,k)
      avg2=avg2+tmp*s1(i,j,k)
 250  continue

      tp=avg1
      call mpi_allreduce(tp,avg1,1,mpi_real8,mpi_sum,mpi_comm_world,ierr)
      tp=avg2
      call mpi_allreduce(tp,avg2,1,mpi_real8,mpi_sum,mpi_comm_world,ierr)
      tp=areaxy
      call mpi_allreduce(tp,areaxy,1,mpi_real8,mpi_sum,mpi_comm_world,ierr)


      write(*,251) nm,k,smin,smax,avg2/areaxy
 251  format('nm,k,smin,smax,savg=',2i3,1x,3f5.1)

!      if(areaxy.eq.0.) print*,mylon,'_',mylat,'k=',k,'areaxy=',areaxy

      txycli(k,nm)=avg1/areaxy
 260  sxycli(k,nm)=avg2/areaxy

      write(*,261) nm,(sxycli(k,nm),k=1,k1)
 261  format(i2/(21f5.1))

! sponge layer climatology
! k=1 is included in ssurf, tssurf
! le
! generally, j0 will not less than 10, so there requires to fix tssp etc.,
! if j0 is less than 10.
      if (mylat .eq. 0) then
!      do 270 k=1,k1
!      do 270 j=1,j2t/18
      do 270 k=2,k1
      do 270 j=1,10
      do 270 i=2,i1
!      tssp(i-1,j,k,nm)=t1(i,j+1,k)
! 270  sssp(i-1,j,k,nm)=s1(i,j+1,k)
      tssp(i-1,j,k-1,nm)=t1(i,j+1,k)
 270  sssp(i-1,j,k-1,nm)=s1(i,j+1,k) 
      endif

      if (mylat .eq. ngy0-1) then
!      do 280 k=1,k1
!      do 280 j=1,j2t/18
      do 280 k=2,k1
      do 280 j=1,10
      do 280 i=2,i1
!      tnsp(i-1,j,k,nm)=t1(i,j0-j,k)
! 280  snsp(i-1,j,k,nm)=s1(i,j0-j,k)
      tnsp(i-1,j,k-1,nm)=t1(i,j0-j,k)
 280  snsp(i-1,j,k-1,nm)=s1(i,j0-j,k)
      endif

 7000 continue


! although the following is superfluous when using seasonal climatology,
! (because the model uses linear interpolations anyway),
! we will need 12-month arrays when we use the monthly climatology

! interpolate tsurf, ssurf in time
      tmp=1./3.
      tem=2./3.
      do 600 n=1,10,3
      nm=n+3
      if (nm.eq.13) nm=1
      do 600 j=1,j2
      do 600 i=1,i2
      tsurf(i,j,n+1)=tsurf(i,j,n)+tmp*(tsurf(i,j,nm)-tsurf(i,j,n))
      tsurf(i,j,n+2)=tsurf(i,j,n)+tem*(tsurf(i,j,nm)-tsurf(i,j,n))
      ssurf(i,j,n+1)=ssurf(i,j,n)+tmp*(ssurf(i,j,nm)-ssurf(i,j,n))
 600  ssurf(i,j,n+2)=ssurf(i,j,n)+tem*(ssurf(i,j,nm)-ssurf(i,j,n))

! interpolate seasonal sponge climatologies
      do 650 n=1,10,3
      nm=n+3
      if (nm.eq.13) nm=1
      do 648 k=1,k1
      txycli(k,n+1)=txycli(k,n)+tmp*(txycli(k,nm)-txycli(k,n))
      txycli(k,n+2)=txycli(k,n)+tem*(txycli(k,nm)-txycli(k,n))
      sxycli(k,n+1)=sxycli(k,n)+tmp*(sxycli(k,nm)-sxycli(k,n))
 648  sxycli(k,n+2)=sxycli(k,n)+tem*(sxycli(k,nm)-sxycli(k,n))
!      do 650 k=1,k1
!      do 650 j=1,j2t/18
      do 650 k=1,k2
      do 650 j=1,10
      do 650 i=1,i2
      tssp(i,j,k,n+1)=tssp(i,j,k,n)+tmp*(tssp(i,j,k,nm)-tssp(i,j,k,n))
      tssp(i,j,k,n+2)=tssp(i,j,k,n)+tem*(tssp(i,j,k,nm)-tssp(i,j,k,n))
      tnsp(i,j,k,n+1)=tnsp(i,j,k,n)+tmp*(tnsp(i,j,k,nm)-tnsp(i,j,k,n))
      tnsp(i,j,k,n+2)=tnsp(i,j,k,n)+tem*(tnsp(i,j,k,nm)-tnsp(i,j,k,n))
      sssp(i,j,k,n+1)=sssp(i,j,k,n)+tmp*(sssp(i,j,k,nm)-sssp(i,j,k,n))
      sssp(i,j,k,n+2)=sssp(i,j,k,n)+tem*(sssp(i,j,k,nm)-sssp(i,j,k,n))
      snsp(i,j,k,n+1)=snsp(i,j,k,n)+tmp*(snsp(i,j,k,nm)-snsp(i,j,k,n))
 650  snsp(i,j,k,n+2)=snsp(i,j,k,n)+tem*(snsp(i,j,k,nm)-snsp(i,j,k,n))

! write surface climatology t,s and specified inflow velocities
!      write(62) ssurf,tsurf,sssp,snsp,tssp,tnsp
!      write(62) sxycli,txycli
      call wr8nc_boundaries !dog
end subroutine
! ----------------------------------------------------------------------

subroutine initfs
! initfs initializes all controlling arrays and derived scalars
!(scalar control parameters are set in block data at start of this file)
! ----------------------------------------------------------------------

      dt=172800./daodt
      odt=1./dt
      orzmx=1./rzmx
      ofltw=1.-2.*fltw

      mvi=mvfreq*daodt
      prn=dm0/de0

! ====================================================
! set background vertical eddy viscosity & diffusivity
! ====================================================
      do 90 k=1,k2
! augment molecular viscosity & diffusivity by parameterized synopti! wind
! events, and by breaking (u,v,t,s, near surface only) and non-breaking
! (u,v only, at all rmetss) internal waves.
! synopti! wind forced mixing having 20m e-folding scale.
! bigger momentum mixing due to pressure-xfers
! associated with internal waves that have no counterpart in t,s xfers.
! after yr 7
! 10m scale height for t,s mixing due to breaking internal waves
       tmpts=exp(-.001*(z(2*k+1)-z(3)))
!
! deeper tmpuv enhancement is motivated by findings of hurlburt & hogan
! that deep bathymetry and deep transient eddies influences gs path.
! tmpts is not enhanced (layer model has none across layer interface!)
! 1000m scale height for u,v mixing (due to pressure-xfers even without
! internal wave breaking). this enhances deep bathymetry effe!cts on
! upper level flow (e.g., gs separation, kuroshio meanders)
       tmpuv=exp(-.00001*(z(2*k+1)-z(3)))
      vbk(k)=.01+.05*tmpuv
 90   hbk(k)=.002+.02*tmpts
!      write(8) vbk,hbk
!      close(8)
      call wr8nc_winmx    !dog

      if(myid.eq.0)then
      write(*,91) (k,int(.01*z(2*k+1)),vbk(k),hbk(k),k=1,k2)
!      write(7,91) (k,int(.01*z(2*k+1)),vbk(k),hbk(k),k=1,k2)
 91   format(' K    Z    VBK    HBK'/(i2,i5,2f7.3))
      endif

! wenien
!      write(9) case,dscrib,dt,tlz,b,g,dm0,de0,dmz0,f,tanphi,rzmx,drag,fltw,daodt,odt,prn,orzmx,ofltw,y,yv,yvdeg,ydeg,cs,csv,ocs,ocsv,dx,dxv,odx,odxv,dy,dyv,ody,odyv,ktrm,lrstrt,mxit,isav,mxsav,lopen,lmovi,mvi,lfsrf,lsolver,xdeg,lprecon,rawf

!      close(9)
      call wr8nc_rundata !dog
      call mpi_barrier(mpi_comm_world,ierr)

      call bounds

      if (lsolver.eq.0) then
      print*,'preprocessor running without solver setting has finished.'
      stop 0
      endif

      do 250 n=1,nb0
 250  ie(n)=min(1+n0*n,j1)
!      write(*,251) ie
! 251  format('ie'/(20i4))
      if (ie(nb0).gt.ie(nb1).and.j1.lt.ie(nb1)+9) go to 260
!      write(17,252) ie(nb1),j1
      write(*,252) ie(nb1),j1
 252  format('ie(nb1),j1=',2i5,' not properly defined. Fix above do 250 loop and restart.')

 260  ie(nb0)=j1


! preprocessor for ellipti! solver used to maintain
! incompressibility for rigid-lid hydrostati! flows
! first, do top layer (all water to regularize evp domain)
      tmp=1./odz(1)
      do 280 j=1,j2
      dzx=odx(j+1)**2*tmp
      dzyb=csv(j)*ocs(j+1)*odyv(j)*ody(j+1)*tmp
      dzyt=csv(j+1)*ocs(j+1)*odyv(j+1)*ody(j+1)*tmp
      do 280 i=1,i2
      al(i,j)=dzx
      ar(i,j)=dzx
      ab(i,j)=dzyb
 280  at(i,j)=dzyt
! remaining layers
      do 282 k=2,k1
      tmp=1./odz(k)
      do 282 j=1,j2
      dzx=odx(j+1)**2*tmp
      dzyb=csv(j)*ocs(j+1)*odyv(j)*ody(j+1)*tmp
      dzyt=csv(j+1)*ocs(j+1)*odyv(j+1)*ody(j+1)*tmp
      do 282 i=1,i2
      al(i,j)=al(i,j)+dzx*iu(i,j+1,k)
      ar(i,j)=ar(i,j)+dzx*iu(i+1,j+1,k)
      ab(i,j)=ab(i,j)+dzyb*iv(i+1,j,k)
 282  at(i,j)=at(i,j)+dzyt*iv(i+1,j+1,k)

      if (mylat .eq. 0) then
      do 285 i=1,i2
 285  ab(i,1)=0.
      endif
      if (mylat .eq. ngy1) then
      do 286 i=1,i2
 286  at(i,j2)=0.
      endif

      gamma_tmp=0.
      if (lfsrf .eq. 1) gamma_tmp=1./(0.5*dt**2)
      do 287 j=1,j2
      do 287 i=1,i2
      al(i,j)=-al(i,j)
      ar(i,j)=-ar(i,j)
      ab(i,j)=-ab(i,j)
      at(i,j)=-at(i,j)
 287  ac(i,j)=(-al(i,j)-ar(i,j)-ab(i,j)-at(i,j))+gamma_tmp

      call mpi_barrier(mpi_comm_world,ierr)

      if(lsolver .eq.1) then
        if(myid.eq.0) print*,'enter evp preprocessor',myid
        call pre(al,ab,ac,ar,at,rinv,rinv1,hh,ie)
!        write(99) al,ar,ab,at,ac,rinv,rinv1,ie
!        close(99)
        call wr8nc_evp1 !dog
        if(myid.eq.0) print*,'finished evp preprocessor',myid
      elseif(lsolver .eq.2) then
        if(myid.eq.0) print*,'enter preconditioner',myid
        if(lprecon .le. 1) then
          if(myid.eq.0) print*,'sip precondition'
          alpha=.5d0
	  call sipinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct,alpha)
        elseif(lprecon .eq. 2) then
          if(myid.eq.0) print*,'lu precondition'
          call luinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct)
        elseif(lprecon .eq. 3) then
          if(myid.eq.0) print*,'chol precondition'
          call cholinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct)
        else
          if(myid.eq.0) print*,'lprecon=',lprecon,' err'
          stop 0
        endif
!	write(99) al,ar,ab,at,ac,cl,cb,cc,cr,ct
!        close(99)
        call wr8nc_evp2 !dog
      endif    
  
      print*,'finished preconditioner',myid

end subroutine

subroutine wr8nc_evp1
use const, only: i2t,j2t,i0t,ibir,ngy0,nbg0,nbg1,nb0
use comm,  only: myid,mylon,mylat,m_cart
implicit none
!declare pnetcdf var
character :: fn*11,desp*30
integer :: ierr,ncid,nt(3),nd(2),ns(1)
integer :: i2id,j2id,i0id,ibid,ibg0,ibg1,j2idd,inb0
integer :: ial,iar,iab,iac,iat,irin,irin1,ieid
integer(kind=mpi_offset_kind) :: dim_nn(7),leng
integer(kind=mpi_offset_kind) :: ns_ijk(3),nc_ijk(3),ns_ij(2),nc_ij(2),ns_i,nc_i
      !derive the size of each mpiid task,320,288

      !output file name
      fn = 'data/evp.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) = i0*ngx0
      dim_nn(4) = ibir*ngy0
      dim_nn(5) = nbg0
      dim_nn(6) = nbg1
      dim_nn(7) = nb0
      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,"i0_global",dim_nn(3),i0id)
      ierr=nfmpi_def_dim(ncid,"j2+1_global",dim_nn(2)+ngy0,j2idd)
      ierr=nfmpi_def_dim(ncid,"ibir_global",dim_nn(4),ibid)
      ierr=nfmpi_def_dim(ncid,"nbg0",dim_nn(5),ibg0)
      ierr=nfmpi_def_dim(ncid,"nbg1",dim_nn(6),ibg1)
      ierr=nfmpi_def_dim(ncid,"nb0",dim_nn(7),inb0)

      !define variables and assign varID
      nd(1) = j2id
      nd(2) = i2id
      ierr=nfmpi_def_var(ncid,"al",nf_double,2,nd,ial)
      ierr=nfmpi_def_var(ncid,"ar",nf_double,2,nd,iar)
      ierr=nfmpi_def_var(ncid,"ab",nf_double,2,nd,iab)
      ierr=nfmpi_def_var(ncid,"ac",nf_double,2,nd,iac)

      nd(1) = i2id
      nd(2) = j2idd
      ierr=nfmpi_def_var(ncid,"at",nf_double,2,nd,iat)

      nt(1) = ibid
      nt(2) = i0id
      nt(3) = ibg0
      ierr=nfmpi_def_var(ncid,"rinv",nf_double,3,nt,irin)

      nt(1) = ibid
      nt(2) = i0id
      nt(3) = ibg1
      ierr=nfmpi_def_var(ncid,"rinv1",nf_double,3,nt,irin1)

      ns(1) = inb0
      ierr=nfmpi_def_var(ncid,"ie",nf_int,1,ns,ieid)

      desp='input files for lsolver = 1'
      leng=len(trim(desp))
      ierr=nfmpi_put_att_text(ncid,NF_GLOBAL,'desp',leng,trim(desp))
      ierr=nfmpi_enddef(ncid)

      !put variables
      ns_ij(1) = mylon*i2+1
      ns_ij(2) = mylat*j2+1
      nc_ij(1) = i2
      nc_ij(2) = j2
      ierr=nfmpi_put_vara_double_all(ncid,ial,ns_ij,nc_ij,al)
      ierr=nfmpi_put_vara_double_all(ncid,iar,ns_ij,nc_ij,ar)
      ierr=nfmpi_put_vara_double_all(ncid,iab,ns_ij,nc_ij,ab)
      ierr=nfmpi_put_vara_double_all(ncid,iac,ns_ij,nc_ij,ac)
      ns_ij(1) = mylon*j2+1
      ns_ij(2) = mylat*(j2+1)+1
      nc_ij(1) = i2
      nc_ij(2) = j2+1
      ierr=nfmpi_put_vara_double_all(ncid,iat,ns_ij,nc_ij,at)

      ns_ijk(1) = mylat*ibir+1
      ns_ijk(2) = mylon*i0+1
      ns_ijk(3) = 1
      nc_ijk(1) = ibir
      nc_ijk(2) = i0
      nc_ijk(3) = nbg0
      ierr=nfmpi_put_vara_double_all(ncid,irin,ns_ijk,nc_ijk,rinv)
      ns_ijk(1) = mylat*ibir+1
      ns_ijk(2) = mylon*i0+1
      ns_ijk(3) = 1
      nc_ijk(1) = ibir
      nc_ijk(2) = i0
      nc_ijk(3) = nbg1
      ierr=nfmpi_put_vara_double_all(ncid,irin1,ns_ijk,nc_ijk,rinv1)

      ns_i = 1
      nc_i = nb0
      ierr=nfmpi_put_vara_int_all(ncid,ieid,ns_i,nc_i,ie)
      ierr=nfmpi_close(ncid)
end subroutine wr8nc_evp1

subroutine wr8nc_evp2
use const, only: i2,j2,ngy0,ngx0
use comm,  only: myid,mylon,mylat,m_cart
use var,   only: grd
implicit none
!declare pnetcdf var
character :: fn*11,desp*120
integer :: ierr,ncid,nd(2)
integer :: i2id,j2id,j2idd
integer :: ial,iar,iab,iac,iat,icl,icb,icc,icr,ict
integer(kind=mpi_offset_kind) :: dim_nn(3),leng
integer(kind=mpi_offset_kind) :: ns_ij(2),nc_ij(2)
      !derive the size of each mpiid task

      !output file name
      fn = 'data/evp.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) = (j2+1)*ngy0
      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,"j2+1_global",dim_nn(3),j2idd)
      !define variables and assign varID
      nd(1) = i2id
      nd(2) = j2id
      ierr=nfmpi_def_var(ncid,"al",nf_double,2,nd,ial)
      ierr=nfmpi_def_var(ncid,"ar",nf_double,2,nd,iar)
      ierr=nfmpi_def_var(ncid,"ab",nf_double,2,nd,iab)
      ierr=nfmpi_def_var(ncid,"ac",nf_double,2,nd,iac)

      nd(2) = j2idd
      ierr=nfmpi_def_var(ncid,"at",nf_double,2,nd,iat)

      nd(2) = j2id
      ierr=nfmpi_def_var(ncid,"cl",nf_double,2,nd,icl)
      ierr=nfmpi_def_var(ncid,"cb",nf_double,2,nd,icb)
      ierr=nfmpi_def_var(ncid,"cc",nf_double,2,nd,icc)
      ierr=nfmpi_def_var(ncid,"cr",nf_double,2,nd,icr)
      ierr=nfmpi_def_var(ncid,"ct",nf_double,2,nd,ict)

      desp='input files for lsolver = 2, note that the global domain &
            is decomposed into ngx0*ngy0 subdomains'
      leng=len(trim(desp))
      ierr=nfmpi_put_att_text(ncid,NF_GLOBAL,'desp',leng,trim(desp))
      ierr=nfmpi_enddef(ncid)

      !put variables
      ns_ij(1) = mylon*i2+1
      ns_ij(2) = mylat*j2+1
      nc_ij(1) = i2
      nc_ij(2) = j2
      ierr=nfmpi_put_vara_double_all(ncid,ial,ns_ij,nc_ij,al)
      ierr=nfmpi_put_vara_double_all(ncid,iar,ns_ij,nc_ij,ar)
      ierr=nfmpi_put_vara_double_all(ncid,iab,ns_ij,nc_ij,ab)
      ierr=nfmpi_put_vara_double_all(ncid,iac,ns_ij,nc_ij,ac)
      ns_ij(1) = mylon*i2+1
      ns_ij(2) = mylat*(j2+1)+1
      nc_ij(1) = i2
      nc_ij(2) = j2+1
      ierr=nfmpi_put_vara_double_all(ncid,iat,ns_ij,nc_ij,at)

      ns_ij(1) = mylon*i2+1
      ns_ij(2) = mylat*j2+1
      nc_ij(1) = i2
      nc_ij(2) = j2
      ierr=nfmpi_put_vara_double_all(ncid,icl,ns_ij,nc_ij,cl)
      ierr=nfmpi_put_vara_double_all(ncid,icb,ns_ij,nc_ij,cb)
      ierr=nfmpi_put_vara_double_all(ncid,icc,ns_ij,nc_ij,cc)
      ierr=nfmpi_put_vara_double_all(ncid,icr,ns_ij,nc_ij,cr)
      ierr=nfmpi_put_vara_double_all(ncid,ict,ns_ij,nc_ij,ct)

      ierr=nfmpi_close(ncid)
end subroutine wr8nc_evp2

subroutine wr8nc_winmx
use const,  only: k2
use var,    only: z
use comm,  only: m_cart
implicit none

!declare pnetcdf var
character :: fn*15,desp*100
integer :: ierr,ncid,ns
integer :: k2id,i
integer :: ivbk,ihbk,idep
integer(kind=mpi_offset_kind) :: dim_nn(1),leng
integer(kind=mpi_offset_kind) :: sns,cns
integer :: dep(k2)
      !derive the size of each mpiid task
      sns = 1
      cns = k2

      do i = 1,k2
      dep(i) = int(.01*z(2*i+1))
      end do

      !output file name
      fn = 'data/windmix.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) = k2
      ierr=nfmpi_def_dim(ncid,"k2",dim_nn(1),k2id)

      !define variables and assign varID
      ns = k2id
      ierr=nfmpi_def_var(ncid,"depth",nf_int,1,ns,idep)
      ierr=nfmpi_def_var(ncid,"vbk",nf_double,1,ns,ivbk)
      ierr=nfmpi_def_var(ncid,"hbk",nf_double,1,ns,ihbk)

      desp='input files for vertical eddy viscosity and diffusivity '
      leng=len(trim(desp))
      ierr=nfmpi_put_att_text(ncid,NF_GLOBAL,'desp',leng,trim(desp))

      ierr=nfmpi_enddef(ncid)

      !put variables
      ierr=nfmpi_put_vara_int_all(ncid,idep,sns,cns,dep)
      ierr=nfmpi_put_vara_double_all(ncid,ivbk,sns,cns,vbk)
      ierr=nfmpi_put_vara_double_all(ncid,ihbk,sns,cns,hbk)

      ierr=nfmpi_close(ncid)
end subroutine wr8nc_winmx

subroutine wr8nc_rundata
use const, only: j2,j1,j0,i0,j2t,j1t,j0t,i0t,ngx0,ngy0
use comm,  only: myid,mylat,mylon,m_cart
implicit none

!declare pnetcdf var
character :: fn*15
integer :: ncid,ierr,ns
integer :: siid,i0id,j0id,j1id,j2id,cid,dsid
integer :: icase,idsp,idt,itlz,ib,ig,idm0,ide0,idmz0
integer :: irzmx,idrag,ifltw,idadt,iodt,iprn,iorz,ioflt,irawf
integer :: iktrm,irst,imxit,isv,imxsv,ilop,ilmov,imvi,ilfs,ilsol,ilpre
integer :: ixdeg,ifcol,iphi,icsv,iocs,idxv,iodx,idyv,iody,iydeg
integer :: iy,iyv,iyvd,ics,iocs2,idx,iodxx,idy,iodyy
integer(kind=mpi_offset_kind) :: dim_nn(5),leng
integer(kind=mpi_offset_kind) :: ns_i,nc_i
      !output file name
      fn = "data/rundata.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) = 1
      dim_nn(2) = i0*ngx0
      dim_nn(3) = j0*ngy0
      dim_nn(4) = j1*ngy0
      dim_nn(5) = j2*ngy0
      ierr=nfmpi_def_dim(ncid,"char_case",dim_nn(1)*5,cid)
      ierr=nfmpi_def_dim(ncid,"char_desp",dim_nn(1)*63,dsid)
      ierr=nfmpi_def_dim(ncid,"pt_val",dim_nn(1),siid)
      ierr=nfmpi_def_dim(ncid,"i0_global",dim_nn(2),i0id)
      ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn(3),j0id)
      ierr=nfmpi_def_dim(ncid,"j1_global",dim_nn(4),j1id)
      ierr=nfmpi_def_dim(ncid,"j2_global",dim_nn(5),j2id)
      !define variables and assign varID
      ns = cid
      ierr=nfmpi_def_var(ncid,"case",nf_char,1,ns,icase)
      ns = dsid
      ierr=nfmpi_def_var(ncid,"describ",nf_char,1,ns,idsp)
      ns = siid
      ierr=nfmpi_def_var(ncid,"dt",nf_double,1,ns,idt)
      ierr=nfmpi_def_var(ncid,"tlz",nf_double,1,ns,itlz)
      ierr=nfmpi_def_var(ncid,"b",nf_double,1,ns,ib)
      ierr=nfmpi_def_var(ncid,"g",nf_double,1,ns,ig)
      ierr=nfmpi_def_var(ncid,"dm0",nf_double,1,ns,idm0)
      ierr=nfmpi_def_var(ncid,"de0",nf_double,1,ns,ide0)
      ierr=nfmpi_def_var(ncid,"dmz0",nf_double,1,ns,idmz0)
      ierr=nfmpi_def_var(ncid,"rzmx",nf_double,1,ns,irzmx)
      ierr=nfmpi_def_var(ncid,"drag",nf_double,1,ns,idrag)
      ierr=nfmpi_def_var(ncid,"fltw",nf_double,1,ns,ifltw)
      ierr=nfmpi_def_var(ncid,"daodt",nf_double,1,ns,idadt)
      ierr=nfmpi_def_var(ncid,"odt",nf_double,1,ns,iodt)
      ierr=nfmpi_def_var(ncid,"prn",nf_double,1,ns,iprn)
      ierr=nfmpi_def_var(ncid,"orzmx",nf_double,1,ns,iorz)
      ierr=nfmpi_def_var(ncid,"ofltw",nf_double,1,ns,ioflt)
      ierr=nfmpi_def_var(ncid,"rawf",nf_double,1,ns,irawf)
      ierr=nfmpi_def_var(ncid,"ktrm",nf_int,1,ns,iktrm)
      ierr=nfmpi_def_var(ncid,"lrstrt",nf_int,1,ns,irst)
      ierr=nfmpi_def_var(ncid,"mxit",nf_int,1,ns,imxit)
      ierr=nfmpi_def_var(ncid,"isav",nf_int,1,ns,isv)
      ierr=nfmpi_def_var(ncid,"mxsav",nf_int,1,ns,imxsv)
      ierr=nfmpi_def_var(ncid,"lopen",nf_int,1,ns,ilop)
      ierr=nfmpi_def_var(ncid,"lmovi",nf_int,1,ns,ilmov)
      ierr=nfmpi_def_var(ncid,"mvi",nf_int,1,ns,imvi)
      ierr=nfmpi_def_var(ncid,"lfsrf",nf_int,1,ns,ilfs)
      ierr=nfmpi_def_var(ncid,"lsolver",nf_int,1,ns,ilsol)
      ierr=nfmpi_def_var(ncid,"lprecon",nf_int,1,ns,ilpre)

      ns = j2id
      ierr=nfmpi_def_var(ncid,"f",nf_double,1,ns,ifcol)
      ierr=nfmpi_def_var(ncid,"tanphi",nf_double,1,ns,iphi)
      ns = j1id ! csv,ocsv,dxv,odxv,dyv,odyv
      ierr=nfmpi_def_var(ncid,"csv",nf_double,1,ns,icsv)
      ierr=nfmpi_def_var(ncid,"ocsv",nf_double,1,ns,iocs)
      ierr=nfmpi_def_var(ncid,"dxv",nf_double,1,ns,idxv)
      ierr=nfmpi_def_var(ncid,"odxv",nf_double,1,ns,iodx)
      ierr=nfmpi_def_var(ncid,"dyv",nf_double,1,ns,idyv)
      ierr=nfmpi_def_var(ncid,"odyv",nf_double,1,ns,iody)
      ns = j0id ! y,yv,yvdeg,cs,ocs,dx,odx,dy,ody
      ierr=nfmpi_def_var(ncid,"y",nf_double,1,ns,iy)
      ierr=nfmpi_def_var(ncid,"yv",nf_double,1,ns,iyv)
      ierr=nfmpi_def_var(ncid,"yvdeg",nf_double,1,ns,iyvd)
      ierr=nfmpi_def_var(ncid,"ydeg",nf_double,1,ns,iydeg)
      ierr=nfmpi_def_var(ncid,"cs",nf_double,1,ns,ics)
      ierr=nfmpi_def_var(ncid,"ocs",nf_double,1,ns,iocs2)
      ierr=nfmpi_def_var(ncid,"dx",nf_double,1,ns,idx)
      ierr=nfmpi_def_var(ncid,"odx",nf_double,1,ns,iodxx)
      ierr=nfmpi_def_var(ncid,"dy",nf_double,1,ns,idy)
      ierr=nfmpi_def_var(ncid,"ody",nf_double,1,ns,iodyy)
      ns = i0id
      ierr=nfmpi_def_var(ncid,"xdeg",nf_double,1,ns,ixdeg)
      ierr=nfmpi_enddef(ncid)
      !put variables in collective data mode
      ns_i = mylat*j2+1
      nc_i = j2
      ierr=nfmpi_put_vara_double_all(ncid,ifcol,ns_i,nc_i,f)
      ierr=nfmpi_put_vara_double_all(ncid,iphi,ns_i,nc_i,tanphi)
      ns_i = mylat*j1+1
      nc_i = j1
      ierr=nfmpi_put_vara_double_all(ncid,icsv,ns_i,nc_i,csv)
      ierr=nfmpi_put_vara_double_all(ncid,iocs,ns_i,nc_i,ocsv)
      ierr=nfmpi_put_vara_double_all(ncid,idxv,ns_i,nc_i,dxv)
      ierr=nfmpi_put_vara_double_all(ncid,iodx,ns_i,nc_i,odxv)
      ierr=nfmpi_put_vara_double_all(ncid,idyv,ns_i,nc_i,dyv)
      ierr=nfmpi_put_vara_double_all(ncid,iody,ns_i,nc_i,odyv)
      ns_i = mylat*j0+1
      nc_i = j0
      ierr=nfmpi_put_vara_double_all(ncid,iy,ns_i,nc_i,y)
      ierr=nfmpi_put_vara_double_all(ncid,iyv,ns_i,nc_i,yv)
      ierr=nfmpi_put_vara_double_all(ncid,iydeg,ns_i,nc_i,ydeg)
      ierr=nfmpi_put_vara_double_all(ncid,iyvd,ns_i,nc_i,yvdeg)
      ierr=nfmpi_put_vara_double_all(ncid,ics,ns_i,nc_i,cs)
      ierr=nfmpi_put_vara_double_all(ncid,iocs2,ns_i,nc_i,ocs)
      ierr=nfmpi_put_vara_double_all(ncid,idx,ns_i,nc_i,dx)
      ierr=nfmpi_put_vara_double_all(ncid,iodxx,ns_i,nc_i,odx)
      ierr=nfmpi_put_vara_double_all(ncid,idy,ns_i,nc_i,dy)
      ierr=nfmpi_put_vara_double_all(ncid,iodyy,ns_i,nc_i,ody)
      ns_i = mylon*i0+1
      nc_i = i0
      ierr=nfmpi_put_vara_double_all(ncid,ixdeg,ns_i,nc_i,xdeg)
      !put variables in independent data mode
      ierr=nfmpi_begin_indep_data(ncid)
      ierr=nfmpi_put_var_text(ncid,icase,case)
      ierr=nfmpi_put_var_text(ncid,idsp,dscrib)
      ierr=nfmpi_put_var_double(ncid,idt,dt)
      ierr=nfmpi_put_var_double(ncid,itlz,tlz)
      ierr=nfmpi_put_var_double(ncid,ib,b)
      ierr=nfmpi_put_var_double(ncid,ig,g)
      ierr=nfmpi_put_var_double(ncid,idm0,dm0)
      ierr=nfmpi_put_var_double(ncid,ide0,de0)
      ierr=nfmpi_put_var_double(ncid,idmz0,dmz0)
      ierr=nfmpi_put_var_double(ncid,irzmx,rzmx)
      ierr=nfmpi_put_var_double(ncid,idrag,drag)
      ierr=nfmpi_put_var_double(ncid,ifltw,fltw)
      ierr=nfmpi_put_var_double(ncid,idadt,daodt)
      ierr=nfmpi_put_var_double(ncid,iodt,odt)
      ierr=nfmpi_put_var_double(ncid,iprn,prn)
      ierr=nfmpi_put_var_double(ncid,iorz,orzmx)
      ierr=nfmpi_put_var_double(ncid,ioflt,ofltw)
      ierr=nfmpi_put_var_double(ncid,irawf,rawf)
      ierr=nfmpi_put_var_int(ncid,iktrm,ktrm)
      ierr=nfmpi_put_var_int(ncid,irst,lrstrt)
      ierr=nfmpi_put_var_int(ncid,imxit,mxit)
      ierr=nfmpi_put_var_int(ncid,isv,isav)
      ierr=nfmpi_put_var_int(ncid,imxsv,mxsav)
      ierr=nfmpi_put_var_int(ncid,ilop,lopen)
      ierr=nfmpi_put_var_int(ncid,ilmov,lmovi)
      ierr=nfmpi_put_var_int(ncid,imvi,mvi)
      ierr=nfmpi_put_var_int(ncid,ilfs,lfsrf)
      ierr=nfmpi_put_var_int(ncid,ilsol,lsolver)
      ierr=nfmpi_put_var_int(ncid,ilpre,lprecon)
      ierr=nfmpi_end_indep_data(ncid)

      ierr=nfmpi_close(ncid)
end subroutine wr8nc_rundata

subroutine wr8nc_zkb
use comm, only: myid,mylat,mylon,m_cart
use const, only: i0,j0,k0,i01,k0,k1,ngx0,ngy0
implicit none
!declare pnetcdf var
character :: fn*11
integer :: ncid,ierr,ns,nd(2),nt(3)
integer :: siid,k0id,k1id,kkid,i01id,i0id,j0id,i11id,j11id
integer :: itlz,iz,iodz,iodzw,ikb,irmt,iiv0,iiu0
integer :: iiu,iiv,iiw,iin
integer(kind=mpi_offset_kind) :: dim_nn(9),ns_ijk(3),nc_ijk(3),ns_ij(2),nc_ij(2)
integer(kind=mpi_offset_kind) :: ns_i,nc_i

      !output file name
      fn = 'data/zkb.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) = 1
      dim_nn(2) = k0
      dim_nn(3) = k1
      dim_nn(4) = k0+k1
      dim_nn(5) = i01*ngx0
      dim_nn(6) = i0*ngx0
      dim_nn(7) = j0*ngy0
      dim_nn(8) = (i01+1)*ngx0
      dim_nn(9) = (j0+1)*ngy0
      ierr=nfmpi_def_dim(ncid,"pt_val",dim_nn(1),siid)
      ierr=nfmpi_def_dim(ncid,"k0",dim_nn(2),k0id)
      ierr=nfmpi_def_dim(ncid,"k1",dim_nn(3),k1id)
      ierr=nfmpi_def_dim(ncid,"k0+k1",dim_nn(4),kkid)
      ierr=nfmpi_def_dim(ncid,"i01_global",dim_nn(5),i01id)
      ierr=nfmpi_def_dim(ncid,"i0_global",dim_nn(6),i0id)
      ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn(7),j0id)
      ierr=nfmpi_def_dim(ncid,"i01+1_global",dim_nn(8),i11id)
      ierr=nfmpi_def_dim(ncid,"j0+1_global",dim_nn(9),j11id)
      !define variables and assign varID
      ns = siid
      ierr=nfmpi_def_var(ncid,"tlz",nf_double,1,ns,itlz)
      ns = kkid
      ierr=nfmpi_def_var(ncid,"z",nf_double,1,ns,iz)
      ns = k1id
      ierr=nfmpi_def_var(ncid,"odz",nf_double,1,ns,iodz)
      ns = k0id
      ierr=nfmpi_def_var(ncid,"odzw",nf_double,1,ns,iodzw)
      nd(1) = i0id
      nd(2) = j0id
      ierr=nfmpi_def_var(ncid,"kb",nf_short,2,nd,ikb)
      ierr=nfmpi_def_var(ncid,"rmets",nf_double,2,nd,irmt)
      nd(1) = i0id
      nd(2) = j11id
      ierr=nfmpi_def_var(ncid,"iv0",nf_short,2,nd,iiv0)
      nd(1) = i11id
      nd(2) = j0id
      ierr=nfmpi_def_var(ncid,"iu0",nf_short,2,nd,iiu0)
      nt(1) = i0id
      nt(2) = j0id
      nt(3) = k0id
      ierr=nfmpi_def_var(ncid,"iw",nf_short,3,nt,iiw)! (i0,j0,k0)
      nt(3) = k1id
      ierr=nfmpi_def_var(ncid,"in",nf_short,3,nt,iin)! (i0,j0,k1)
      nt(2) = j11id
      ierr=nfmpi_def_var(ncid,"iv",nf_short,3,nt,iiv)! (i0,0:j0,k1)
      nt(1) = i11id
      nt(2) = j0id
      ierr=nfmpi_def_var(ncid,"iu",nf_short,3,nt,iiu)
      ierr=nfmpi_enddef(ncid)

      !put variables
      ns_i = 1
      nc_i = 1
      ierr=nfmpi_put_vara_double_all(ncid,itlz,ns_i,nc_i,tlz)
      nc_i = k0+k1
      ierr=nfmpi_put_vara_double_all(ncid,iz,ns_i,nc_i,z)
      nc_i = k1
      ierr=nfmpi_put_vara_double_all(ncid,iodz,ns_i,nc_i,odz)
      nc_i = k0
      ierr=nfmpi_put_vara_double_all(ncid,iodzw,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_put_vara_int2_all(ncid,ikb,ns_ij,nc_ij,kb)
      ierr=nfmpi_put_vara_double_all(ncid,irmt,ns_ij,nc_ij,rmets)
      ns_ij(2) = mylat*(j0+1)+1
      nc_ij(2) = j0+1
      ierr=nfmpi_put_vara_int2_all(ncid,iiv0,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_put_vara_int2_all(ncid,iiu0,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_put_vara_int2_all(ncid,iiw,ns_ijk,nc_ijk,iw)
      nc_ijk(3) = k1
      ierr=nfmpi_put_vara_int2_all(ncid,iin,ns_ijk,nc_ijk,in)
      ns_ijk(2) = mylat*(j0+1)+1
      nc_ijk(2) = j0+1
      ierr=nfmpi_put_vara_int2_all(ncid,iiv,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_put_vara_int2_all(ncid,iiu,ns_ijk,nc_ijk,iu)


      ierr=nfmpi_close(ncid)
end subroutine wr8nc_zkb

subroutine wr8nc_boundaries
use const,  only: i2,j2,ngx0,ngy0,k1,k2
use comm,  only: m_cart,mylat,mylon
implicit none

!declare pnetcdf var
character :: fn*18
integer :: ierr,ncid,nd(2),nt(3),nf(4)
integer :: i2id,j2id,spid,moid,k1id,k2id
integer :: isxy,itxy,itsu,issu,isnp,issp,itnp,itsp
integer(kind=mpi_offset_kind) :: dim_nn(6)
integer(kind=mpi_offset_kind) :: ns_ijkt(4),nc_ijkt(4),ns_ijk(3),nc_ijk(3)
integer(kind=mpi_offset_kind) :: ns_ij(2),nc_ij(2)

      !output file name
      fn = 'data/boundaries.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) = 10*ngy0
      dim_nn(4) = 12
      dim_nn(5) = k1
      dim_nn(6) = k2
      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,"j0-sponge_global",dim_nn(3),spid)
      ierr=nfmpi_def_dim(ncid,"month",dim_nn(4),moid)
      ierr=nfmpi_def_dim(ncid,"k1",dim_nn(5),k1id)
      ierr=nfmpi_def_dim(ncid,"k2",dim_nn(6),k2id)

      !define variables and assign varID
      nd(1) = k1id
      nd(2) = moid
      ierr=nfmpi_def_var(ncid,"sxycli",nf_double,2,nd,isxy)
      ierr=nfmpi_def_var(ncid,"txycli",nf_double,2,nd,itxy)
      nt(1) = i2id
      nt(2) = j2id
      nt(3) = moid
      ierr=nfmpi_def_var(ncid,"tsurf",nf_double,3,nt,itsu)
      ierr=nfmpi_def_var(ncid,"ssurf",nf_double,3,nt,issu)
      nf(1) = i2id
      nf(2) = spid
      nf(3) = k2id
      nf(4) = moid
      ierr=nfmpi_def_var(ncid,"snsp",nf_double,4,nf,isnp)
      ierr=nfmpi_def_var(ncid,"sssp",nf_double,4,nf,issp)
      ierr=nfmpi_def_var(ncid,"tnsp",nf_double,4,nf,itnp)
      ierr=nfmpi_def_var(ncid,"tssp",nf_double,4,nf,itsp)
      ierr=nfmpi_enddef(ncid)

      !put variables
      ns_ij(1) = 1
      ns_ij(2) = 1
      nc_ij(1) = k1
      nc_ij(2) = 12
      ierr=nfmpi_put_vara_double_all(ncid,isxy,ns_ij,nc_ij,sxycli)
      ierr=nfmpi_put_vara_double_all(ncid,itxy,ns_ij,nc_ij,txycli)
      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,issu,ns_ijk,nc_ijk,ssurf)
      ierr=nfmpi_put_vara_double_all(ncid,itsu,ns_ijk,nc_ijk,tsurf)
      ns_ijkt(1) = mylon*i2+1
      ns_ijkt(2) = mylat*10+1
      ns_ijkt(3) = 1
      ns_ijkt(4) = 1
      nc_ijkt(1) = i2
      nc_ijkt(2) = 10
      nc_ijkt(3) = k2
      nc_ijkt(4) = 12
      ierr=nfmpi_put_vara_double_all(ncid,isnp,ns_ijkt,nc_ijkt,snsp)
      ierr=nfmpi_put_vara_double_all(ncid,issp,ns_ijkt,nc_ijkt,sssp)
      ierr=nfmpi_put_vara_double_all(ncid,itnp,ns_ijkt,nc_ijkt,tnsp)
      ierr=nfmpi_put_vara_double_all(ncid,itnp,ns_ijkt,nc_ijkt,tssp)
      ierr=nfmpi_close(ncid)
end subroutine wr8nc_boundaries

subroutine wr8nc_preini
use const, only: i0,j0,k1,ngx0,ngy0
use comm, only: m_cart,mylat,mylon
implicit none
!declare pnetcdf var
character :: fn*14
integer :: ierr,ncid,nt(3)
integer :: i0id,j0id,k1id
integer :: isd,itd
integer(kind=mpi_offset_kind) :: dim_nn(3)
integer(kind=mpi_offset_kind) :: ns_ijk(3),nc_ijk(3)

      !output file name
      fn = 'data/preini.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
      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)

      !define variables and assign varID
      nt(1) = i0id
      nt(2) = j0id
      nt(3) = k1id
      ierr=nfmpi_def_var(ncid,"s1",nf_double,3,nt,isd)
      ierr=nfmpi_def_var(ncid,"t1",nf_double,3,nt,itd)
      ierr=nfmpi_enddef(ncid)

      !put variables
      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) = k1
      ierr=nfmpi_put_vara_double_all(ncid,isd,ns_ijk,nc_ijk,s1)
      ierr=nfmpi_put_vara_double_all(ncid,itd,ns_ijk,nc_ijk,t1)

      ierr=nfmpi_close(ncid)
end subroutine wr8nc_preini

subroutine rdnc_levitus(nm)
use comm, only: m_cart,myid,mylat,mylon
use const, only: i0,j0,k1,ngy0,ngx0
implicit none
character :: fn*20
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_ijkt(4),nc_ijkt(4)
integer,intent(in) :: nm
      fn = '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,s1)
      ierr=nfmpi_inq_varid(ncid,'t1_data',itmp)
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ijkt,nc_ijkt,t1)
      ierr=nfmpi_close(ncid)
endsubroutine rdnc_levitus

subroutine rdnc_depth(rmets)
use comm, only: m_cart,myid,mylat,mylon
use const, only: i0,j0
implicit none
character :: fn*14
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_ij(2),nc_ij(2)
real*8,intent(inout) :: rmets(i0,j0)
      fn = 'data/depth.nc'
      ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)
      ierr=nfmpi_inq_varid(ncid,'rmets',itmp)
      ns_ij(1) = mylon*i0+1
      ns_ij(2) = mylat*j0+1
      nc_ij(1) = i0
      nc_ij(2) = j0
      ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,rmets)
      ierr=nfmpi_close(ncid)
endsubroutine rdnc_depth

subroutine rdnc_yzgrid2(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_yzgrid2

end module