subroutine initialize(u,v,w,s1,s2,ambient_density,u_cplx,v_cplx,w_cplx,s1_cplx,s2_cplx,
     *                      nx,ny,nz,num_dims,efactor,
     *                      Lx,Ly,Lz,dx,dy,dz,U0,DGRAD,s1_scale,s2_scale,rho_0,f,g,kappa_1,kappa_2,nu,
     *                      bc_flag,ifax,ifay,ifaz,trigx,trigy,trigz,wn,
     *                      wnx,wny,wnz,amp,p,T_diss,diss_flag,N,NM1,M,MM1,MM2,step_flag,
     *                      Rot,Ri,Re,Pr_1,Pr_2,istart,iend,i_all,i_selective,i_stat,
     *                      dt,t_start,t_end,t_write_selective,
     *                      t_stat,xt,work,scratch,myid,numprocs,locnx,locnz,
     *                      comm,twoslice,subslice,temp,netcdf_file,nyplanes,
     *                      force_cplx,force_flag,work_P,scalars)
     
c
c     Routine to set up the initial conditions, scale
c     all variables (i.e. nondimensionalize variables) and initialize
c     various arrays, flags and counters.
c

      implicit none
      integer nyplanes,nx,ny,nz,num_dims,p,iswitch,itrig,i_stat,N,NM1,M,MM1,MM2
      integer myid,numprocs,ierr,locnx,locnz,i,j,k,iwrite
      integer comm,twoslice,subslice
      integer iseed,npts
      parameter (npts=4*2048)  ! >=nx
      double precision urv(npts),vrv(npts)
      real u(nx+2,ny+1,locnz+1),v(nx+2,ny+1,locnz+1),w(nx+2,ny+1,locnz+1)
      real s1(nx+2,ny+1,locnz+1),s2(nx+2,ny+1,locnz+1)
      real ambient_density(nx+2,ny+1,locnz+1)
      real temp(nx,nyplanes,locnz+1)
      complex u_cplx(locnx,ny+1,nz+1,2),v_cplx(locnx,ny+1,nz+1,2)
      complex w_cplx(locnx,ny+1,nz+1,2),s1_cplx(locnx,ny+1,nz+1,2)
      complex s2_cplx(locnx,ny+1,nz+1,2)
      real Lx,Ly,Lz,U0,DGRAD,s1_scale,s2_scale,rho_0,f,g,kappa_1,kappa_2,nu
      real Rot,Ri,Re,Pr_1,Pr_2,dt,t_start,t_end,t_write_selective
      character*80 bc_flag,step_flag,strat_flag,read_netcdf,netcdf_file
      character*80 diss_flag,force_flag
      character*2 scalars
      character tswitch
      integer ifax(13),ifay(13),ifaz(13)
      integer jstart(4),icount(4),data_dims(4),ncid,RCODE,var_id
      real wnx(locnx),wny(ny+1),wnz(nz+1)
      real amp(locnx,ny+1,nz+1,2),T_diss
      real trigx(3*(nx+2)/2+1),trigy(2*(ny+1)),trigz(2*nz)
      real work((nx+2)*(ny+1)*(locnz+1)),efactor
      real work_P(nx+2,ny+1,locnz+1)
      complex scratch(locnx,ny+1,nz+1)
      complex force_cplx(locnx,ny+1,nz+1)
      complex wn(nz/4+1),xt(nx+2, nz/2+1)

      real pi,dx,dy,dz,x,y,z,time_scale,t_stat,total_density,N2,delta_u
      real xpos,ypos,zpos,s1_bar,s2_bar,junk,rho_bar,Tgrad
      double precision Temperature,Salinity,pMpa,density
      integer istart,iend,i_all,i_selective,kglobal
      integer dim_length,ii,jj

!
! netcdf version 3 fortran interface:
!

!
! external netcdf data types:
!
      integer nf_byte
      integer nf_int1
      integer nf_char
      integer nf_short
      integer nf_int2
      integer nf_int
      integer nf_float
      integer nf_real
      integer nf_double

      parameter (nf_byte = 1)
      parameter (nf_int1 = nf_byte)
      parameter (nf_char = 2)
      parameter (nf_short = 3)
      parameter (nf_int2 = nf_short)
      parameter (nf_int = 4)
      parameter (nf_float = 5)
      parameter (nf_real = nf_float)
      parameter (nf_double = 6)

!
! default fill values:
!
      integer           nf_fill_byte
      integer           nf_fill_int1
      integer           nf_fill_char
      integer           nf_fill_short
      integer           nf_fill_int2
      integer           nf_fill_int
      real              nf_fill_float
      real              nf_fill_real
      doubleprecision   nf_fill_double

      parameter (nf_fill_byte = -127)
      parameter (nf_fill_int1 = nf_fill_byte)
      parameter (nf_fill_char = 0)
      parameter (nf_fill_short = -32767)
      parameter (nf_fill_int2 = nf_fill_short)
      parameter (nf_fill_int = -2147483647)
      parameter (nf_fill_float = 9.9692099683868690e+36)
      parameter (nf_fill_real = nf_fill_float)
      parameter (nf_fill_double = 9.9692099683868690e+36)

!
! mode flags for opening and creating a netcdf dataset:
!
      integer nf_nowrite
      integer nf_write
      integer nf_clobber
      integer nf_noclobber
      integer nf_fill
      integer nf_nofill
      integer nf_lock
      integer nf_share

      parameter (nf_nowrite = 0)
      parameter (nf_write = 1)
      parameter (nf_clobber = 0)
      parameter (nf_noclobber = 4)
      parameter (nf_fill = 0)
      parameter (nf_nofill = 256)
      parameter (nf_lock = 1024)
      parameter (nf_share = 2048)

!
! size argument for defining an unlimited dimension:
!
      integer nf_unlimited
      parameter (nf_unlimited = 0)

!
! global attribute id:
!
      integer nf_global
      parameter (nf_global = 0)

