C***********************************************************************
C
C     NetCDF subroutines for pom2k.f
C
C     For version number and date, see start of subroutine write_netcdf
C
C***********************************************************************
C
      subroutine def_var_netcdf(ncid,name,nvdims,vdims,varid,
     $                          n_long_name,long_name,n_units,units,
     $                          n_coords,coords,lcoords,
     $                          nf_float,nf_noerr)
C **********************************************************************
C *                                                                    *
C * FUNCTION    :  Defines a netCDF variable and its attributes.       *
C *                                                                    *
C *                ncid ........ the netCDF I.D.                       *
C *                name ........ the variable name                     *
C *                nvdims ...... the number of dimensions of name      *
C *                vdims ....... a vector of length at least nvdims,   *
C *                              containing the dimension I.D.s        *
C *                varid ....... the variable I.D. returned            *
C *                n_long_name . the number of characters in long_name *
C *                long_name ... the long name for the variable        *
C *                n_units ..... the number of characters in units     *
C *                units ....... the units of the variable             *
C *                n_coords .... the number of characters in coords    *
C *                              (if applicable)                       *
C *                coords ...... the names of variables for the        *
C *                              "coordinates" attribute               *
C *                              (if applicable)                       *
C *                lcoords ..... .true. if a "coordinates" attribute   *
C *                              required, otherwise .false.           *
c *                nf_float .... an integer defining the netCDF        *
C *                              variable type                         *
C *                nf_noerr .... an integer defining the netCDF "no    *
C *                              error" status                         *
C *                                                                    *
C *                (nf_float and nf_noerr are declared and defined in  *
C *                 the file netcdf.inc "included" in the subroutine   *
C *                 write_netcdf.)                                     *
C *                                                                    *
C **********************************************************************
C
      integer vdims(4)
C
      integer ncid,nf_noerr,nvdims,varid,n_long_name,n_units,n_coords
C
      logical lcoords
C
      character*(*) name,long_name,units,coords
C
      integer status
C
      status=nf_def_var(ncid,name,nf_float,nvdims,vdims,varid)
C
      call handle_netcdf_error('nf_def_var          ',
     $                         status,nf_noerr)
C
C     write(6,1) name,varid
C   1 format('Variable ID returned by nf_def_var for variable ',
C    $       a7,'  = ',i5)
C
      status=nf_put_att_text(ncid,varid,'long_name',
     $                       n_long_name,long_name)
      call handle_netcdf_error('nf_put_att_text     ',
     $                         status,nf_noerr)
C
      status=nf_put_att_text(ncid,varid,'units',
     $                       n_units,units)
      call handle_netcdf_error('nf_put_att_text     ',
     $                         status,nf_noerr)
C
C     Add coordinates attribute, if necessary:
C
      if(lcoords) then
C
        status=nf_put_att_text(ncid,varid,'coordinates',
     $                         n_coords,coords)
        call handle_netcdf_error('nf_put_att_text     ',
     $                           status,nf_noerr)
C
      endif
C
      return
C
      end
C
      subroutine handle_netcdf_error(routine,status,nf_noerr)
C **********************************************************************
C *                                                                    *
C * FUNCTION    :  Checks for netCDF error.                            *
C *                                                                    *
C **********************************************************************
C
      implicit none
C
      integer status,nf_noerr
C
      character*20 routine
C
      character*80 nf_strerror
C
      if (status.ne.nf_noerr) then
C
        write(6,1) routine,nf_strerror(status)
    1   format('NetCDF routine ',a20,' terminated with error:'/2x,a80)
        stop
C
      else
C
        return
C
      endif
C
      end
C
      subroutine write_netcdf(netcdf_file,option)
C **********************************************************************
C *                                                                    *
C * FUNCTION    :  Initialises, writes to and closed netCDF file.      *
C *                                                                    *
C *                netcdf_file ....... name of output file             *
C *                option ............ 1 to initialise file            *
C *                                    2 to write a set of data        *
C *                                    3 to close file                 *
C *                                                                    *
C *                In order to include another variable:               *
C *                                                                    *
C *                1. Within the loop "if(option.eq.1) then":          *
C *                                                                    *
C *                   Insert appropriate "call def_var_netcdf(.....)", *
C *                   followed by "status=nf_put_att_<type>(.....)"    *
C *                   statements for any additional attributes. Each   *
C *                   of these latter statements should be followed by *
C *                   a "call handle_netcdf_error(.....)" statement.   *
C *                                                                    *
C *                2. Within the loop "if(option.eq.2) then":          *
C *                                                                    *
C *                   Insert appropriate                               *
C *                   "status=nf_put_var_<type>(.....)" statement,     *
C *                   followed by a "call handle_netcdf_error(.....)"  *
C *                   statement.                                       *
C *                                                                    *
C **********************************************************************
C
      implicit none
