#include "options.inc"
    

c=======================================================================
c       Bolus velocity and isopycnal diffusion
c       linked in the code in driver.F
#ifdef enable_ekecalc
c       linked also in the code in eke_closure.F 
#endif
#ifdef enable_ekecalc_3D
c       linked also in the code in eke_closure_3D.F 
#endif
c=======================================================================

#ifdef enable_bolus_velocity

c write output to netcdf file
#define enable_bolus_diag

c only do isopycnal diffusion, bolus part is done elsewhere
#define enable_bolus_isopycnal_diffusion_only

c=======================================================================
c prevent spurious under shooting by advection of passive tracer
c due to isopycnal diffusion
c works only for positive tracer like concentrations
c#define delimit_tracer_fluxes
c=======================================================================

      module bolus_module
c-----------------------------------------------------------------------
c     module for bolus velocities
c-----------------------------------------------------------------------
      implicit none
      real, allocatable :: ub(:,:,:)
      real, allocatable :: vb(:,:,:)
      real, allocatable :: wb(:,:,:)
      real, allocatable :: psix(:,:,:)
      real, allocatable :: psiy(:,:,:)
      real, allocatable :: ahthk(:,:,:)
      real, parameter :: epsln = 1.0e-20 ! a small parameter
      real, parameter :: K_const=1e3 ! thickness diffusivity in m^2/s
      real, parameter :: K_min     = 10.0, K_max     = 5e4  ! threshold values bounding K_const
      real, parameter :: dslope=0.0008, slopec=0.005  ! parameters controlling max allowed isopycnal slopes
      end module bolus_module
 

      subroutine init_bolus
c-----------------------------------------------------------------------
c     initialize module
c-----------------------------------------------------------------------
      use cpflame_module
      use bolus_module
      implicit none
      real :: fxa

      if (my_pe==0) print*,' Initializing GM module '
      allocate( ub(imt,jmt,km) ); ub=0.
      allocate( vb(imt,jmt,km) ); vb=0.
      allocate( wb(imt,jmt,km) ); wb=0.
      allocate( psix(imt,jmt,km) ); psix=0.
      allocate( psiy(imt,jmt,km) ); psiy=0.
      allocate( ahthk(imt,jmt,km) ); ahthk=0.
      ahthk = K_const

#ifdef enable_bolus_diag
      if (my_pe==0) call init_bolus_diag
#else
      if (my_pe==0) print*,' no diagnostic output '
#endif
c
c-----------------------------------------------------------------------
c     Compute the grid factors which set the maximum slopes available
c     for the mixing schemes. 
c-----------------------------------------------------------------------
c
      fxa = dz*dx/(4*K_const*dt_in+epsln)
      if (my_pe==0) then
       print*,' steepest slope possible for linear stab = ',fxa
       print*,' critical slope specified as               ',slopec
#ifdef enable_bolus_isopycnal_diffusion_only
       print*,' does isopycnal diffusion only '
       print*,' assuming that bolus part is done elsewhere '
#endif
      print*,' K_const = ',k_const,' m^2/s'
      print*,' K_min = ',k_min,' m^2/s'
      print*,' K_max = ',k_max,' m^2/s'
      endif
      if (my_pe==0) print*,' done initializing GM module '
      end subroutine init_bolus




      subroutine bolus_velocity
c-----------------------------------------------------------------------
c     Verteiler funktion
c-----------------------------------------------------------------------
      use cpflame_module
      use bolus_module
#ifdef enable_diag_tracer
      use tracer_module
#endif
      implicit none
      integer :: n
#ifndef enable_bolus_isopycnal_diffusion_only
      call calculate_bolus
      call add_bolus(b)
#ifdef enable_diag_tracer
      do n=1,nt
       call add_bolus(tr(:,:,:,:,n))
      enddo
#endif
#endif
#ifdef enable_diag_tracer
      do n=1,nt
       call isopycnal_diffusion(n,tr(:,:,:,:,n))
      enddo
#endif
#ifdef enable_bolus_diag
      if (snapshot_time_step.or.initial_time==current_time) then
        call diag_bolus
      endif
#endif
      end subroutine bolus_velocity


      subroutine isopycnal_diffusion(n,tr)