!
! implementation limits:
!
      integer nf_max_dims
      integer nf_max_attrs
      integer nf_max_vars
      integer nf_max_name
      integer nf_max_var_dims

      parameter (nf_max_dims = 100)
      parameter (nf_max_attrs = 2000)
      parameter (nf_max_vars = 2000)
      parameter (nf_max_name = 128)
      parameter (nf_max_var_dims = nf_max_dims)

!
! error codes:
!
      integer nf_noerr
      integer nf_ebadid
      integer nf_eexist
      integer nf_einval
      integer nf_eperm
      integer nf_enotindefine
      integer nf_eindefine
      integer nf_einvalcoords
      integer nf_emaxdims
      integer nf_enameinuse
      integer nf_enotatt
      integer nf_emaxatts
      integer nf_ebadtype
      integer nf_ebaddim
      integer nf_eunlimpos
      integer nf_emaxvars
      integer nf_enotvar
      integer nf_eglobal
      integer nf_enotnc
      integer nf_ests
      integer nf_emaxname
      integer nf_eunlimit
      integer nf_enorecvars
      integer nf_echar
      integer nf_eedge
      integer nf_estride
      integer nf_ebadname
      integer nf_erange

      parameter (nf_noerr = 0)
      parameter (nf_ebadid = -33)
      parameter (nf_eexist = -35)
      parameter (nf_einval = -36)
      parameter (nf_eperm = -37)
      parameter (nf_enotindefine = -38)
      parameter (nf_eindefine = -39)
      parameter (nf_einvalcoords = -40)
      parameter (nf_emaxdims = -41)
      parameter (nf_enameinuse = -42)
      parameter (nf_enotatt = -43)
      parameter (nf_emaxatts = -44)
      parameter (nf_ebadtype = -45)
      parameter (nf_ebaddim = -46)
      parameter (nf_eunlimpos = -47)
      parameter (nf_emaxvars = -48)
      parameter (nf_enotvar = -49)
      parameter (nf_eglobal = -50)
      parameter (nf_enotnc = -51)
      parameter (nf_ests = -52)
      parameter (nf_emaxname = -53)
      parameter (nf_eunlimit = -54)
      parameter (nf_enorecvars = -55)
      parameter (nf_echar = -56)
      parameter (nf_eedge = -57)
      parameter (nf_estride = -58)
      parameter (nf_ebadname = -59)
      parameter (nf_erange = -60)

!
! error handling modes:
!
      integer  nf_fatal
      integer nf_verbose

      parameter (nf_fatal = 1)
      parameter (nf_verbose = 2)

!
! miscellaneous routines:
!
      character*80   nf_inq_libvers
      external       nf_inq_libvers

      character*80   nf_strerror
!                         (integer             ncerr)
      external       nf_strerror

      logical        nf_issyserr
!                         (integer             ncerr)
      external       nf_issyserr

!
! control routines:
!

      integer         nf_create
!                         (character*(*)       path,
!                          integer             cmode,
!                          integer             ncid)
      external        nf_create

      integer         nf__create
!                         (character*(*)       path,
!                          integer             cmode,
!                          integer             initialsz,
!                          integer             chunksizehint,
!                          integer             ncid)
      external        nf__create

      integer         nf_open
!                         (character*(*)       path,
!                          integer             mode,
!                          integer             ncid)
      external        nf_open

      integer         nf__open
!                         (character*(*)       path,
!                          integer             mode,
!                          integer             chunksizehint,
!                          integer             ncid)
      external        nf__open

      integer         nf_set_fill
!                         (integer             ncid,
!                          integer             fillmode,
!                          integer             old_mode)
      external        nf_set_fill

      integer         nf_redef
!                         (integer             ncid)
      external        nf_redef

      integer         nf_enddef
!                         (integer             ncid)
      external        nf_enddef

      integer         nf__enddef
!                         (integer             ncid,
!                          integer             h_minfree,
!                          integer             v_align,
!                          integer             v_minfree,
!                          integer             r_align)
      external        nf__enddef

      integer         nf_sync
!                         (integer             ncid)
      external        nf_sync

      integer         nf_abort
!                         (integer             ncid)
      external        nf_abort

      integer         nf_close
!                         (integer             ncid)
      external        nf_close

      integer         nf_delete
!                         (character*(*)       ncid)
      external        nf_delete

!
! general inquiry routines:
!

      integer         nf_inq
!                         (integer             ncid,
!                          integer             ndims,
!                          integer             nvars,
!                          integer             ngatts,
!                          integer             unlimdimid)
      external        nf_inq

      integer         nf_inq_ndims
!                         (integer             ncid,
!                          integer             ndims)
      external        nf_inq_ndims

      integer         nf_inq_nvars
!                         (integer             ncid,
!                          integer             nvars)
      external        nf_inq_nvars

      integer         nf_inq_natts
!                         (integer             ncid,
!                          integer             ngatts)
      external        nf_inq_natts

      integer         nf_inq_unlimdim
!                         (integer             ncid,
!                          integer             unlimdimid)
      external        nf_inq_unlimdim

!
! dimension routines:
!

      integer         nf_def_dim
!                         (integer             ncid,
!                          character(*)        name,
!                          integer             len,
!                          integer             dimid)
      external        nf_def_dim

      integer         nf_inq_dimid
!                         (integer             ncid,
!                          character(*)        name,
!                          integer             dimid)
      external        nf_inq_dimid

      integer         nf_inq_dim
!                         (integer             ncid,
!                          integer             dimid,
!                          character(*)        name,
!                          integer             len)
      external        nf_inq_dim

      integer         nf_inq_dimname
!                         (integer             ncid,
!                          integer             dimid,
!                          character(*)        name)
      external        nf_inq_dimname

      integer         nf_inq_dimlen
!                         (integer             ncid,
!                          integer             dimid,
!                          integer             len)
      external        nf_inq_dimlen

      integer         nf_rename_dim
!                         (integer             ncid,
!                          integer             dimid,
!                          character(*)        name)
      external        nf_rename_dim

