#include "options.inc"



#ifdef netcdf_diagnostics

      subroutine store_info_cdf(ncid)
c
c-----------------------------------------------------------------------
c     store some documentation about this run in a netcdf file
c-----------------------------------------------------------------------
c
      use spflame_module
      implicit none
      integer, intent(in) :: ncid
#include "netcdf.inc"
      integer :: iret,len
      character (len=200) :: text,text2

c      iret= nf_put_att_real(ncid,nf_global,'',nf_float,1,sp)

      write(text,'("SPFLAME Version ",f5.3)') version
#ifdef no_mpp
      text=text(1:len_trim(text))//' without MPP'
#else
      text=text(1:len_trim(text))//' with MPI'
#endif
      len=len_trim(text)
      iret= nf_put_att_text(ncid,nf_global,'history',len,text)
#ifndef no_mpp
c
c     domain decomposition
c
      iret= nf_put_att_int(ncid,nf_global,'n_pes_i',nf_int,1,n_pes_i)
      iret= nf_put_att_int(ncid,nf_global,'n_pes_j',nf_int,1,n_pes_j)
      iret= nf_put_att_int(ncid,nf_global,'i_blk',nf_int,1,i_blk)
      iret= nf_put_att_int(ncid,nf_global,'j_blk',nf_int,1,j_blk)
c
c     sub domain setup
c
      iret= nf_put_att_int(ncid,nf_global,'mother',nf_int,1,mother)
      iret= nf_put_att_int(ncid,nf_global,'nr_childs',nf_int,
     &                     1,nr_childs)
      if (nr_childs>0) then
       iret= nf_put_att_int(ncid,nf_global,'childs',
     &                      nf_int,nr_childs,childs(1:nr_childs))
      endif
      if (mother>0) then
       iret= nf_put_att_int(ncid,nf_global,'zoom_region',
     &                      nf_int,4,is_zoom,ie_zoom,js_zoom,je_zoom)
       iret= nf_put_att_int(ncid,nf_global,'zoom_fac',nf_int,
     &                      1,zoom_fac)
       iret= nf_put_att_int(ncid,nf_global,'zoom_fac_k',nf_int,
     &                      1,zoom_fac_k)
       iret= nf_put_att_int(ncid,nf_global,'topo_smooth_iterations',
     &                      nf_int,1,topo_smooth_iterations)
       text='no'
# ifdef sub_use_ref_topo_in_sub_domain
       text='yes'
# endif
       len=len_trim(text)
       iret= nf_put_att_text(ncid,nf_global,
     &                       'ref_topo_in_sub_domain',len,text)
       text='no'
# ifdef sub_hand_tuning_topo_in_sub_domain
       text='yes'
# endif
       iret= nf_put_att_text(ncid,nf_global,
     &              'hand_tuning_topo_in_sub_domain',len,text)
       text2 = 'normal'
#ifdef sub_wide_obc_friction_zone
       text2='wide'
#endif
#ifdef sub_narrow_obc_friction_zone
       text2='narrow'
#endif
       text='extent '//text2;text2='normal'
#ifdef sub_high_obc_friction_zone
       text2='high'
#endif
#ifdef sub_quadratic_obc_friction_zone
       text2='quadratic'
#endif
       text=text(1:len_trim(text))//' amplitude '//text2
       len=len_trim(text)
       iret= nf_put_att_text(ncid,nf_global,'obc_friction',len,text)
      endif ! mother > 0
#endif
c
c     general setup
c
#if defined linear_density
      text='linear'
#elif defined flame_density
      text='3. order polynomial taken from FLAME'
#elif defined polynomial_fit_density
      text='3. order polynomial fit'
#elif defined mcdougall_density
      text='Jackett and McDougall'
#elif defined  MWJ_density
      text='Jackett, Wright and McDougall'
#else
       error
#endif
      len=len_trim(text)
      iret= nf_put_att_text(ncid,nf_global,'eq_of_state',len,text)
      text='spherical earth'
      if (enable_beta_plane)   text='beta plane'
      if (enable_f_plane)      text='f plane'
      if (enable_rotated_grid) text='spherical earth, but rotated grid'
      len=len_trim(text)
      iret= nf_put_att_text(ncid,nf_global,'coriolis_setup',len,text)
      text = 'no'
      if (enable_blue) text = 'yes'
      len=len_trim(text)
      iret= nf_put_att_text(ncid,nf_global,'with_blue',len,text)
      text='Central differences'
      if (enable_quicker_advection) text='Quicker'
      if (enable_flux_delimiter)
     &  text=text(1:len_trim(text))//
     &    ' and flux delimiter for passive tracer'
      len=len_trim(text)
      iret= nf_put_att_text(ncid,nf_global,'advection_scheme',len,text)
      text = 'no'