c=======================================================================
c     Isopycnal diffusion
c=======================================================================
      use cpflame_module
      use bolus_module
      implicit none
      real :: tr(imt,jmt,km,0:2)
      integer :: i,j,k,js,je,n
      real :: diff_fn(imt,jmt,km),K_iso(imt,jmt,km)
      real :: diff_fe(imt,jmt,km),diff_ft(imt,jmt,km)
      real :: Tz(imt,jmt,km),bz(imt,jmt,km),sl,fxa
      real :: Tx(imt,jmt,km),bx(imt,jmt,km)
      real :: Ty(imt,jmt,km),by(imt,jmt,km)

      diff_fe(:,js_pe:je_pe,:)=0; diff_fn(:,js_pe:je_pe,:)=0
      diff_ft(:,js_pe:je_pe,:)=0

      js=max(2,js_pe); je = min(je_pe,jmt-1)
c-----------------------------------------------------------------------
c     derivatives of buoyancy
c-----------------------------------------------------------------------
      bz(:,js_pe:je_pe,:)=0; bx(:,js_pe:je_pe,:)=0
      by(:,js_pe:je_pe,:)=0
      do k=1,km-1
       do j=js,je
        do i=1,imt-1
         bx(i,j,k)= (b(i+1,j,k,taum1)-b(i,j,k,taum1))/dx*maskU(i,j,k) 
         by(i,j,k)= (b(i,j+1,k,taum1)-b(i,j,k,taum1))/dx*maskV(i,j,k) 
         bz(i,j,k)= (b(i,j,k+1,taum1)-b(i,j,k,taum1))/dz*maskW(i,j,k) 
        enddo
       enddo
      enddo
      call border_exchg3D(bx,1); call setcyclic3D(bx)
      call border_exchg3D(by,1); call setcyclic3D(by)
      call border_exchg3D(bz,1); call setcyclic3D(bz)
c-----------------------------------------------------------------------
c     derivatives of tracer
c-----------------------------------------------------------------------
      Tz(:,js_pe:je_pe,:)=0; Tx(:,js_pe:je_pe,:)=0
      Ty(:,js_pe:je_pe,:)=0
      do k=1,km-1
       do j=js,je
        do i=1,imt-1
         Tx(i,j,k)= (tr(i+1,j,k,taum1)-tr(i,j,k,taum1))/dx*maskU(i,j,k) 
         Ty(i,j,k)= (tr(i,j+1,k,taum1)-tr(i,j,k,taum1))/dx*maskV(i,j,k) 
         Tz(i,j,k)= (tr(i,j,k+1,taum1)-tr(i,j,k,taum1))/dz*maskW(i,j,k) 
        enddo
       enddo
      enddo
      call border_exchg3D(Tx,1); call setcyclic3D(Tx)
      call border_exchg3D(Ty,1); call setcyclic3D(Ty)
      call border_exchg3D(Tz,1); call setcyclic3D(Tz)

c-----------------------------------------------------------------------
c     diff_fe =  K (T_x - b_x/b_z T_z)
c-----------------------------------------------------------------------
      do k=2,km-1
       do j=js,je
        do i=1,imt-1
         fxa=(bz(i,j,k)+bz(i,j,k-1)+bz(i+1,j,k)+bz(i+1,j,k-1))
     &    /(maskW(i,j,k  )+maskW(i+1,j,k)
     &     +maskW(i,j,k-1)+maskW(i+1,j,k-1)+epsln)
         sl=bx(i,j,k)*sign(1.,fxa)/(abs(fxa)+epsln)
         sl=sl*0.5*(1.+tanh((-abs(sl)+slopec)/dslope))   
         fxa=(Tz(i,j,k)+Tz(i,j,k-1)+Tz(i+1,j,k)+Tz(i+1,j,k-1))
     &    /(maskW(i,j,k  )+maskW(i+1,j,k)
     &     +maskW(i,j,k-1)+maskW(i+1,j,k-1)+epsln)
         diff_fe(i,j,k)=(Tx(i,j,k)-sl*fxa)
     &                  *maskU(i,j,k)*ahthk(i,j,k)
        enddo
       enddo
      enddo
      call setcyclic3D(diff_fe)