!
! general attribute routines:
!

      integer         nf_inq_att
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             xtype,
!                          integer             len)
      external        nf_inq_att

      integer         nf_inq_attid
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             attnum)
      external        nf_inq_attid

      integer         nf_inq_atttype
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             xtype)
      external        nf_inq_atttype

      integer         nf_inq_attlen
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             len)
      external        nf_inq_attlen

      integer         nf_inq_attname
!                         (integer             ncid,
!                          integer             varid,
!                          integer             attnum,
!                          character(*)        name)
      external        nf_inq_attname

      integer         nf_copy_att
!                         (integer             ncid_in,
!                          integer             varid_in,
!                          character(*)        name,
!                          integer             ncid_out,
!                          integer             varid_out)
      external        nf_copy_att

      integer         nf_rename_att
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        curname,
!                          character(*)        newname)
      external        nf_rename_att

      integer         nf_del_att
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name)
      external        nf_del_att

!
! attribute put/get routines:
!

      integer         nf_put_att_text
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             len,
!                          character(*)        text)
      external        nf_put_att_text

      integer         nf_get_att_text
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          character(*)        text)
      external        nf_get_att_text

      integer         nf_put_att_int1
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             xtype,
!                          integer             len,
!                          nf_int1_t           i1vals(1))
      external        nf_put_att_int1

      integer         nf_get_att_int1
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          nf_int1_t           i1vals(1))
      external        nf_get_att_int1

      integer         nf_put_att_int2
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             xtype,
!                          integer             len,
!                          nf_int2_t           i2vals(1))
      external        nf_put_att_int2

      integer         nf_get_att_int2
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          nf_int2_t           i2vals(1))
      external        nf_get_att_int2

      integer         nf_put_att_int
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             xtype,
!                          integer             len,
!                          integer             ivals(1))
      external        nf_put_att_int

      integer         nf_get_att_int
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             ivals(1))
      external        nf_get_att_int

      integer         nf_put_att_real
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             xtype,
!                          integer             len,
!                          real                rvals(1))
      external        nf_put_att_real

      integer         nf_get_att_real
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          real                rvals(1))
      external        nf_get_att_real

      integer         nf_put_att_double
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             xtype,
!                          integer             len,
!                          double              dvals(1))
      external        nf_put_att_double

      integer         nf_get_att_double
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          double              dvals(1))
      external        nf_get_att_double

!
! general variable routines:
!

      integer         nf_def_var
!                         (integer             ncid,
!                          character(*)        name,
!                          integer             datatype,
!                          integer             ndims,
!                          integer             dimids(1),
!                          integer             varid)
      external        nf_def_var

      integer         nf_inq_var
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name,
!                          integer             datatype,
!                          integer             ndims,
!                          integer             dimids(1),
!                          integer             natts)
      external        nf_inq_var

      integer         nf_inq_varid
!                         (integer             ncid,
!                          character(*)        name,
!                          integer             varid)
      external        nf_inq_varid

      integer         nf_inq_varname
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name)
      external        nf_inq_varname

      integer         nf_inq_vartype
!                         (integer             ncid,
!                          integer             varid,
!                          integer             xtype)
      external        nf_inq_vartype

      integer         nf_inq_varndims
!                         (integer             ncid,
!                          integer             varid,
!                          integer             ndims)
      external        nf_inq_varndims

      integer         nf_inq_vardimid
!                         (integer             ncid,
!                          integer             varid,
!                          integer             dimids(1))
      external        nf_inq_vardimid

      integer         nf_inq_varnatts
!                         (integer             ncid,
!                          integer             varid,
!                          integer             natts)
      external        nf_inq_varnatts

      integer         nf_rename_var
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        name)
      external        nf_rename_var

      integer         nf_copy_var
!                         (integer             ncid_in,
!                          integer             varid,
!                          integer             ncid_out)
      external        nf_copy_var

!
! entire variable put/get routines:
!

      integer         nf_put_var_text
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        text)
      external        nf_put_var_text

      integer         nf_get_var_text
!                         (integer             ncid,
!                          integer             varid,
!                          character(*)        text)
      external        nf_get_var_text

      integer         nf_put_var_int1
!                         (integer             ncid,
!                          integer             varid,
!                          nf_int1_t           i1vals(1))
      external        nf_put_var_int1

      integer         nf_get_var_int1
!                         (integer             ncid,
!                          integer             varid,
!                          nf_int1_t           i1vals(1))
      external        nf_get_var_int1

      integer         nf_put_var_int2
!                         (integer             ncid,
!                          integer             varid,
!                          nf_int2_t           i2vals(1))
      external        nf_put_var_int2

      integer         nf_get_var_int2
!                         (integer             ncid,
!                          integer             varid,
!                          nf_int2_t           i2vals(1))
      external        nf_get_var_int2

      integer         nf_put_var_int
!                         (integer             ncid,
!                          integer             varid,
!                          integer             ivals(1))
      external        nf_put_var_int

      integer         nf_get_var_int
!                         (integer             ncid,
!                          integer             varid,
!                          integer             ivals(1))
      external        nf_get_var_int

      integer         nf_put_var_real
!                         (integer             ncid,
!                          integer             varid,
!                          real                rvals(1))
      external        nf_put_var_real

      integer         nf_get_var_real
!                         (integer             ncid,
!                          integer             varid,
!                          real                rvals(1))
      external        nf_get_var_real

      integer         nf_put_var_double
!                         (integer             ncid,
!                          integer             varid,
!                          doubleprecision     dvals(1))
      external        nf_put_var_double

      integer         nf_get_var_double
!                         (integer             ncid,
!                          integer             varid,
!                          doubleprecision     dvals(1))
      external        nf_get_var_double

