#include "options.inc"



      subroutine diag_press
c
c-----------------------------------------------------------------------
c     solve for the surface pressure, compute SSH and bottom pressure
c-----------------------------------------------------------------------
c
      use spflame_module
      implicit none
      real, dimension(is_pe:ie_pe,js_pe:je_pe,-1:1,-1:1) :: coef
      real, dimension(is_pe:ie_pe,js_pe:je_pe)           :: forc,hsp
      real, dimension(is_pe-1:ie_pe+1,js_pe-1:je_pe+1)   :: res
      real, dimension(is_pe-1:ie_pe+1,js_pe-1:je_pe+1,2) :: nabla_p
      real, dimension(is_pe-1:ie_pe+1,js_pe-1:je_pe+1)   :: sfp

      integer :: i,j,k,n,is,ie,js,je
      real :: spval=-9.9e12,tt
#ifdef netcdf_diagnostics
#include "netcdf.inc"
      integer :: npe,ncid,iret,sfpid,hspid,timedim,timeid
      integer :: corner(4),edges(4)
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:,:)
#else
      real         , allocatable :: v4(:,:)
#endif
#else
      integer :: io
      real :: buf(imt,jmt)
#endif

      if (my_pe==0) then
       print*,' diagnosing surface pressure '
      endif

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      nabla_p = 0.
      do n=1,2
       do j=js,je
        do k=1,km
         do i=is,ie
          nabla_p(i,j,n) = nabla_p(i,j,n) +
     &       (u(i,k,j,n,taup1)-u(i,k,j,n,taum1))/c2dt 
#ifdef partial_cell
     &                       *dhu(i,k,j)
#else
     &                       *dzt(k)
#endif
         enddo
        enddo
       enddo
       do j=js,je
        do i=is,ie
          nabla_p(i,j,n) = nabla_p(i,j,n)*hr(i,j)-zu(i,j,n)
        enddo
       enddo
      enddo
      call border_exchg(nabla_p(:,:,1),1,1)
      call set_cyclic(  nabla_p(:,:,1),1,1)
      call border_exchg(nabla_p(:,:,2),1,1)
      call set_cyclic(  nabla_p(:,:,2),1,1)

      coef=0.; forc=0.; res=0.; sfp=0.
      call spforc (nabla_p, forc )
      call spc9pt (coef)
      call congr (sfp, res, forc, coef, eps_surf_press,.false.)
c     ssh=sfp/(rho0*grav)/100.
      sfp=sfp*10./100000. ! cgs -> SI ,pascal -> bar
      if (my_pe==0) print*,' iterations of the solver : ',itts_solver
c
c     remove mean and checkerboard
c
      call checkerboard(sfp)
      call border_exchg(sfp,1,1)
      call set_cyclic(  sfp,1,1)

c      call zero_level(sfp,' surface pressure')
      call border_exchg(sfp,1,1)
      call set_cyclic(  sfp,1,1)
c
c     hydrostatic pressure
c
      hsp=0.
      do j=js,je
       do k=1,km
        do i=is,ie
         hsp(i,j) = hsp(i,j) - grav*rho(i,k,j)*tmask(i,k,j)
#ifdef partial_cell
     &                       *dhu(i,k,j)
#else
     &                       *dzt(k)
#endif
        enddo
       enddo
      enddo
      hsp = hsp*10 /100000. ! cgs -> SI, pascal -> bar

      if (my_pe==0) then
       print*,' --> writing diagnostic pressure to file ',
     &        diag_press_file(1:len_trim(diag_press_file))
      endif