#ifdef partial_cell
      text = 'yes'
#endif
      len=len_trim(text)
      iret= nf_put_att_text(ncid,nf_global,'with_partial_cell',len,text)
c
c     parameterisations
c
      text=''
      if (enable_friction_harmonic) then
        write(text2,'(" harmonic, am=",e12.6," cm^2/s")') am
        text=text(1:len_trim(text))//text2
      endif
      if (enable_friction_biharmonic) then
        write(text2,'(" biharmonic, ambi=",e12.6," cm^4/s")') ambi
        text=text(1:len_trim(text))//text2
      endif
      if (enable_friction_cosine_scal) 
     & text=text(1:len_trim(text))//' scaled with cosine of latitude'
      len=len_trim(text)
      iret= nf_put_att_text(ncid,nf_global,'hor_friction',len,text)

      text=''
      if (enable_diffusion_harmonic) then
        write(text2,'(" harmonic, ah=",e12.6," cm^2/s")') ah
        text=text(1:len_trim(text))//text2
      endif
      if (enable_diffusion_biharmonic) then
        write(text2,'(" biharmonic, ahbi=",e12.6," cm^4/s")') ahbi
        text=text(1:len_trim(text))//text2
      endif
      if (enable_diffusion_cosine_scal)  then
       if (enable_diffusion_harmonic .or.
     &    enable_diffusion_biharmonic) then
        text=text(1:len_trim(text))//' scaled with cosine of latitude'
       endif
      endif
      if (enable_diffusion_isoneutral) then
        write(text2,'(" isoneutral, ahisop=",e12.6," cm^2/s")') ahisop
        text=text(1:len_trim(text))//text2
        write(text2,'(" ,athkdf=",e12.6," cm^2/s")') athkdf
        text=text(1:len_trim(text))//text2
        write(text2,'(" ,ahsteep=",e12.6," cm^2/s")') ahsteep
        text=text(1:len_trim(text))//text2
        write(text2,'(" ,del_dm=",e12.6," s_dem=",e)') del_dm,s_dm
        text=text(1:len_trim(text))//text2
      endif
      if (enable_diffusion_isopycnic) then
        write(text2,'(" isopycnal, ahisop=",e12.6," cm^2/s")') ahisop
        text=text(1:len_trim(text))//text2
        write(text2,'(" ,athkdf=",e12.6," cm^2/s")') athkdf
        text=text(1:len_trim(text))//text2
        write(text2,'(" ,del_dm=",e12.6," s_dem=",e)') del_dm,s_dm
        text=text(1:len_trim(text))//text2
      endif
      len=len_trim(text)
      iret= nf_put_att_text(ncid,nf_global,'hor_diffusion',len,text)

      if (enable_const_vert_mixing) then
        write(text,'("constant, kappa_m=",e12.6," kappa_h=",e12.6)') 
     &            kappa_m,kappa_h
      endif
      if (enable_cgh_vert_mixing) then
        text=' CGH vertical mixing'
      endif
      if (enable_tkemix) then
        text=' TKE vertical mixing'
      endif
      if (enable_implicit_vert_fric) 
     &   text=text(1:len_trim(text))//' implicit momentum mixing'
      if (enable_implicit_vert_diff) 
     &   text=text(1:len_trim(text))//' implicit tracer mixing'
      len=len_trim(text)
      iret= nf_put_att_text(ncid,nf_global,'vert_mixing',len,text)

      end subroutine store_info_cdf



      subroutine dvcdf(ncid,ivarid,name,iname,unit,iunit,spval)
c-----------------------------------------------------------------------
c     define some standard attributes of variable ivarid in NetCDF file ncid 
c-----------------------------------------------------------------------
      implicit none
      integer ncid,ivarid,iname,iunit,iret
      character (len=*) :: name, unit
      real :: spval
#ifdef netcdf_real4
      real (kind=4) :: vv
#else
      real :: vv
#endif
#include "netcdf.inc"
      vv=spval
      call ncaptc(ncid,ivarid, 'long_name', NCCHAR,iname , name, iret) 
        if (iret.ne.0) print*,nf_strerror(iret)
      call ncaptc(ncid,ivarid, 'units',     NCCHAR,iunit, unit, iret) 
        if (iret.ne.0) print*,nf_strerror(iret)
      call ncapt (ncid,ivarid, 'missing_value',NCFLOAT,1,vv,iret)
        if (iret.ne.0) print*,nf_strerror(iret)
      call ncapt (ncid,ivarid, '_FillValue', NCFLOAT, 1,vv, iret)
        if (iret.ne.0) print*,nf_strerror(iret)
      end subroutine dvcdf
#else
      subroutine store_info_cdf(ncid)
      integer :: ncid
      end subroutine store_info_cdf
#endif