C
      include 'pom2k.c'
C
C     Filename of netCDF "include" file:
C
      include '/usr/local/include/netcdf.inc'
C
      integer option
C
      character*120 netcdf_file
C
      integer count(4),start(4),vdims(4)
C
      integer dum_varid,dvm_varid
      integer dx_varid,dy_varid
      integer east_c_varid,east_e_varid,east_u_varid,east_v_varid
      integer elb_varid
      integer fsm_varid
      integer h_varid
      integer iout
      integer ncid
      integer north_c_varid,north_e_varid,north_u_varid,north_v_varid
      integer rho_varid,rmean_varid,rot_varid
      integer status,s_varid
      integer time_dimid,time_varid,t_varid
      integer uab_varid,u_varid
      integer vab_varid,v_varid
      integer w_varid
      integer x_dimid
      integer y_dimid
      integer z_dimid,z_varid,zz_varid
C
      save dum_varid,dvm_varid
      save dx_varid,dy_varid
      save east_c_varid,east_e_varid,east_u_varid,east_v_varid
      save elb_varid
      save fsm_varid
      save h_varid
      save iout
      save ncid
      save north_c_varid,north_e_varid,north_u_varid,north_v_varid
      save rho_varid,rmean_varid,rot_varid
      save s_varid
      save time_dimid,time_varid,t_varid
      save uab_varid,u_varid
      save vab_varid,v_varid
      save w_varid
      save x_dimid
      save y_dimid
      save z_dimid,z_varid,zz_varid
C
C***********************************************************************
C
C     source_n should agree with source defined in pom2k.f. 
C
      character*40 source_n
      parameter(source_n='pom2k  2006-05-03')
C
c     if(source.ne.source_n) then
C
c       write(6,4)
c   4   format(/'Incompatible versions of program and include files ',
c    $          '..... program terminated'/)
c       stop
C
      endif
C
C***********************************************************************
C
      if(option.eq.1) then
C
C     Initialise netCDF ouput:
C
C     Initialise time index:
C
      iout=0
C
C     Create netcdf file:
C
        status=nf_create(netcdf_file,nf_clobber,ncid)
        call handle_netcdf_error('nf_create           ',
     $                           status,nf_noerr)
C       write(6,1) ncid
C   1   format('ncid returned by nf_create = ',i5)
C
C     Define global attributes:
C
        status=nf_put_att_text(ncid,nf_global,'source',
     $                         40,source)
        call handle_netcdf_error('nf_put_att_text     ',
     $                           status,nf_noerr)
C
        status=nf_put_att_text(ncid,nf_global,'title',
     $                         40,title)
        call handle_netcdf_error('nf_put_att_text     ',
     $                           status,nf_noerr)
C
C     Define dimensions:
C
        status=nf_def_dim(ncid,'time',nf_unlimited,time_dimid)
        call handle_netcdf_error('nf_def_dim          ',
     $                           status,nf_noerr)
C       write(6,2) 'time',time_dimid
C   2   format('Dimension ID returned by nf_def_dim for variable ',
C    $         a4,' = ',i5)
C
        status=nf_def_dim(ncid,'z',kb,z_dimid)
        call handle_netcdf_error('nf_def_dim          ',
     $                           status,nf_noerr)
C       write(6,2) 'z',z_dimid
C
        status=nf_def_dim(ncid,'y',jm,y_dimid)
        call handle_netcdf_error('nf_def_dim          ',
     $                           status,nf_noerr)
C       write(6,2) 'y',y_dimid
C
        status=nf_def_dim(ncid,'x',im,x_dimid)
        call handle_netcdf_error('nf_def_dim          ',
     $                           status,nf_noerr)
C       write(6,2) 'x',x_dimid
C
C     Define variables and their attributes:
C
        vdims(1)=time_dimid
C
        call def_var_netcdf(ncid,'time',1,vdims,time_varid,
     $                      4,'time',37,'days since '//time_start,
     $                      1,' ',.false.,
     $                      nf_float,nf_noerr)
C
        vdims(1)=z_dimid
C
        call def_var_netcdf(ncid,'z',1,vdims,z_varid,
     $                      18,'sigma of cell face',11,'sigma_level',
     $                      1,' ',.false.,
     $                      nf_float,nf_noerr)
        status=nf_put_att_text(ncid,z_varid,'standard_name',
     $                         22,'ocean_sigma_coordinate')
        call handle_netcdf_error('nf_put_att_text     ',
     $                           status,nf_noerr)
        status=nf_put_att_text(ncid,z_varid,'formula_terms',
     $                         26,'sigma: z eta: elb depth: h')
        call handle_netcdf_error('nf_put_att_text     ',
     $                           status,nf_noerr)