!
! single variable put/get routines:
!

      integer         nf_put_var1_text
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          character*1         text)
      external        nf_put_var1_text

      integer         nf_get_var1_text
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          character*1         text)
      external        nf_get_var1_text

      integer         nf_put_var1_int1
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          nf_int1_t           i1val)
      external        nf_put_var1_int1

      integer         nf_get_var1_int1
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          nf_int1_t           i1val)
      external        nf_get_var1_int1

      integer         nf_put_var1_int2
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          nf_int2_t           i2val)
      external        nf_put_var1_int2

      integer         nf_get_var1_int2
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          nf_int2_t           i2val)
      external        nf_get_var1_int2

      integer         nf_put_var1_int
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          integer             ival)
      external        nf_put_var1_int

      integer         nf_get_var1_int
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          integer             ival)
      external        nf_get_var1_int

      integer         nf_put_var1_real
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          real                rval)
      external        nf_put_var1_real

      integer         nf_get_var1_real
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          real                rval)
      external        nf_get_var1_real

      integer         nf_put_var1_double
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          doubleprecision     dval)
      external        nf_put_var1_double

      integer         nf_get_var1_double
!                         (integer             ncid,
!                          integer             varid,
!                          integer             index(1),
!                          doubleprecision     dval)
      external        nf_get_var1_double

!
! variable array put/get routines:
!

      integer         nf_put_vara_text
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          character(*)        text)
      external        nf_put_vara_text

      integer         nf_get_vara_text
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          character(*)        text)
      external        nf_get_vara_text

      integer         nf_put_vara_int1
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          nf_int1_t           i1vals(1))
      external        nf_put_vara_int1

      integer         nf_get_vara_int1
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          nf_int1_t           i1vals(1))
      external        nf_get_vara_int1

      integer         nf_put_vara_int2
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          nf_int2_t           i2vals(1))
      external        nf_put_vara_int2

      integer         nf_get_vara_int2
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          nf_int2_t           i2vals(1))
      external        nf_get_vara_int2

      integer         nf_put_vara_int
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             ivals(1))
      external        nf_put_vara_int

      integer         nf_get_vara_int
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             ivals(1))
      external        nf_get_vara_int

      integer         nf_put_vara_real
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          real                rvals(1))
      external        nf_put_vara_real

      integer         nf_get_vara_real
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          real                rvals(1))
      external        nf_get_vara_real

      integer         nf_put_vara_double
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          doubleprecision     dvals(1))
      external        nf_put_vara_double

      integer         nf_get_vara_double
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          doubleprecision     dvals(1))
      external        nf_get_vara_double

!
! strided variable put/get routines:
!

      integer         nf_put_vars_text
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          character(*)        text)
      external        nf_put_vars_text

      integer         nf_get_vars_text
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          character(*)        text)
      external        nf_get_vars_text

      integer         nf_put_vars_int1
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          nf_int1_t           i1vals(1))
      external        nf_put_vars_int1

      integer         nf_get_vars_int1
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          nf_int1_t           i1vals(1))
      external        nf_get_vars_int1

      integer         nf_put_vars_int2
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          nf_int2_t           i2vals(1))
      external        nf_put_vars_int2

      integer         nf_get_vars_int2
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          nf_int2_t           i2vals(1))
      external        nf_get_vars_int2

      integer         nf_put_vars_int
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             ivals(1))
      external        nf_put_vars_int

      integer         nf_get_vars_int
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             ivals(1))
      external        nf_get_vars_int

      integer         nf_put_vars_real
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          real                rvals(1))
      external        nf_put_vars_real

      integer         nf_get_vars_real
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          real                rvals(1))
      external        nf_get_vars_real

      integer         nf_put_vars_double
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          doubleprecision     dvals(1))
      external        nf_put_vars_double

      integer         nf_get_vars_double
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          doubleprecision     dvals(1))
      external        nf_get_vars_double

!
! mapped variable put/get routines:
!

      integer         nf_put_varm_text
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          character(*)        text)
      external        nf_put_varm_text

      integer         nf_get_varm_text
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          character(*)        text)
      external        nf_get_varm_text

      integer         nf_put_varm_int1
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          nf_int1_t           i1vals(1))
      external        nf_put_varm_int1

      integer         nf_get_varm_int1
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          nf_int1_t           i1vals(1))
      external        nf_get_varm_int1

      integer         nf_put_varm_int2
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          nf_int2_t           i2vals(1))
      external        nf_put_varm_int2

      integer         nf_get_varm_int2
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          nf_int2_t           i2vals(1))
      external        nf_get_varm_int2

      integer         nf_put_varm_int
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          integer             ivals(1))
      external        nf_put_varm_int

      integer         nf_get_varm_int
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          integer             ivals(1))
      external        nf_get_varm_int

      integer         nf_put_varm_real
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          real                rvals(1))
      external        nf_put_varm_real

      integer         nf_get_varm_real
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          real                rvals(1))
      external        nf_get_varm_real

      integer         nf_put_varm_double
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          doubleprecision     dvals(1))
      external        nf_put_varm_double

      integer         nf_get_varm_double
!                         (integer             ncid,
!                          integer             varid,
!                          integer             start(1),
!                          integer             count(1),
!                          integer             stride(1),
!                          integer             imap(1),
!                          doubleprecision     dvals(1))
      external        nf_get_varm_double

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! begin netcdf 2.4 backward compatibility:
!

!      
! functions in the fortran interface
!
      integer nccre
      integer ncopn
      integer ncddef
      integer ncdid
      integer ncvdef
      integer ncvid
      integer nctlen
      integer ncsfil

      external nccre
      external ncopn
      external ncddef
      external ncdid
      external ncvdef
      external ncvid
      external nctlen
      external ncsfil


      integer ncrdwr
      integer nccreat
      integer ncexcl
      integer ncindef
      integer ncnsync
      integer nchsync
      integer ncndirty
      integer nchdirty
      integer nclink
      integer ncnowrit
      integer ncwrite
      integer ncclob
      integer ncnoclob
      integer ncglobal
      integer ncfill
      integer ncnofill
      integer maxncop
      integer maxncdim
      integer maxncatt
      integer maxncvar
      integer maxncnam
      integer maxvdims
      integer ncnoerr
      integer ncebadid
      integer ncenfile
      integer nceexist
      integer nceinval
      integer nceperm
      integer ncenotin
      integer nceindef
      integer ncecoord
      integer ncemaxds
      integer ncename
      integer ncenoatt
      integer ncemaxat
      integer ncebadty
      integer ncebadd
      integer ncests
      integer nceunlim
      integer ncemaxvs
      integer ncenotvr
      integer nceglob
      integer ncenotnc
      integer ncfoobar
      integer ncsyserr
      integer ncfatal
      integer ncverbos
      integer ncentool