#ifdef netcdf_diagnostics
c     write to netcdf format here
      do npe=0,n_pes
       call barrier
       if (my_pe==npe) then
        iret=nf_open(diag_press_file,NF_WRITE,ncid)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        iret=nf_inq_varid(ncid,'sfp',sfpid)
        iret=nf_inq_varid(ncid,'hsp',hspid)
        iret=nf_inq_varid(ncid,'Time',timeid)
        iret=nf_inq_dimid(ncid,'Time',timedim)
        iret=nf_inq_dimlen(ncid, timedim,n)
        allocate( v4(is_pe:ie_pe,js_pe:je_pe) )
        if (my_pe==0) then
         n=n+1
         corner(1)=n
         edges(1)=1
         call read_stamp(current_stamp,tt)
         print*,' at stamp=',current_stamp,
     &          ' (days since origin : ',tt,')',
     &          ' (time steps in file : ',n,')'
         v4(is_pe,js_pe)=tt
         iret= nf_put_vara_real (ncid,timeid, corner,edges,v4)
        endif
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        v4=sfp(is_pe:ie_pe,js_pe:je_pe)
        where( tmask(is_pe:ie_pe,1,js_pe:je_pe)==0) v4=spval
        iret= nf_put_vara_real (ncid,sfpid ,corner, edges,v4)
        v4=hsp(is_pe:ie_pe,js_pe:je_pe)
        where( tmask(is_pe:ie_pe,1,js_pe:je_pe)==0) v4=spval
        iret= nf_put_vara_real (ncid,hspid ,corner, edges,v4)
        call ncclos (ncid, iret)
       endif
       call barrier
      enddo
#else
      if (my_pe==0) then
       call getunit(io,diag_press_file,'usa ieee')
       write (io) current_stamp, itt
       write (io) imt,jmt
       write (io) spval
       write (io) xt,yt
      endif

      buf(is_pe:ie_pe,js_pe:je_pe)=sfp(is_pe:ie_pe,js_pe:je_pe)
      where(tmask(is_pe:ie_pe,1,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf

      buf(is_pe:ie_pe,js_pe:je_pe)=hsp(is_pe:ie_pe,js_pe:je_pe)
      where(tmask(is_pe:ie_pe,1,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf

      if (my_pe==0) close(io)

#endif
      end subroutine diag_press




      subroutine init_diag_press
      use spflame_module
      implicit none
#ifdef netcdf_diagnostics

#include "netcdf.inc"
      integer ncid,iret,i,j,k
      real :: spval=-9.9e12
      integer lat_tdim,itimedim, lon_tdim
      integer lat_tid,itimeid,lon_tid,sfpid,hspid
      integer dims(4), corner(4), edges(4)
      character name*24, unit*16, text*80
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:)
#else
      real         , allocatable :: v4(:)
#endif

      if (my_pe == 0) then
      print*,' initializing NetCDF output file ',
     &        diag_press_file(1:len_trim(diag_press_file))
      ncid = nccre (diag_press_file, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      call store_info_cdf(ncid)
c     dimensions
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
      iTimedim  = ncddef(ncid, 'Time', nf_unlimited, iret)
c     grid variables
      dims(1)  = Lon_tdim
      Lon_tid  = ncvdef (ncid,'Longitude_t', NCFLOAT,1,dims,iret)
      dims(1)  = Lat_tdim
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)

      dims=(/lon_tdim,lat_tdim,itimedim,1/)
      sfpid=ncvdef (ncid,'sfp', NCFLOAT,3,dims,iret)
      dims=(/lon_tdim,lat_tdim,itimedim,1/)
      hspid=ncvdef (ncid,'hsp', NCFLOAT,3,dims,iret)
c     attributes of the grid
      name = 'Latitude on T grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Longitude on T grid      '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Time                    '; unit = 'days            '
      call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
      call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
c     attributes of variables
      name = 'Sea surface pressure'; unit = 'bar'
      call dvcdf(ncid,sfpid,name,24,unit,16,spval)
      name = 'Hydrostatic pressure'; unit = 'bar'
      call dvcdf(ncid,hspid,name,24,unit,16,spval)
      call ncendf(ncid, iret)

      allocate( v4(max(imt,jmt)) )
      corner(1) = 1; edges(1) = jmt
      v4(1:jmt)=yt
      call ncvpt(ncid, Lat_tid, corner, edges,v4, iret)
      corner(1) = 1; edges(1) = imt
      v4(1:imt)=xt
      call ncvpt(ncid, lon_tid, corner, edges,v4, iret)
      deallocate(v4)
      call ncclos (ncid, iret)
      print*,' done'
      endif ! my_pe==0
#else
c     needs no intialisaztion
#endif
      end subroutine init_diag_press