C
        call def_var_netcdf(ncid,'zz',1,vdims,zz_varid,
     $                      20,'sigma of cell centre',11,'sigma_level',
     $                      1,' ',.false.,
     $                      nf_float,nf_noerr)
        status=nf_put_att_text(ncid,zz_varid,'standard_name',
     $                         22,'ocean_sigma_coordinate')
        call handle_netcdf_error('nf_put_att_text     ',
     $                           status,nf_noerr)
        status=nf_put_att_text(ncid,zz_varid,'formula_terms',
     $                         27,'sigma: zz eta: elb depth: h')
        call handle_netcdf_error('nf_put_att_text     ',
     $                           status,nf_noerr)
C
        vdims(1)=x_dimid
        vdims(2)=y_dimid
C
        call def_var_netcdf(ncid,'dx',2,vdims,dx_varid,
     $                      19,'grid increment in x',5,'metre',
     $                      14,'east_e north_e',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'dy',2,vdims,dy_varid,
     $                      19,'grid increment in y',5,'metre',
     $                      14,'east_e north_e',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'east_u',2,vdims,east_u_varid,
     $                      19,'easting of u-points',5,'metre',
     $                      14,'east_u north_u',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'east_v',2,vdims,east_v_varid,
     $                      19,'easting of v-points',5,'metre',
     $                      14,'east_v north_v',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'east_e',2,vdims,east_e_varid,
     $                      27,'easting of elevation points',5,'metre',
     $                      14,'east_e north_e',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'east_c',2,vdims,east_c_varid,
     $                      23,'easting of cell corners',5,'metre',
     $                      14,'east_c north_c',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'north_u',2,vdims,north_u_varid,
     $                      20,'northing of u-points',5,'metre',
     $                      14,'east_u north_u',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'north_v',2,vdims,north_v_varid,
     $                      20,'northing of v-points',5,'metre',
     $                      14,'east_v north_v',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'north_e',2,vdims,north_e_varid,
     $                      28,'northing of elevation points',5,'metre',
     $                      14,'east_e north_e',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'north_c',2,vdims,north_c_varid,
     $                      24,'northing of cell corners',5,'metre',
     $                      14,'east_c north_c',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'rot',2,vdims,rot_varid,
     $               34,'Rotation angle of x-axis wrt. east',6,'degree',
     $                      14,'east_e north_e',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'h',2,vdims,h_varid,
     $                      23,'undisturbed water depth',5,'metre',
     $                      14,'east_e north_e',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'fsm',2,vdims,fsm_varid,
     $                      17,'free surface mask',13,'dimensionless',
     $                      14,'east_e north_e',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'dum',2,vdims,dum_varid,
     $                      15,'u-velocity mask',13,'dimensionless',
     $                      14,'east_u north_u',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'dvm',2,vdims,dvm_varid,
     $                      15,'v-velocity mask',13,'dimensionless',
     $                      14,'east_v north_v',.true.,
     $                      nf_float,nf_noerr)
C
        vdims(1)=x_dimid
        vdims(2)=y_dimid
        vdims(3)=z_dimid
C
        call def_var_netcdf(ncid,'rmean',3,vdims,rmean_varid,
     $                25,'horizontally-averaged rho',13,'dimensionless',
     $                      17,'east_e north_e zz',.true.,
     $                      nf_float,nf_noerr)
C
        vdims(1)=x_dimid
        vdims(2)=y_dimid
        vdims(3)=time_dimid
C
        call def_var_netcdf(ncid,'uab',3,vdims,uab_varid,
     $                      16,'depth-averaged u',9,'metre/sec',
     $                      14,'east_u north_u',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'vab',3,vdims,vab_varid,
     $                      16,'depth-averaged v',9,'metre/sec',
     $                      14,'east_v north_v',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'elb',3,vdims,elb_varid,
     $                      17,'surface elevation',5,'metre',
     $                      14,'east_e north_e',.true.,
     $                      nf_float,nf_noerr)
C
        vdims(1)=x_dimid
        vdims(2)=y_dimid
        vdims(3)=z_dimid
        vdims(4)=time_dimid