!
! netcdf data types:
!
      integer ncbyte
      integer ncchar
      integer ncshort
      integer nclong
      integer ncfloat
      integer ncdouble

      parameter(ncbyte = 1)
      parameter(ncchar = 2)
      parameter(ncshort = 3)
      parameter(nclong = 4)
      parameter(ncfloat = 5)
      parameter(ncdouble = 6)

!     
!     masks for the struct nc flag field; passed in as 'mode' arg to
!     nccreate and ncopen.
!     

!     read/write, 0 => readonly 
      parameter(ncrdwr = 1)
!     in create phase, cleared by ncendef 
      parameter(nccreat = 2)
!     on create destroy existing file 
      parameter(ncexcl = 4)
!     in define mode, cleared by ncendef 
      parameter(ncindef = 8)
!     synchronise numrecs on change (x'10')
      parameter(ncnsync = 16)
!     synchronise whole header on change (x'20')
      parameter(nchsync = 32)
!     numrecs has changed (x'40')
      parameter(ncndirty = 64)  
!     header info has changed (x'80')
      parameter(nchdirty = 128)
!     prefill vars on endef and increase of record, the default behavior
      parameter(ncfill = 0)
!     do not fill vars on endef and increase of record (x'100')
      parameter(ncnofill = 256)
!     isa link (x'8000')
      parameter(nclink = 32768)

!     
!     'mode' arguments for nccreate and ncopen
!     
      parameter(ncnowrit = 0)
      parameter(ncwrite = ncrdwr)
      parameter(ncclob = nf_clobber)
      parameter(ncnoclob = nf_noclobber)

!     
!     'size' argument to ncdimdef for an unlimited dimension
!     
      integer ncunlim
      parameter(ncunlim = 0)

!     
!     attribute id to put/get a global attribute
!     
      parameter(ncglobal  = 0)

!     
!     advisory maximums:
!     
      parameter(maxncop = 32)
      parameter(maxncdim = 100)
      parameter(maxncatt = 2000)
      parameter(maxncvar = 2000)
!     not enforced 
      parameter(maxncnam = 128)
      parameter(maxvdims = maxncdim)

!     
!     global netcdf error status variable
!     initialized in error.c
!     

!     no error 
      parameter(ncnoerr = nf_noerr)
!     not a netcdf id 
      parameter(ncebadid = nf_ebadid)
!     too many netcdfs open 
      parameter(ncenfile = -31)   ! nc_syserr
!     netcdf file exists && ncnoclob
      parameter(nceexist = nf_eexist)
!     invalid argument 
      parameter(nceinval = nf_einval)
!     write to read only 
      parameter(nceperm = nf_eperm)
!     operation not allowed in data mode 
      parameter(ncenotin = nf_enotindefine )   
!     operation not allowed in define mode 
      parameter(nceindef = nf_eindefine)   
!     coordinates out of domain 
      parameter(ncecoord = nf_einvalcoords)
!     maxncdims exceeded 
      parameter(ncemaxds = nf_emaxdims)
!     string match to name in use 
      parameter(ncename = nf_enameinuse)   
!     attribute not found 
      parameter(ncenoatt = nf_enotatt)
!     maxncattrs exceeded 
      parameter(ncemaxat = nf_emaxatts)
!     not a netcdf data type 
      parameter(ncebadty = nf_ebadtype)
!     invalid dimension id 
      parameter(ncebadd = nf_ebaddim)
!     ncunlimited in the wrong index 
      parameter(nceunlim = nf_eunlimpos)
!     maxncvars exceeded 
      parameter(ncemaxvs = nf_emaxvars)
!     variable not found 
      parameter(ncenotvr = nf_enotvar)
!     action prohibited on ncglobal varid 
      parameter(nceglob = nf_eglobal)
!     not a netcdf file 
      parameter(ncenotnc = nf_enotnc)
      parameter(ncests = nf_ests)
      parameter (ncentool = nf_emaxname) 
      parameter(ncfoobar = 32)
      parameter(ncsyserr = -31)

!     
!     global options variable. used to determine behavior of error handler.
!     initialized in lerror.c
!     
      parameter(ncfatal = 1)
      parameter(ncverbos = 2)

!
!     default fill values.  these must be the same as in the c interface.
!
      integer filbyte
      integer filchar
      integer filshort
      integer fillong
      real filfloat
      doubleprecision fildoub

      parameter (filbyte = -127)
      parameter (filchar = 0)
      parameter (filshort = -32767)
      parameter (fillong = -2147483647)
      parameter (filfloat = 9.9692099683868690e+36)
      parameter (fildoub = 9.9692099683868690e+36)

      N=1   ! put startup fields on plane 1 
      NM1=2 ! the other field
      pi = 4.0*atan(1.0)
      dx = Lx/float(nx)
      dz = Lz/float(nz)
      time_scale = Lz/U0

      if( num_dims .eq. 2 ) then
       dy = 1
       jj=1   ! j index for writing out rho_bar
      elseif( num_dims .eq. 3 ) then
       dy = Ly/float(ny)
       jj=nyplanes/2   ! j index for writing out rho_bar
      endif