c-----------------------------------------------------------------------
c     diff_fn =  K (T_y - b_y/b_z T_z)
c-----------------------------------------------------------------------
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         fxa=(bz(i,j,k)+bz(i,j,k-1)+bz(i,j+1,k)+bz(i,j+1,k-1))
     &    /(maskW(i,j,k  )+maskW(i,j+1,k)
     &     +maskW(i,j,k-1)+maskW(i,j+1,k-1)+epsln)
         sl=by(i,j,k)*sign(1.,fxa)/(abs(fxa)+epsln)
         sl=sl*0.5*(1.+tanh((-abs(sl)+slopec)/dslope))   
         fxa=(Tz(i,j,k)+Tz(i,j,k-1)+Tz(i,j+1,k)+Tz(i,j+1,k-1))
     &    /(maskW(i,j,k  )+maskW(i,j+1,k)
     &     +maskW(i,j,k-1)+maskW(i,j+1,k-1)+epsln)
         diff_fn(i,j,k)=(Ty(i,j,k)-sl*fxa)
     &                  *maskV(i,j,k)*ahthk(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(diff_fn,1); call setcyclic3D(diff_fn)
c-----------------------------------------------------------------------
c     diff_fb =  - K(b_x/b_z T_x +b_y/b_z T_y) 
c                  + T_z K((bx/bz)**2+(by/bz))**2
c      second term is done implicitly
c-----------------------------------------------------------------------
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         fxa=(bx(i,j,k)+bx(i-1,j,k)+bx(i,j,k+1)+bx(i-1,j,k+1))
     &    /(maskU(i,j,k  )+maskU(i-1,j,k  )
     &     +maskU(i,j,k+1)+maskU(i-1,j,k+1)+epsln)
         sl=fxa*sign(1.,bz(i,j,k))/(abs(bz(i,j,k))+epsln)
         sl=sl*0.5*(1.+tanh((-abs(sl)+slopec)/dslope))   
         fxa=(Tx(i,j,k)+Tx(i-1,j,k)+Tx(i,j,k+1)+Tx(i-1,j,k+1))
     &    /(maskU(i,j,k  )+maskU(i-1,j,k  )
     &     +maskU(i,j,k+1)+maskU(i-1,j,k+1)+epsln)
         diff_ft(i,j,k)=(-sl*fxa)*maskW(i,j,k)*ahthk(i,j,k)
         K_iso(i,j,k)=sl**2*maskW(i,j,k)*ahthk(i,j,k)

         fxa=(by(i,j,k)+by(i,j-1,k  )+by(i,j,k+1)+by(i,j-1,k+1))
     &    /(maskV(i,j,k  )+maskV(i,j-1,k  )
     &     +maskV(i,j,k+1)+maskV(i,j-1,k+1)+epsln)
         sl=fxa*sign(1.,bz(i,j,k))/(abs(bz(i,j,k))+epsln)
         sl=sl*0.5*(1.+tanh((-abs(sl)+slopec)/dslope))   
         fxa=(Ty(i,j,k)+Ty(i,j-1,k)+Ty(i,j,k+1)+Ty(i,j-1,k+1))
     &    /(maskV(i,j,k  )+maskV(i,j-1,k  )
     &     +maskV(i,j,k+1)+maskV(i,j-1,k+1)+epsln)
         diff_ft(i,j,k)=diff_ft(i,j,k)
     &                  +(-sl*fxa)*maskW(i,j,k)*ahthk(i,j,k)
         K_iso(i,j,k)=K_iso(i,j,k)+
     &                    sl**2*maskW(i,j,k)*ahthk(i,j,k)
        enddo
       enddo
      enddo

#if defined enable_diag_tracer & defined delimit_tracer_fluxes
       diff_fe(:,js_pe:je_pe,:)=-diff_fe(:,js_pe:je_pe,:)
       diff_fn(:,js_pe:je_pe,:)=-diff_fn(:,js_pe:je_pe,:)
       diff_ft(:,js_pe:je_pe,:)=-diff_ft(:,js_pe:je_pe,:)
       call border_exchg3D(diff_fn,1);
       call delimit_adv_flux_taup1(n,diff_fe,diff_fn,diff_ft)
       diff_fe(:,js_pe:je_pe,:)=-diff_fe(:,js_pe:je_pe,:)
       diff_fn(:,js_pe:je_pe,:)=-diff_fn(:,js_pe:je_pe,:)
       diff_ft(:,js_pe:je_pe,:)=-diff_ft(:,js_pe:je_pe,:)
       call border_exchg3D(diff_fn,1);
#endif
c---------------------------------------------------------------------------------
c     add explicit part 
c---------------------------------------------------------------------------------
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
          tr(i,j,k,taup1)=tr(i,j,k,taup1)+maskT(i,j,k)*c2dt*( 
     &          +(diff_ft(i,j,k)-diff_ft(i,j,k-1))/dz
     &          +(diff_fe(i,j,k)-diff_fe(i-1,j,k))/dx
     &          +(diff_fn(i,j,k)-diff_fn(i,j-1,k))/dx  
     &              )  
        enddo
       enddo
      enddo
c---------------------------------------------------------------------------------
c     add implicit part 
c---------------------------------------------------------------------------------
      call implicit_tracer_mix(tr,K_iso)
      call border_exchg3D(tr(:,:,:,taup1),2)
      call setcyclic3D(tr(:,:,:,taup1) )
      end subroutine isopycnal_diffusion


      subroutine bolus_set_ahthk(diff)
c-----------------------------------------------------------------------
c     set thickness diffusivity from other module
c-----------------------------------------------------------------------
      use cpflame_module
      use bolus_module
      implicit none
      real :: diff(imt,jmt,km)
      ahthk(:,js_pe:je_pe,:)=min(max(diff(:,js_pe:je_pe,:),K_min),K_max)
      call border_exchg3D(ahthk,1); call setcyclic3D(ahthk)
      end subroutine bolus_set_ahthk



      subroutine calculate_bolus
c-----------------------------------------------------------------------
c     calculate bolus velocities using diffusivity ahthk
c-----------------------------------------------------------------------
      use cpflame_module
      use bolus_module
      implicit none
      integer :: i,j,k,js,je
      real :: fxa,fxb,fxc

      real :: sx(imt,jmt,km),bz(imt,jmt,km),sy(imt,jmt,km)
      real :: a(imt,jmt,km)

      js=max(2,js_pe); je = min(je_pe,jmt-1)
c-----------------------------------------------------------------------
c     vertical derivative of buoyancy
c-----------------------------------------------------------------------
      bz(:,js_pe:je_pe,:)=0
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         bz(i,j,k)= (b(i,j,k+1,taum1)-b(i,j,k,taum1))/dz*maskW(i,j,k) 
        enddo
       enddo
      enddo
      call border_exchg3D(bz,1); call setcyclic3D(bz)
c-----------------------------------------------------------------------
c     zonal component of streamfunction
c-----------------------------------------------------------------------
      a(:,js_pe:je_pe,:)=0
      sx(:,js_pe:je_pe,:)=0
      psix(:,js_pe:je_pe,:)=0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         a(i,j,k)=(b(i+1,j,k,taum1)-b(i,j,k,taum1))/dx*maskU(i,j,k)
        enddo
       enddo
      enddo
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         fxa=(a(i,j,k+1)+a(i,j,k))/(maskU(i,j,k)+maskU(i,j,k+1)+epsln)
         fxb=(bz(i+1,j,k)+bz(i,j,k))/(maskW(i+1,j,k)+maskW(i,j,k)+epsln)
         fxa=fxa*sign(1.,fxb)/(abs(fxb)+epsln)
c         sx(i,j,k)=fxa*0.5*(1.+tanh((-abs(fxa)+slopec)/dslope))   
         sx(i,j,k)= max(min(fxa,slopec),-slopec)
        enddo
       enddo
      enddo
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         fxa=ahthk(i,j,k)+ahthk(i+1,j,k)+ahthk(i,j,k+1)+ahthk(i+1,j,k+1)
         psix(i,j,k)=fxa/4.*sx(i,j,k)*maskW(i,j,k)*maskU(i,j,k)
     &        *maskU(i,j,k+1)*maskW(i+1,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(psix,1); call setcyclic3D(psix)
c-----------------------------------------------------------------------
c     meridional component of streamfunction
c-----------------------------------------------------------------------
      a(:,js_pe:je_pe,:)=0
      sy(:,js_pe:je_pe,:)=0
      psiy(:,js_pe:je_pe,:)=0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         a(i,j,k)=(b(i,j+1,k,taum1)-b(i,j,k,taum1))/dx*maskV(i,j,k)
        enddo
       enddo
      enddo
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         fxa=(a(i,j,k+1)+a(i,j,k))/(maskV(i,j,k)+maskV(i,j,k+1)+epsln)
         fxb=(bz(i,j+1,k)+bz(i,j,k))/(maskW(i,j+1,k)+maskW(i,j,k)+epsln)
         fxa=fxa*sign(1.,fxb)/(abs(fxb)+epsln)
c         sy(i,j,k)=fxa*0.5*(1.+tanh((-abs(fxa)+slopec)/dslope))   
         sy(i,j,k)= max(min(fxa,slopec),-slopec)
        enddo
       enddo
      enddo
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         fxa=ahthk(i,j,k)+ahthk(i,j+1,k)+ahthk(i,j,k+1)+ahthk(i,j+1,k+1)
         psiy(i,j,k)=fxa/4.0*sy(i,j,k)*maskW(i,j,k)*maskV(i,j,k)
     &        *maskV(i,j,k+1)*maskW(i,j+1,k)
        enddo
       enddo
      enddo
      call border_exchg3D(psiy,1); call setcyclic3D(psiy)
c-----------------------------------------------------------------------
c        bolus velocities
c-----------------------------------------------------------------------
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         ub(i,j,k)=(psix(i,j,k)-psix(i,j,k-1))/dz*maskU(i,j,k)
         vb(i,j,k)=(psiy(i,j,k)-psiy(i,j,k-1))/dz*maskV(i,j,k)
         wb(i,j,k)=(-(psiy(i,j,k)-psiy(i,j-1,k))/dx
     &              -(psix(i,j,k)-psix(i-1,j,k))/dx)*maskW(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(ub,1); call setcyclic3D(ub)
      call border_exchg3D(vb,1); call setcyclic3D(vb)
      call border_exchg3D(wb,1); call setcyclic3D(wb)
      end subroutine calculate_bolus


      subroutine add_bolus(tr)
c-----------------------------------------------------------------------
c     add divergence of advective eddy fluxes to tracer tr
c-----------------------------------------------------------------------
      use cpflame_module
      use bolus_module
      implicit none
      real :: tr(imt,jmt,km,0:2)
      real :: adv_fe(imt,jmt,km), adv_ft(imt,jmt,km)
      real :: adv_fn(imt,jmt,km)
      integer :: i,j,k,js,je


      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do k=2,km-1
       do j=js,je
        do i=1,imt-1
         adv_fe(i,j,k)=
     &  0.5*(tr(i,j,k,tau)+tr(i+1,j,k,tau) )*ub(i,j,k)*maskU(i,j,k)
        enddo
       enddo
      enddo
      adv_fn(:,1,:)=0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         adv_fn(i,j,k)=
     &  0.5*(tr(i,j,k,tau)+tr(i,j+1,k,tau) )*vb(i,j,k)*maskV(i,j,k)
        enddo
       enddo
      enddo
      do k=1,km-1
       do j=js,je
        do i=2,imt-1
         adv_ft(i,j,k)=
     &  0.5*(tr(i,j,k,tau)+tr(i,j,k+1,tau) )*wb(i,j,k)*maskW(i,j,k)
        enddo
       enddo
      enddo
      call border_exchg3D(adv_fn,1)
      call setcyclic3D(adv_fe)
      call setcyclic3D(adv_fn)
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
         tr(i,j,k,taup1)=tr(i,j,k,taup1)+maskT(i,j,k)*c2dt*( 
     &                    -(adv_fe(i,j,k)-adv_fe(i-1,j,k))/dx
     &                    -(adv_fn(i,j,k)-adv_fn(i,j-1,k))/dx
     &                    -(adv_ft(i,j,k)-adv_ft(i,j,k-1))/dz )
        enddo
       enddo
      enddo
      call border_exchg3D(tr(:,:,:,taup1),2);
      call setcyclic3D(tr(:,:,:,taup1) )
      end subroutine add_bolus




      subroutine init_bolus_diag
c-----------------------------------------------------------------------
c     initialize NetCDF snapshot file
c-----------------------------------------------------------------------
      use cpflame_module
      use bolus_module
      implicit none
#include "netcdf.inc"
      integer :: ncid,iret,i,j,k,n
      integer :: lon_tdim,lon_udim,z_tdim,z_udim,itimedim
      integer :: lat_tdim,lat_udim,varid
      integer :: dims(4), corner(4), edges(4)
      character :: name*24, unit*16

      call def_grid_cdf('bolus.cdf')
      iret=nf_open('bolus.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      call ncredf(ncid, iret)
      iret=nf_inq_dimid(ncid,'xt',lon_tdim)
      iret=nf_inq_dimid(ncid,'xu',lon_udim)
      iret=nf_inq_dimid(ncid,'yt',lat_tdim)
      iret=nf_inq_dimid(ncid,'yu',lat_udim)
      iret=nf_inq_dimid(ncid,'zt',z_tdim)
      iret=nf_inq_dimid(ncid,'zu',z_udim)
      iret=nf_inq_dimid(ncid,'Time',itimedim)

c     2 dim variables on t grid
      dims = (/Lon_tdim,lat_udim, z_udim, iTimedim/)
      varid  = ncvdef (ncid,'psiy', NCFLOAT,4,dims,iret)
      name = 'Meridional streamfct'; unit = 'm^2/s'
      call dvcdf(ncid,varid,name,24,unit,16,spval)

      dims = (/Lon_udim,lat_tdim, z_udim, iTimedim/)
      varid  = ncvdef (ncid,'psix', NCFLOAT,4,dims,iret)
      name = 'Zonal streamfct'; unit = 'm^2/s'
      call dvcdf(ncid,varid,name,24,unit,16,spval)

      dims = (/Lon_tdim,lat_tdim, z_tdim, iTimedim/)
      varid  = ncvdef (ncid,'ahthk', NCFLOAT,4,dims,iret)
      name = 'thickness diffusivity'; unit = 'm^2/s'
      call dvcdf(ncid,varid,name,24,unit,16,spval)

      call ncclos (ncid, iret)
      end subroutine init_bolus_diag


      subroutine diag_bolus
c-----------------------------------------------------------------------
c     write to NetCDF snapshot file
c-----------------------------------------------------------------------
      use cpflame_module
      use bolus_module
      implicit none
#include "netcdf.inc"
      integer :: ncid,iret,n,npe, corner(4), edges(4)
      real :: a(imt,js_pe:je_pe,km)
      integer :: itdimid,ilen,rid,itimeid
      integer :: i,j,is,ie,js,je
      real :: fxa,ut,vt
      type(time_type) :: time

      js=max(2,js_pe); je = min(je_pe,jmt-1)

      do npe=0,n_pes
       if (my_pe==npe) then

        iret=nf_open('bolus.cdf',NF_WRITE,ncid)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        iret=nf_inq_dimid(ncid,'Time',itdimid)
        iret=nf_inq_dimlen(ncid, itdimid,ilen)
        iret=nf_inq_varid(ncid,'Time',itimeid)
        if (my_pe==0) then
         ilen=ilen+1
         time = current_time-initial_time
         fxa = time%days + time%seconds/86400.
         iret= nf_put_vara_double(ncid,itimeid,ilen,1,fxa)
        endif

        Corner = (/1,js_pe,1,ilen/); 
        edges  = (/imt,je_pe-js_pe+1,km,1/)
        iret=nf_inq_varid(ncid,'psiy',rid)
        a=psiy(:,js_pe:je_pe,:)
        iret= nf_put_vara_double(ncid,rid,corner,edges,a)

        iret=nf_inq_varid(ncid,'psix',rid)
        a=psix(:,js_pe:je_pe,:)
        iret= nf_put_vara_double(ncid,rid,corner,edges,a)

        iret=nf_inq_varid(ncid,'ahthk',rid)
        a=ahthk(:,js_pe:je_pe,:)
        iret= nf_put_vara_double(ncid,rid,corner,edges,a)

        call ncclos (ncid, iret)
       endif
       call barrier
      enddo
      end subroutine diag_bolus

#else
      subroutine bolus_dummy
      end subroutine bolus_dummy

#endif

    