C
        call def_var_netcdf(ncid,'u',4,vdims,u_varid,
     $                      10,'x-velocity',9,'metre/sec',
     $                      17,'east_u north_u zz',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'v',4,vdims,v_varid,
     $                      10,'y-velocity',9,'metre/sec',
     $                      17,'east_v north_v zz',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'w',4,vdims,w_varid,
     $                      10,'z-velocity',9,'metre/sec',
     $                      16,'east_e north_e z',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'t',4,vdims,t_varid,
     $                      21,'potential temperature',1,'K',
     $                      17,'east_e north_e zz',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'s',4,vdims,s_varid,
     $                      23,'salinity x rho / rhoref',3,'PSS',
     $                      17,'east_e north_e zz',.true.,
     $                      nf_float,nf_noerr)
C
        call def_var_netcdf(ncid,'rho',4,vdims,rho_varid,
     $                    21,'(density-1000)/rhoref',13,'dimensionless',
     $                      17,'east_e north_e zz',.true.,
     $                      nf_float,nf_noerr)
C
C     End definitions:
C
        status=nf_enddef(ncid)
        call handle_netcdf_error('nf_enddef           ',
     $                           status,nf_noerr)
C
C     Write initial data:
C
        status=nf_put_var_real(ncid,z_varid,z)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,zz_varid,zz)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,dx_varid,dx)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,dy_varid,dy)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,east_u_varid,east_u)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,east_v_varid,east_v)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,east_e_varid,east_e)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,east_c_varid,east_c)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,north_u_varid,north_u)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,north_v_varid,north_v)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,north_e_varid,north_e)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,north_c_varid,north_c)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,rot_varid,rot)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,h_varid,h)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,fsm_varid,fsm)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,dum_varid,dum)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,dvm_varid,dvm)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
        status=nf_put_var_real(ncid,rmean_varid,rmean)
        call handle_netcdf_error('nf_put_var_real     ',
     $                           status,nf_noerr)
C
C-----------------------------------------------------------------------
C
      else if(option.eq.2) then
C
C     Write a set of data:
C
        iout=iout+1
        start(1)=iout
        count(1)=1
C
        status=nf_put_vara_real(ncid,time_varid,start,count,time)
        call handle_netcdf_error('nf_put_vara_real    ',
     $                           status,nf_noerr)
C
        start(1)=1
        start(2)=1
        start(3)=iout
        count(1)=im
        count(2)=jm
        count(3)=1
C
        status=nf_put_vara_real(ncid,uab_varid,start,count,uab)
        call handle_netcdf_error('nf_put_vara_real    ',
     $                           status,nf_noerr)
C
        status=nf_put_vara_real(ncid,vab_varid,start,count,vab)
        call handle_netcdf_error('nf_put_vara_real    ',
     $                           status,nf_noerr)
C
        status=nf_put_vara_real(ncid,elb_varid,start,count,elb)
        call handle_netcdf_error('nf_put_vara_real    ',
     $                           status,nf_noerr)
        start(1)=1
        start(2)=1
        start(3)=1
        start(4)=iout
        count(1)=im
        count(2)=jm
        count(3)=kb
        count(4)=1
C
        status=nf_put_vara_real(ncid,u_varid,start,count,u)
        call handle_netcdf_error('nf_put_vara_real    ',
     $                           status,nf_noerr)
C
        status=nf_put_vara_real(ncid,v_varid,start,count,v)
        call handle_netcdf_error('nf_put_vara_real    ',
     $                           status,nf_noerr)
C
        status=nf_put_vara_real(ncid,w_varid,start,count,w)
        call handle_netcdf_error('nf_put_vara_real    ',
     $                           status,nf_noerr)
C
        status=nf_put_vara_real(ncid,t_varid,start,count,t)
        call handle_netcdf_error('nf_put_vara_real    ',
     $                           status,nf_noerr)
C
        status=nf_put_vara_real(ncid,s_varid,start,count,s)
        call handle_netcdf_error('nf_put_vara_real    ',
     $                           status,nf_noerr)
C
        status=nf_put_vara_real(ncid,rho_varid,start,count,rho)
        call handle_netcdf_error('nf_put_vara_real    ',
     $                           status,nf_noerr)
C
C-----------------------------------------------------------------------
C
      else if(option.eq.3) then
C
C     Close file:
C
        status=nf_close(ncid)
        call handle_netcdf_error('nf_close            ',
     $                           status,nf_noerr)
C
C-----------------------------------------------------------------------
C
      else
C
        write(6,3)
    3   format(/'Invalid option for subroutine write_netcdf ..... ',
     $          'program terminated'/)
        stop
C
      endif
C
C-----------------------------------------------------------------------
C
      return
C
      end
C
C     End of source code
C
C-----------------------------------------------------------------------
C