c     Will write out rho_bar(z) at Lx/2,Ly/2
      ii=nx/2
      if(myid.eq.0) open(1,file='output/rho_bar')

       do k=1,locnz+1
        kglobal = k + myid*locnz
        zpos=(kglobal-1.)*dz/Lz     ! use dimensionless values here
        do j=1,nyplanes
         do i=1,nx
           xpos=(i-1.)*dx/Lz
           ypos=(j-1.)*dy/Lz

          call s1_bar_func(xpos,ypos,zpos,s1_bar,
     *                     junk,junk,junk,junk,
     *                     Lz,s1_scale)
          call s2_bar_func(xpos,ypos,zpos,s2_bar,
     *                     junk,junk,junk,junk,
     *                     Lz,s2_scale)
          if( scalars .eq. 'TS' ) then
           Temperature = (s1_bar*s1_scale)       ! [deg C]
           Salinity = (s2_bar*s2_scale)   ! [psu]
           Salinity=Salinity/1000.     ! fortran routine rho uses concentration units
           pmPa=0.d0
           call rho(Salinity,Temperature,pMpa,density)
           rho_bar = (sngl(density) - rho_0)   ! [kg/m3]
           ambient_density(i, j, k) = rho_bar/(Lz*DGRAD)   ! store dimensionless values
          elseif( scalars .eq. 'TP' .or. scalars .eq. 'T0' ) then
           Temperature = (s1_bar*s1_scale)       ! [deg C]
           Salinity = 0.  ! [psu]
           pmPa=0.d0
           call rho(Salinity,Temperature,pMpa,density)
           rho_bar = (sngl(density) - rho_0)   ! [kg/m3]
           ambient_density(i, j, k) = rho_bar/(Lz*DGRAD)   ! store dimensionless values
          elseif( scalars .eq. 'SP' .or. scalars .eq. 'S0' ) then
           Temperature = 0.       ! [deg C]
           Salinity = (s2_bar*s2_scale)   ! [psu]
           Salinity=Salinity/1000.     ! fortran routine rho uses concentration units
           pmPa=0.d0
           call rho(Salinity,Temperature,pMpa,density)
           rho_bar = (sngl(density) - rho_0)   ! [kg/m3]
           ambient_density(i, j, k) = rho_bar/(Lz*DGRAD)   ! store dimensionless values
          elseif( scalars .eq. 'RP' .or. scalars .eq. 'R0' ) then
           ambient_density(i, j, k) = s1_bar/s1_scale      ! store dimensionless values
          elseif( scalars .eq. 'PP' ) then
           ambient_density(i, j, k) = 0.0                  ! store dimensionless values
          endif

          if(myid.eq.0 .and. i.eq.ii .and. j.eq.jj) write(1,*) zpos*Lz,rho_bar
         enddo
        enddo
       enddo
       if(myid.eq.0) close(10)
       call MPI_BARRIER(comm,ierr)   ! all processors loop and compute but only writes file
 

      if( netcdf_file .eq. 'none' .or. netcdf_file .eq. 'NONE' ) then
       do k=1,locnz+1
        kglobal = k + myid*locnz
        zpos=(kglobal-1.)*dz     ![m]

        do j=1,nyplanes
         do i=1,nx
      
          xpos=(i-1.)*dx
          ypos=(j-1.)*dy

          call user_defined_ics(xpos,ypos,zpos,Lx,Ly,Lz,
     *         u(i,j,k),v(i,j,k),w(i,j,k),s1(i,j,k),s2(i,j,k) )

         enddo
        enddo
       enddo
 
      else  ! we will read ICS from an existing netcdf file
c     for example, to read ICs from a netcdf file generated via GMwaves.x

c      define size of data blocks for reading from netcdf file 
       jstart(1)=1
       jstart(2)=1
       jstart(3)=(myid*locnz)+1
       jstart(4)=1
       icount(1)=nx
       data_dims(1)=nx
       icount(2)=nyplanes 
       data_dims(2)=nyplanes 
       if( myid .eq. numprocs-1 .and. bc_flag .eq. 'zslip') then  
        icount(3)=locnz+1
        data_dims(3)=numprocs*locnz+1
       else
        icount(3)=locnz
        data_dims(3)=numprocs*locnz
       endif
       icount(4)=1
       data_dims(4)=1

c      ACCESS THE FILE SEQUENTIALLY ACROSS PROCESSORS

       DO IWRITE=0,NUMPROCS-1
        IF( MYID .NE. IWRITE ) GOTO 999    ! wait for IWRITE to finish w/ file

       write(0,*) 'Processor ',myid,' accessing ICs file ',netcdf_file
c      processor IWRITE opens the specified file
       rcode=NF_OPEN(netcdf_file,NF_NOWRITE,ncid)

c      read in the dimensional values from the file
       do i=0, 4
        rcode=NF_INQ_DIMLEN(ncid, i, dim_length)
        if(dim_length .ne. data_dims(i)) then
         write (0,*) "dimensional values from the file are inconsistant
     * with global dimensions of the problem given."
         exit
        endif
       enddo

c      U
       rcode=NF_INQ_VARID(ncid,'uVar',var_id)
       rcode=NF_GET_VARA_REAL(ncid,var_id,jstart,icount,temp)
       do k=1,icount(3)
        do j=1,nyplanes
         do i=1,nx
          u(i,j,k)=temp(i,j,k)
         enddo
        enddo
       enddo

c      V
       rcode=NF_INQ_VARID(ncid,'vVar',var_id)
       rcode=NF_GET_VARA_REAL(ncid,var_id,jstart,icount,temp)
        do k=1,icount(3)
         do j=1,nyplanes
          do i=1,nx
           v(i,j,k)=temp(i,j,k)
          enddo
         enddo
        enddo

c      W
       rcode=NF_INQ_VARID(ncid,'wVar',var_id)
       rcode=NF_GET_VARA_REAL(ncid,var_id,jstart,icount,temp)
       do k=1,icount(3)
        do j=1,nyplanes
         do i=1,nx
          w(i,j,k)=temp(i,j,k)
         enddo
        enddo
       enddo


c      s1
       rcode=NF_INQ_VARID(ncid,'s1Var',var_id)
       rcode=NF_GET_VARA_REAL(ncid,var_id,jstart,icount,temp)

       do k=1,icount(3)      
        kglobal = k + myid*locnz
        zpos=(kglobal-1.)*dz/Lz     ![dless]
        do j=1,nyplanes 
         ypos=(j-1.)*dy/Lz
         do i=1,nx
          xpos=(i-1.)*dx/Lz
          call s1_bar_func(xpos,ypos,zpos,s1_bar,junk,junk,junk,junk,Lz,s1_scale)
          s1(i,j,k) = temp(i,j,k) - s1_bar*(s1_scale)
         enddo
        enddo
       enddo 

c      s2
       rcode=NF_INQ_VARID(ncid,'s2Var',var_id)
       rcode=NF_GET_VARA_REAL(ncid,var_id,jstart,icount,temp)
       do k=1,icount(3)      
        kglobal = k + myid*locnz
        zpos=(kglobal-1.)*dz/Lz     ![dless]
        do j=1,nyplanes 
         ypos=(j-1.)*dy/Lz
         do i=1,nx
          xpos=(i-1.)*dx/Lz
          call s2_bar_func(xpos,ypos,zpos,s2_bar,junk,junk,junk,junk,Lz,s2_scale)
          s2(i,j,k) = temp(i,j,k) - s2_bar*(s2_scale)
         enddo
        enddo
       enddo

       rcode=NF_CLOSE(ncid)
999    call MPI_BARRIER(comm,ierr)    ! all processes wait here before reading next block
      ENDDO

      endif  ! end if block for netcdf or user specified ICS

c     write out parameters characterizing simulation
c
      if( myid .eq. 0 ) then
       write(6,*) ' Number of dimensions:  ',num_dims
       write(6,*) ' nx, ny, nz : ',nx,ny,nz
       write(6,*) ' boundary conditions:  ',bc_flag
       write(6,*) ' '
       write(6,*) ' '
       write(6,*) ' domain lengths Lx, Ly, Lz [m]:  ',Lx,Ly,Lz
       write(6,*) ' grid spacings dx, dy, dz [m]:  ',dx,dy,dz
       write(6,*) ' bulk velocity scale U [m/s]:  ',U0
       write(6,*) ' characteristic density gradient scale [kg/m4]:  ',DGRAD
       write(6,*) ' characteristic range of scalar 1 [e.g. deg C]  ',s1_scale
       write(6,*) ' characteristic range of scalar 2 [e.g. psu]  ',s2_scale
       write(6,*) 
     &   ' characteristic density of Boussinesq fluid [kg/m3]:  ',rho_0
       write(6,*) ' coriolis parameter f [1/s]:  ',f
       write(6,*) ' gravitational acceleration g [m/s2]:  ',g
       write(6,*) ' molecular diffusivity of scalar 1 [m2/s]:  ',kappa_1
       write(6,*) ' molecular diffusivity of scalar 2 [m2/s]:  ',kappa_2
       write(6,*) ' viscosity [m2/s]:  ',nu
c      write(6,*) ' minimum BV period [s]:  ',per_min
c      write(6,*) ' CFL_x = Umax dt/dx:  ', UMAX*dt/dx
c      write(6,*) ' CFL_y = Vmax dt/dy:  ', VMAX*dt/dy
c      write(6,*) ' CFL_z = Wmax dt/dz:  ', WMAX*dt/dz

c     nondimensionalize the initial fields, the stratification parameters 
c     and time variables
       write(6,*) ' time step dt [s]    :  ',dt
       write(6,*) ' initial time  [s]   :  ',t_start
       write(6,*) ' final time    [s]   :  ',t_end
       write(6,*) ' Lz/U0        [s]    :  ',time_scale
      endif

      istart = 0
      iend = ( t_end - t_start ) / dt
      i_selective = t_write_selective/dt
      i_stat = t_stat/dt

      dt = dt/time_scale
      t_start = t_start/time_scale
      t_end = t_end/time_scale
      t_write_selective = t_write_selective/time_scale
      t_stat = t_stat/time_scale

      if( myid .eq. 0 ) then
       write(6,*) ' '
       write(6,*) ' time variables initialized via problem_parameters.h'
       write(6,*) ' non-dimensionalized in initilize.f'
      endif
      
      if( myid .eq. 0 ) then
       write(6,*) ' '
       write(6,*) ' time increment counters set'
       write(6,*) ' number of steps:  ',iend-istart+1
       write(6,*) ' number of steps between energetics calculations:  ',i_stat
       write(6,*) ' number of steps between full writes:  ',i_all
       write(6,*) ' number of steps between selective writes:  ',i_selective
      endif
      
      do i=0,numprocs-1
       if( myid .eq. i .and. myid .lt. numprocs-1 ) then

        open(1,file='output/u_profile',status='unknown',
     &       position='append')  ! f90 syntax

        do k=1,locnz
         write(1,*) 0.,0.,(k-1+myid*locnz)*dz,u(1,1,k)
        enddo
        close(1)
       elseif( myid .eq. i .and. myid .eq. numprocs-1 ) then

        open(1,file='output/u_profile',status='unknown',
     &       position='append')   ! f90 syntax

c '
       do k=1,locnz+1
         write(1,*) 0.,0.,(k-1+myid*locnz)*dz,u(1,1,k)
        enddo
        close(1)
        endif
         call MPI_BARRIER(comm,ierr)
        enddo
      
c     Nondimensionalize dependent variables:
      do k=1,locnz+1
       do j=1,nyplanes
        do i=1,nx
         u(i,j,k) = u(i,j,k)/U0
         v(i,j,k) = v(i,j,k)/U0
         w(i,j,k) = w(i,j,k)/U0
         s1(i,j,k) = s1(i,j,k)/(s1_scale)
         s2(i,j,k) = s2(i,j,k)/(s2_scale)
        enddo
       enddo
      enddo

      if( myid .eq. 0 ) then
       write(6,*) ' '
       write(6,*) 'dependent variables non-dimensionalized'
      endif
      
c     compute the dimensionless parameters
      Rot= f*Lz/U0
      N2 = (g/rho_0)*DGRAD
      Ri = N2*Lz**2/(U0**2)
      Re = U0*Lz/nu
      Pr_1 = nu/kappa_1
      Pr_2 = nu/kappa_2
      
      if( myid .eq. 0 ) then
       write(6,*) ' '
       write(6,*) ' '
       write(6,*) ' inverse Rossby number fL/U:  ',Rot
       write(6,*) ' bulk Richardson number (NL/U)^2 : ',Ri
       write(6,*) ' bulk Reynolds number UL/nu:  ',Re
       write(6,*) ' Prandtl number  nu/kappa_1:  ',Pr_1
       write(6,*) ' Prandtl number  nu/kappa_2:  ',Pr_2
       write(6,*) ' '
       write(6,*) ' '
      endif

c
c     set up dimensionless wavenumber vectors, 
c     N.B. for fill_wave_vectors, input length & time scales, 
c     including wavenumber truncation radius,
c     are dimensionless quantities, scaled by vertical domain size Lz & U0
c
      call fill_wave_vectors (wnx, wny, wnz, nx, ny, nz, Lx/Lz, Ly/Lz, amp,
     *                        dt, Re, Pr_1, Pr_2, p,T_diss/time_scale,diss_flag,num_dims,bc_flag,
     *                        myid, numprocs, locnx)
     
      if( myid .eq. 0 ) then
       write(6,*) ' '
       write(6,*) ' dimensionless wavenumber vectors and '
       write(6,*) ' integrating factors initialized'
       write(6,*) ' '
      
       write(6,*) 'TRANSFORMING DIMENSIONLESS INITIAL CONDITIONS P-->S'
      endif

c
c     take the 3d transforms of the dimensionless, physical space ICS
c
	iswitch = +1   ! Forward transforms: physical to wavenumber domain
	itrig = 1      ! compute & return trig tables with first call
                       ! to transform3d
	
C
C     set tswitch to indicate type of z transform, and transform u & v & s1
C
      if( bc_flag .eq. 'zperiodic' ) then
       tswitch='f'  ! indicating fourier in xy, fourier in z
      elseif( bc_flag .eq. 'zslip' ) then
       tswitch = 'c'   ! free slip, du/dz=dv/dz=ds1d/dz=0 at top & bottom
      endif
      
      call transform3d (u,u_cplx(1,1,1,N),nx,ny,nz,iswitch,tswitch,itrig,
     *                   trigx,trigy,trigz,wn,ifax,ifay,ifaz,
     *                   num_dims,work,scratch,xt,myid,numprocs,locnx,locnz,
     *                   comm,twoslice,subslice)
     
      if( myid .eq. 0 ) then
       write(6,*) ' u transformed to wavenumber domain: ', tswitch,iswitch,itrig
      endif
     
      itrig = 0      ! trig tables now computed and stored
      
      call transform3d (v,v_cplx(1,1,1,N),nx,ny,nz,iswitch,tswitch,itrig,
     *                   trigx,trigy,trigz,wn,ifax,ifay,ifaz,
     *                   num_dims,work,scratch,xt,myid,numprocs,locnx,locnz,
     *                   comm,twoslice,subslice)
      if( myid .eq. 0 ) then
       write(6,*) ' v transformed to wavenumber domain: ', tswitch,iswitch,itrig
      endif
     
      if( bc_flag .eq. 'zperiodic' ) then
       tswitch='f'
      elseif( bc_flag .eq. 'zslip' ) then
       tswitch = 's'  
      endif
      
      call transform3d (s1,s1_cplx(1,1,1,N),nx,ny,nz,iswitch,tswitch,itrig,
     *                   trigx,trigy,trigz,wn,ifax,ifay,ifaz,
     *                   num_dims,work,scratch,xt,myid,numprocs,locnx,locnz,
     *                   comm,twoslice,subslice)
      if( myid .eq. 0 ) then
       write(6,*) ' s1 transformed to wavenumber domain: ', tswitch,iswitch,itrig
      endif
      
      if( bc_flag .eq. 'zperiodic' ) then
       tswitch='f'
      elseif( bc_flag .eq. 'zslip' ) then
       tswitch = 's'  
      endif
      
      call transform3d (s2,s2_cplx(1,1,1,N),nx,ny,nz,iswitch,tswitch,itrig,
     *                   trigx,trigy,trigz,wn,ifax,ifay,ifaz,
     *                   num_dims,work,scratch,xt,myid,numprocs,locnx,locnz,
     *                   comm,twoslice,subslice)
      if( myid .eq. 0 ) then
       write(6,*) ' s2 transformed to wavenumber domain: ', tswitch,iswitch,itrig
      endif
      
C
C     set tswitch to indicate type of z transform, and transform u & v & s1 
C
      if( bc_flag .eq. 'zperiodic' ) then
       tswitch='f'
      elseif( bc_flag .eq. 'zslip' ) then
       tswitch = 's'   ! w=0 at top & bottom
      endif
      
       call transform3d (w,w_cplx(1,1,1,N),nx,ny,nz,iswitch,tswitch,itrig,
     *                  trigx,trigy,trigz,wn,ifax,ifay,ifaz,
     *                  num_dims,work,scratch,xt,myid,numprocs,locnx,locnz,
     *                  comm,twoslice,subslice)

      if( myid .eq. 0 ) then
       write(6,*) ' w transformed to wavenumber domain: ', tswitch,iswitch,itrig
      endif
     
      if( efactor .gt. 0 ) then
       call add_noise(u_cplx,v_cplx,efactor,nx,ny,nz,num_dims,
     *                wnx,wny,wnz,N,NM1,myid,numprocs,locnx,comm)
      endif

c
c     Initialize forcing wavenumber array, to be called later in
c
      if( force_flag .eq. 'yesFourier' ) then      
c        if( myid .eq. 0 ) then
c         write(6,*) ''
c         write(6,*) 'Initializing Fourier forcing  '
c         write(6,*) ''
c        endif
        call initialize_user_fforcing (force_cplx,wnx,wny,wnz,f,g,rho_0,DGRAD,Rot,U0,
     *                                    nx,ny,nyplanes,nz,locnx,Lx,Ly,Lz,
     *                                    locnz,myid,numprocs,num_dims,ifax,ifay,ifaz,
     *                                    wn,comm,twoslice,subslice,trigx,trigy,trigz,
     *                                    work_P,scratch,bc_flag,work)
      endif

c
c     initialize array pointers and time step method
c
      N=1
      NM1=2
      M=1
      MM1=2
      MM2=3
      step_flag='euler'   ! startup method, no previous information known

      call MPI_BARRIER(comm,ierr)
      return
      end