enkf_ncdio.f90

Go to the documentation of this file.

!!$/=============================================================================/
!!$Copyright (c) 2007, The University of Massachusetts Dartmouth 
!!$Produced at the School of Marine Science & Technology 
!!$Marine Ecosystem Dynamics Modeling group
!!$All rights reserved.
!!$
!!$The FVCOM Offline Lagrangian Model has been developed by the joint UMASSD-WHOI
!!$research team.   For details of authorship and attribution of credit please see
!!$the FVCOM technical manual or contact the MEDM group.
!!$
!!$ 
!!$This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu/ The
!!$full copyright notice is contained in the file COPYRIGHT located in the root
!!$directory of the FVCOM code. This original header must be maintained in all
!!$distributed versions.
!!$
!!$THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
!!$ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
!!$IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
!!$ARE  DISCLAIMED.  
!!$
!!$/-----------------------------------------------------------------------------/
!!$CVS VERSION INFORMATION
!!$$Id$
!!$$Name$
!!$$Revision$
!!$/=============================================================================/
! VERSION 2.0
! created by Martin Huret
! modified by J. Churchill



!    *** VERSION CREATED BY JHC (04/07) TO GET TIME FROM NETCDF FILE ON NCD_READ CALL ***
!    ***   AND TO OUTPUT SURFACE ELEVATION AND BOTTOM DEPTH ON NCD_WRITE CALL         ***
!      *** modified (08/07) to re-try opening the input cdf file on error  ******

SUBROUTINE ncd_read_grid(INFILE)
  !---------------------------------------------------------------------
  ! READ DIMENSIONS IN A NETCDF FILES
  !---------------------------------------------------------------------
  USE mod_ncd
  USE lims
  IMPLICIT NONE
  !----------------------------------------------------------------------------!
  CHARACTER(LEN=100), INTENT(IN) :: INFILE
  !----------------------------------------------------------------------------!
  INTEGER            :: IERR
  INTEGER            :: N_ELEMS,N_NODES,N_SIG_M1,N_SIG
  REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: TEMP
  !----------------------------------------------------------------------------!

  !--Open NetCDF DATA FILE
  ierr = nf90_open(trim(infile),nf90_nowrite,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR READING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     stop
  END IF

  !--Get Model Dimensions
  n_elems   = getdim(nc_fid,len_trim('nele'),'nele')
  n_nodes   = getdim(nc_fid,len_trim('node'),'node')
  n_sig_m1  = getdim(nc_fid,len_trim('siglay'),'siglay')
  n_sig     = getdim(nc_fid,len_trim('siglev'),'siglev')

  m=n_nodes
  n=n_elems
  kb=n_sig    
  kbm1=n_sig_m1
  kbm2=kb-2

  !--close file
  ierr = nf90_close(nc_fid)

  RETURN
END SUBROUTINE ncd_read_grid

!==============================================================================|

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!

!==============================================================================|

SUBROUTINE ncd_read_shape(INFILE)
  !---------------------------------------------------------------------
  ! READ BATHYMETRY, SIGMA LEVELS AND GRID COEFFICIENTS IN A NETCDF FILES
  !---------------------------------------------------------------------
  USE mod_ncd
  USE all_vars
  IMPLICIT NONE
  !----------------------------------------------------------------------------!
  CHARACTER(LEN=100), INTENT(IN) :: INFILE
  !----------------------------------------------------------------------------!
  INTEGER            :: IERR
  INTEGER            :: K,I
  REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: TEMP

  !--OPEN NETCDF DATA FILE
  ierr = nf90_open(trim(infile),nf90_nowrite,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR READING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     stop
  END IF

  !--Get Node Coordinates
  ALLOCATE(temp(m,1))
  CALL getsvar(nc_fid,len_trim('lon'),'lon',m,1,temp)
  vx(1:m) = temp(1:m,1)
  WHERE(vx < 0.0_sp) vx=360.0_sp+vx
  DEALLOCATE(temp)

  ALLOCATE(temp(m,1))
  CALL getsvar(nc_fid,len_trim('lat'),'lat',m,1,temp)
  vy(1:m) = temp(1:m,1)
  DEALLOCATE(temp)

  !--Get Node Numbering
  ALLOCATE(temp(n,3))
  CALL getsvar(nc_fid,len_trim('nv'),'nv',n,3,temp)
  nv(1:n,1:3) = temp(1:n,1:3)
  DEALLOCATE(temp)
  nv(:,4) = nv(:,1)

  !--Get Bathymetry
  ALLOCATE(temp(n,1))
  CALL getsvar(nc_fid,len_trim('h'),'h',m,1,temp)
  h(1:m) = temp(1:m,1)
  DEALLOCATE(temp)

  !--Get Sigma levels
  ALLOCATE(temp(m,kb))
  CALL getsvar(nc_fid,len_trim('siglev'),'siglev',m,kb,temp)
  z(1:m,1:kb) = temp(1:m,1:kb)
  DEALLOCATE(temp)

  !--Compute derivative and intra-sigma levels
  DO k=1,kbm1
     zz(:,k)=0.5_sp*(z(:,k)+z(:,k+1))
     dz(:,k)=z(:,k)-z(:,k+1)
  END DO
  zz(:,kb)=2.0_sp*zz(:,kbm1)-zz(:,kbm2)

  DO k=1,kbm2
     dzz(:,k)=zz(:,k)-zz(:,k+1)
  END DO
  dzz(:,kb-1)=0.0
  dz(:,kb)=0.0

  DO i=1,n
    z1(i,:)   = (z(nv(i,1),:)+z(nv(i,2),:)+z(nv(i,3),:))/3.0
    zz1(i,:)  = (zz(nv(i,1),:)+zz(nv(i,2),:)+zz(nv(i,3),:))/3.0 
    dz1(i,:)  = (dz(nv(i,1),:)+dz(nv(i,2),:)+dz(nv(i,3),:))/3.0
    dzz1(i,:) = (dzz(nv(i,1),:)+dzz(nv(i,2),:)+dzz(nv(i,3),:))/3.0
  END DO
  !--Get Interpolation Parameters
  ALLOCATE(temp(n,4))
  CALL getsvar(nc_fid,len_trim('a1u'),'a1u',n,4,temp)
  a1u(1:n,:) = temp(1:n,:)
  DEALLOCATE(temp)

  ALLOCATE(temp(n,4))
  CALL getsvar(nc_fid,len_trim('a2u'),'a2u',n,4,temp)
  a2u(1:n,:) = temp(1:n,:)
  DEALLOCATE(temp)

  ALLOCATE(temp(n,3))
  CALL getsvar(nc_fid,len_trim('aw0'),'aw0',n,3,temp)
  aw0(1:n,:) = temp(1:n,:)
  DEALLOCATE(temp)

  ALLOCATE(temp(n,3))
  CALL getsvar(nc_fid,len_trim('awx'),'awx',n,3,temp)
  awx(1:n,:) = temp(1:n,:)
  DEALLOCATE(temp)

  ALLOCATE(temp(n,3))
  CALL getsvar(nc_fid,len_trim('awy'),'awy',n,3,temp)
  awy(1:n,:) = temp(1:n,:)
  DEALLOCATE(temp)

  !--Close file
  ierr = nf90_close(nc_fid)

  RETURN
END SUBROUTINE ncd_read_shape

!==============================================================================|

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
! **************MODIFIED BY JCH (04/07) TO ACQUIRE TIME FROM NETCDF FILE AT "LEVEL' OF READ
SUBROUTINE ncd_find_read_time_enkf(INFILE,Time,HO)
  USE mod_ncd
  USE all_vars, ONLY : nv
  USE lims
  IMPLICIT NONE
  !----------------------------------------------------------------------------!
  REAL(DP), INTENT(IN)                     :: time
  INTEGER, INTENT(OUT)                         :: HO
  CHARACTER(LEN=*), INTENT(IN)              :: INFILE
  !----------------------------------------------------------------------------!
  INTEGER            :: IERR
  INTEGER            :: HT,I
  REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: TEMP,TEMP2
  INTEGER            :: N_TIMES
  ierr = nf90_open(trim(infile),nf90_nowrite,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR READING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     WRITE(*,*)' Second Try'
     ierr = nf90_open(trim(infile),nf90_nowrite,nc_fid)
     IF(ierr /=nf90_noerr)THEN
        WRITE(*,*)' 2nd ERROR READING ',trim(infile)
        WRITE(*,*)trim(nf90_strerror(ierr))
                WRITE(*,*)' Third Try'
        ierr = nf90_open(trim(infile),nf90_nowrite,nc_fid)
        IF(ierr /=nf90_noerr)THEN
           WRITE(*,*)' 3rd ERROR READING ',trim(infile)
           WRITE(*,*)trim(nf90_strerror(ierr))
                   WRITE(*,*)' Yer Out!'
                   stop
                END IF
          END IF
  END IF
  !---------------------------------------------------------------------
  ! Read Data from file INFILE at time level ht
  !---------------------------------------------------------------------
  n_times   = getdim(nc_fid,len_trim('time'),'time') ! number of times in file

  ALLOCATE(temp(n_times,1))

  CALL getsvar_d(nc_fid,len_trim('time'),'time',n_times,1,temp)

  do i=1, n_times
     if (abs(temp(i,1)-time)<1e-6) then
         ho=i-1 ! NEED TO MINUS BECAUSE WE ASSUME THE FIRST ONE IS ZERO
         print *, "FIND THE DATA INDEX MATCHING CURRENT TIME STEP :", ho+1 , temp(ho+1,1)
         RETURN
     end if
  end do
  DEALLOCATE(temp)

  ierr = nf90_close(nc_fid)
  print *, 'dont find the nc index match current time:', time
  stop
  RETURN

END SUBROUTINE ncd_find_read_time_enkf
SUBROUTINE ncd_read_enkf(INFILE,UL,VL,T1L,S1L,ELL,HO)
  !---------------------------------------------------------------------
  ! READ DATA FROM DAILY NETCDF FILES
  !---------------------------------------------------------------------

  USE mod_ncd
  USE all_vars, ONLY : nv
  USE lims
  IMPLICIT NONE

  !----------------------------------------------------------------------------!
  REAL(SP), DIMENSION(0:NGL,KB),INTENT(OUT)   :: UL,VL 
  REAL(SP), DIMENSION(0:MGL,KB),INTENT(OUT)   :: T1L,S1L 
  ! REAL(SP), DIMENSION(0:M,KB),INTENT(OUT) :: T1L
  REAL(SP), DIMENSION(0:MGL),INTENT(OUT)      :: ELL
!  REAL(SP), INTENT(OUT)                     :: time
  INTEGER, INTENT(IN)                         :: HO
  CHARACTER(LEN=*), INTENT(IN)              :: INFILE
  !----------------------------------------------------------------------------!
  INTEGER            :: IERR
  INTEGER            :: HT,I
  REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: TEMP,TEMP2
  INTEGER            :: N_TIMES
  ul=0.0
  vl=0.0
  t1l=0.0
  s1l=0.0
  ell=0.0

  !--Adjustement to read in Netcdf file
  ht=ho+1
  print *, 'nc read from:', infile, " in hour: ",ht  
  !--Open NetCDF Datafile 

!  **** MODIFIED 08/07 BY JHC TO RETRY OPENING THE INPUT FILE TWICE ON ERROR

  ierr = nf90_open(trim(infile),nf90_nowrite,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR READING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     WRITE(*,*)' Second Try'
     ierr = nf90_open(trim(infile),nf90_nowrite,nc_fid)
     IF(ierr /=nf90_noerr)THEN
        WRITE(*,*)' 2nd ERROR READING ',trim(infile)
        WRITE(*,*)trim(nf90_strerror(ierr))
        WRITE(*,*)' Third Try'
        ierr = nf90_open(trim(infile),nf90_nowrite,nc_fid)
        IF(ierr /=nf90_noerr)THEN
           WRITE(*,*)' 3rd ERROR READING ',trim(infile)
           WRITE(*,*)trim(nf90_strerror(ierr))
           WRITE(*,*)' Yer Out!'
           stop
        END IF
      END IF
  END IF

  !---------------------------------------------------------------------
  ! Read Data from file INFILE at time level ht
  !---------------------------------------------------------------------

  ! ----- ADDED BY JHC (04/07) TO GET time (SECONDS AFTER START TIME) OF EACH READ------

!  N_TIMES   = GETDIM(NC_FID,LEN_TRIM('time'),'time') ! number of times in file

!  ALLOCATE(TEMP(HT,1))
  
!  CALL GETSVAR(NC_FID,LEN_TRIM('time'),'time',HT,1,TEMP)
  
  
!  time = TEMP(HT,1)
!  DEALLOCATE(TEMP)

  !--free surface elevation
  print *, 'finish open nc file'
  ALLOCATE(temp(mgl,1))
  CALL getdvar(nc_fid,len_trim('zeta'),'zeta',mgl,1,temp,ht)
  ell(1:mgl) = temp(1:mgl,1)
  DEALLOCATE(temp)
  print *, 'finish reading el'
  !--salinity
  ALLOCATE(temp(mgl,kbm1))
  CALL getdvar(nc_fid,len_trim('salinity'),'salinity',mgl,kbm1,temp,ht)
  s1l(1:mgl,1:kbm1) = temp(1:mgl,1:kbm1)
  DEALLOCATE(temp)
  print *, 'finish reading salinity'
  ! temperature
   ALLOCATE(temp(mgl,kbm1))
   CALL getdvar(nc_fid,len_trim('temp'),'temp',mgl,kbm1,temp,ht)
   t1l(1:mgl,1:kbm1) = temp(1:mgl,1:kbm1)
   DEALLOCATE(temp)
  print *, 'finish reading temp'

  !--U velocity 
  ALLOCATE(temp(ngl,kbm1))
  CALL getdvar(nc_fid,len_trim('u'),'u',ngl,kbm1,temp,ht)
  ul(1:ngl,1:kbm1) = temp(1:ngl,1:kbm1)
  DEALLOCATE(temp)
  print *, 'finish reading u'
  !--V velocity
  ALLOCATE(temp(ngl,kbm1))
  CALL getdvar(nc_fid,len_trim('v'),'v',ngl,kbm1,temp,ht)
  vl(1:ngl,1:kbm1) = temp(1:ngl,1:kbm1)
  DEALLOCATE(temp)
  print *, 'finish reading v'
!----------------------------here you have to check what name and format(node or cell) you output
  !--WW velocity
!  ALLOCATE(TEMP(N,KBM1))
!  CALL GETDVAR(NC_FID,LEN_TRIM('ww'),'ww',N,KBM1,TEMP,HT)
!  WWL(1:N,1:KBM1) = TEMP(1:N,1:KBM1)
!  DEALLOCATE(TEMP)
  !--W velocity (omega)
!  ALLOCATE(TEMP(M,KBM1))
!  ALLOCATE(TEMP2(N,KBM1))
!print *, 'do not read omega, set to zero:'
!temp=0.0
!  CALL GETDVAR(NC_FID,LEN_TRIM('omega'),'omega',M,KBM1,TEMP,HT)  !original
! CALL GETDVAR(NC_FID,LEN_TRIM('wts'),'wts',M,KBM1,TEMP,HT) !pengfei some case people output as wts, maybe from cell(N) or from node(M)
!  DO I = 1,N
!    TEMP2(I,:) = ((TEMP(NV(I,1),:))+(TEMP(NV(I,2),:))+(TEMP(NV(I,3),:)))/3.0
!  END DO  
!  WWL(1:N,1:KBM1) = TEMP2(1:N,1:KBM1)
!  DEALLOCATE(TEMP,TEMP2)

  !--KH
!  ALLOCATE(TEMP(M,KBM1))
!print *, 'do not read km, set to zero'
!temp=0.0
!  CALL GETDVAR(NC_FID,LEN_TRIM('km'),'km',M,KBM1,TEMP,HT)
!  KHL(1:M,1:KBM1) = TEMP(1:M,1:KBM1)
!  DEALLOCATE(TEMP)

  !--Close file
  ierr = nf90_close(nc_fid)
  
  RETURN
END SUBROUTINE ncd_read_enkf

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!

!==============================================================================|


! **************MODIFIED BY JCH (04/07) TO ACQUIRE TIME FROM NETCDF FILE AT "LEVEL' OF READ

SUBROUTINE ncd_read(INFILE,UL,VL,WWL,KHL,ELL,time,HO)
  !---------------------------------------------------------------------
  ! READ DATA FROM DAILY NETCDF FILES
  !---------------------------------------------------------------------

  USE mod_ncd
  USE all_vars, ONLY : nv
  USE lims
  IMPLICIT NONE

  !----------------------------------------------------------------------------!
  REAL(SP), DIMENSION(0:N,KB),INTENT(OUT)   :: UL,VL,WWL 
  REAL(SP), DIMENSION(0:M,KB),INTENT(OUT)   :: KHL 
  ! REAL(SP), DIMENSION(0:M,KB),INTENT(OUT) :: T1L
  REAL(SP), DIMENSION(0:M),INTENT(OUT)      :: ELL
  REAL(SP), INTENT(OUT)                     :: time
  INTEGER, INTENT(IN)                         :: HO
  CHARACTER(LEN=100), INTENT(IN)              :: INFILE
  !----------------------------------------------------------------------------!
  INTEGER            :: IERR
  INTEGER            :: HT,I
  REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: TEMP,TEMP2
  INTEGER            :: N_TIMES


  !--Adjustement to read in Netcdf file
  ht=ho+1
  print *, 'nc read from:', infile, "in hour: ",ht  
  !--Open NetCDF Datafile 

!  **** MODIFIED 08/07 BY JHC TO RETRY OPENING THE INPUT FILE TWICE ON ERROR

  ierr = nf90_open(trim(infile),nf90_nowrite,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR READING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     WRITE(*,*)' Second Try'
     ierr = nf90_open(trim(infile),nf90_nowrite,nc_fid)
     IF(ierr /=nf90_noerr)THEN
        WRITE(*,*)' 2nd ERROR READING ',trim(infile)
        WRITE(*,*)trim(nf90_strerror(ierr))
        WRITE(*,*)' Third Try'
        ierr = nf90_open(trim(infile),nf90_nowrite,nc_fid)
        IF(ierr /=nf90_noerr)THEN
           WRITE(*,*)' 3rd ERROR READING ',trim(infile)
           WRITE(*,*)trim(nf90_strerror(ierr))
           WRITE(*,*)' Yer Out!'
           stop
        END IF
      END IF
  END IF

  !---------------------------------------------------------------------
  ! Read Data from file INFILE at time level ht
  !---------------------------------------------------------------------

  ! ----- ADDED BY JHC (04/07) TO GET time (SECONDS AFTER START TIME) OF EACH READ------

  n_times   = getdim(nc_fid,len_trim('time'),'time') ! number of times in file

  ALLOCATE(temp(ht,1))
  
  CALL getsvar(nc_fid,len_trim('time'),'time',ht,1,temp)
  
  
  time = temp(ht,1)
  DEALLOCATE(temp)

  !--free surface elevation
  ALLOCATE(temp(m,1))
  CALL getdvar(nc_fid,len_trim('zeta'),'zeta',m,1,temp,ht)
  ell(1:m) = temp(1:m,1)
  DEALLOCATE(temp)

  !--salinity
  !ALLOCATE(TEMP(M,KBM1))
  !CALL GETDVAR(NC_FID,LEN_TRIM('salinity'),'salinity',M,KBM1,TEMP,HT)
  !S1L(1:M,1:KBM1) = TEMP(1:M,1:KBM1)
  !DEALLOCATE(TEMP)

  ! temperature
  ! ALLOCATE(TEMP(M,KBM1))
  ! CALL GETDVAR(NC_FID,LEN_TRIM('temp'),'temp',M,KBM1,TEMP,HT)
  ! T1L(1:M,1:KBM1) = TEMP(1:M,1:KBM1)
  ! DEALLOCATE(TEMP)


  !--U velocity 
  ALLOCATE(temp(n,kbm1))
  CALL getdvar(nc_fid,len_trim('u'),'u',n,kbm1,temp,ht)
  ul(1:n,1:kbm1) = temp(1:n,1:kbm1)
  DEALLOCATE(temp)

  !--V velocity
  ALLOCATE(temp(n,kbm1))
  CALL getdvar(nc_fid,len_trim('v'),'v',n,kbm1,temp,ht)
  vl(1:n,1:kbm1) = temp(1:n,1:kbm1)
  DEALLOCATE(temp)
!----------------------------here you have to check what name and format(node or cell) you output
  !--WW velocity
!  ALLOCATE(TEMP(N,KBM1))
!  CALL GETDVAR(NC_FID,LEN_TRIM('ww'),'ww',N,KBM1,TEMP,HT)
!  WWL(1:N,1:KBM1) = TEMP(1:N,1:KBM1)
!  DEALLOCATE(TEMP)
  !--W velocity (omega)
  ALLOCATE(temp(m,kbm1))
  ALLOCATE(temp2(n,kbm1))
print *, 'do not read omega, set to zero:'
temp=0.0
!  CALL GETDVAR(NC_FID,LEN_TRIM('omega'),'omega',M,KBM1,TEMP,HT)  !original
! CALL GETDVAR(NC_FID,LEN_TRIM('wts'),'wts',M,KBM1,TEMP,HT) !pengfei some case people output as wts, maybe from cell(N) or from node(M)
  DO i = 1,n
    temp2(i,:) = ((temp(nv(i,1),:))+(temp(nv(i,2),:))+(temp(nv(i,3),:)))/3.0
  END DO  
  wwl(1:n,1:kbm1) = temp2(1:n,1:kbm1)
  DEALLOCATE(temp,temp2)

  !--KH
  ALLOCATE(temp(m,kbm1))
print *, 'do not read km, set to zero'
temp=0.0
!  CALL GETDVAR(NC_FID,LEN_TRIM('km'),'km',M,KBM1,TEMP,HT)
  khl(1:m,1:kbm1) = temp(1:m,1:kbm1)
  DEALLOCATE(temp)

  !--Close file
  ierr = nf90_close(nc_fid)
  
  RETURN
END SUBROUTINE ncd_read

!==============================================================================|

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!

!==============================================================================|
! MODIFIED BY JHC 04/07 TO OUTPUT SURFACE ELEVATION (EP) AND BOTTOM DEPTH (HP)
!          Modified by JHC 07/07 to include the inwater variable in the output NCD file

!SUBROUTINE NCD_WRITE(INFILE,NPTS,TIME,LABEL,INDOMAIN,XP,YP,ZP,UP,VP,WP,NT)
SUBROUTINE ncd_write_el(INFILE,nnode,nlayer,temp)
USE mod_ncd
IMPLICIT NONE
  INTEGER, INTENT(IN)                   :: nnode,nlayer
  CHARACTER(LEN=120), INTENT(IN)        :: INFILE
  CHARACTER(LEN=100)                  :: varname
  CHARACTER(LEN=100)                  :: TSTRING,NETCDF_CONVENTION
  INTEGER                             :: IERR,varid
  INTEGER                             :: TIME_DID,TIME_VID,NLAG_DID,IND_VID,LAB_VID
  INTEGER                             :: X_VID,Y_VID,Z_VID,U_VID,V_VID,W_VID
  INTEGER                             :: DYNMTIME(1)
  INTEGER                             :: STAT1D(1),DYNM1D(2)
  REAL(SP), DIMENSION(nnode,nlayer)   :: TEMP
      integer :: RecordDimID
      integer :: nRecords
      integer :: dims(2)
      character(len = nf90_max_name) :: RecordDimName

  varname='zeta'

  ierr = nf90_open(trim(infile),nf90_write,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR OPENING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     stop
  END IF
  ierr = nf90_inquire(nc_fid, unlimiteddimid = recorddimid)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire unlimited id for el: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  ierr = nf90_inquire_dimension(nc_fid, recorddimid, &
                                      name = recorddimname, len = nrecords)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire length of unlimited id for el: ' 
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  print *, "length of unlimited dimension are : ",nrecords
  dims(1)=1
  dims(2)=nrecords


  !write time to file
!  CALL PUTDVAR(NC_FID,LEN_TRIM('time'),'time',1,TEMP,NT)
   ierr = nf90_inq_varid(nc_fid,trim(varname),varid)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error getting variable id: ',trim(varname)
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF

    ierr = nf90_put_var(nc_fid,varid,temp,start=dims)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'3 error getting variable: ',trim(varname)
!       print *, temp,dims ! should be mark
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF
     !--Close file
     ierr = nf90_close(nc_fid)
END SUBROUTINE ncd_write_el 

SUBROUTINE ncd_write_u(INFILE,ncell,nlayer,temp)
USE mod_ncd
IMPLICIT NONE
  INTEGER, INTENT(IN)                   :: ncell,nlayer
  CHARACTER(LEN=120), INTENT(IN)        :: INFILE
  CHARACTER(LEN=100)                  :: varname
  CHARACTER(LEN=100)                  :: TSTRING,NETCDF_CONVENTION
  INTEGER                             :: IERR,VARID
  INTEGER                             :: TIME_DID,TIME_VID,NLAG_DID,IND_VID,LAB_VID
  INTEGER                             :: X_VID,Y_VID,Z_VID,U_VID,V_VID,W_VID
  INTEGER                             :: DYNMTIME(1)
  INTEGER                             :: STAT1D(1),DYNM1D(2)
  REAL(SP), DIMENSION(ncell,nlayer)   :: TEMP
  REAL(SP), DIMENSION(ncell,nlayer,1)   :: TEMP3
      integer :: RecordDimID
      integer :: nRecords
      integer :: dims(3)
      character(len = nf90_max_name) :: RecordDimName
! temp3(:,:,1)=temp(:,:) 
 varname='u'
  print *, 'open nc file in write u: ',trim(infile)
  ierr = nf90_open(trim(infile),nf90_write,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR OPENING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     stop
  END IF
  ierr = nf90_inquire(nc_fid, unlimiteddimid = recorddimid)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire unlimited id for u: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  ierr = nf90_inquire_dimension(nc_fid, recorddimid, &
                                      name = recorddimname, len = nrecords)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire length of unlimited id for u: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  print *, "length of unlimited dimension are : ",nrecords
  dims(1)=1
  dims(2)=1
  dims(3)=nrecords

  !write time to file
!  CALL PUTDVAR(NC_FID,LEN_TRIM('time'),'time',1,TEMP,NT)
   ierr = nf90_inq_varid(nc_fid,trim(varname),varid)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error getting variable id: ',trim(varname)
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF

    ierr = nf90_put_var(nc_fid,varid,temp,start=dims)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'3 error getting variable: ',trim(varname)
!       print *, temp,dims ! should be mark
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF
     !--Close file
     ierr = nf90_close(nc_fid)



END SUBROUTINE ncd_write_u


SUBROUTINE ncd_write_v(INFILE,ncell,nlayer,temp)
USE mod_ncd
IMPLICIT NONE
  INTEGER, INTENT(IN)                   :: ncell,nlayer
  CHARACTER(LEN=120), INTENT(IN)        :: INFILE
  CHARACTER(LEN=100)                   :: varname
  CHARACTER(LEN=100)                  :: TSTRING,NETCDF_CONVENTION
  INTEGER                             :: IERR,VARID
  INTEGER                             :: TIME_DID,TIME_VID,NLAG_DID,IND_VID,LAB_VID
  INTEGER                             :: X_VID,Y_VID,Z_VID,U_VID,V_VID,W_VID
  INTEGER                             :: DYNMTIME(1)
  INTEGER                             :: STAT1D(1),DYNM1D(2)
  REAL(SP), DIMENSION(ncell,nlayer)   :: TEMP
      integer :: RecordDimID
      integer :: nRecords
      integer :: dims(3)
      character(len = nf90_max_name) :: RecordDimName
  varname='v'

  ierr = nf90_open(trim(infile),nf90_write,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR OPENING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     stop
  END IF

  !write time to file
!  CALL PUTDVAR(NC_FID,LEN_TRIM('time'),'time',1,TEMP,NT)
   ierr = nf90_inq_varid(nc_fid,trim(varname),varid)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error getting variable id: ',trim(varname)
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF
  ierr = nf90_inquire(nc_fid, unlimiteddimid = recorddimid)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire unlimited id for v: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  ierr = nf90_inquire_dimension(nc_fid, recorddimid, &
                                      name = recorddimname, len = nrecords)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire length of unlimited id for v: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  print *, "length of unlimited dimension are : ",nrecords
  dims(1)=1
  dims(2)=1
  dims(3)=nrecords

    ierr = nf90_put_var(nc_fid,varid,temp,start=dims)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'3 error getting variable: ',trim(varname)
!       print *, temp,dims ! should be mark
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF
     !--Close file
     ierr = nf90_close(nc_fid)


END SUBROUTINE ncd_write_v



!#################################
SUBROUTINE ncd_write_u10(INFILE,ncell,nlayer,temp)
USE mod_ncd
IMPLICIT NONE
  INTEGER, INTENT(IN)                   :: ncell,nlayer
  CHARACTER(LEN=120), INTENT(IN)        :: INFILE
  CHARACTER(LEN=100)                  :: varname
  CHARACTER(LEN=100)                  :: TSTRING,NETCDF_CONVENTION
  INTEGER                             :: IERR,VARID
  INTEGER                             :: TIME_DID,TIME_VID,NLAG_DID,IND_VID,LAB_VID
  INTEGER                             :: X_VID,Y_VID,Z_VID,U_VID,V_VID,W_VID
  INTEGER                             :: DYNMTIME(1)
  INTEGER                             :: STAT1D(1),DYNM1D(2)
  REAL(SP), DIMENSION(ncell,nlayer)   :: TEMP
  REAL(SP), DIMENSION(ncell,nlayer,1)   :: TEMP3
  integer :: num
      integer :: RecordDimID
      integer :: nRecords
      integer :: dims(3),cnts(3)
      character(len = nf90_max_name) :: RecordDimName
! temp3(:,:,1)=temp(:,:) 
 num=10
 varname='u'
  print *, 'open nc file in write u: ',trim(infile)
  ierr = nf90_open(trim(infile),nf90_write,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR OPENING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     stop
  END IF
  ierr = nf90_inquire(nc_fid, unlimiteddimid = recorddimid)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire unlimited id for u: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  ierr = nf90_inquire_dimension(nc_fid, recorddimid, &
                                      name = recorddimname, len = nrecords)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire length of unlimited id for u: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  print *, "length of unlimited dimension are : ",nrecords
  dims(1)=1
  dims(2)=1
  dims(3)=nrecords
  cnts(1)=ncell
  cnts(2)=num
  cnts(3)=1

  !write time to file
!  CALL PUTDVAR(NC_FID,LEN_TRIM('time'),'time',1,TEMP,NT)
   ierr = nf90_inq_varid(nc_fid,trim(varname),varid)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error getting variable id: ',trim(varname)
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF

    ierr = nf90_put_var(nc_fid,varid,temp(:,1:num),start=dims,count=cnts)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'3 error getting variable: ',trim(varname)
!       print *, temp,dims ! should be mark
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF
     !--Close file
     ierr = nf90_close(nc_fid)



END SUBROUTINE ncd_write_u10


SUBROUTINE ncd_write_v10(INFILE,ncell,nlayer,temp)
USE mod_ncd
IMPLICIT NONE
  INTEGER, INTENT(IN)                   :: ncell,nlayer
  CHARACTER(LEN=120), INTENT(IN)        :: INFILE
  CHARACTER(LEN=100)                   :: varname
  CHARACTER(LEN=100)                  :: TSTRING,NETCDF_CONVENTION
  INTEGER                             :: IERR,VARID
  INTEGER                             :: TIME_DID,TIME_VID,NLAG_DID,IND_VID,LAB_VID
  INTEGER                             :: X_VID,Y_VID,Z_VID,U_VID,V_VID,W_VID
  INTEGER                             :: DYNMTIME(1)
  INTEGER                             :: STAT1D(1),DYNM1D(2)
  REAL(SP), DIMENSION(ncell,nlayer)   :: TEMP
  integer :: num
      integer :: RecordDimID
      integer :: nRecords
      integer :: dims(3),cnts(3)
      character(len = nf90_max_name) :: RecordDimName
  varname='v'
  num=10
  ierr = nf90_open(trim(infile),nf90_write,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR OPENING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     stop
  END IF

  !write time to file
!  CALL PUTDVAR(NC_FID,LEN_TRIM('time'),'time',1,TEMP,NT)
   ierr = nf90_inq_varid(nc_fid,trim(varname),varid)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error getting variable id: ',trim(varname)
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF
  ierr = nf90_inquire(nc_fid, unlimiteddimid = recorddimid)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire unlimited id for v: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  ierr = nf90_inquire_dimension(nc_fid, recorddimid, &
                                      name = recorddimname, len = nrecords)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire length of unlimited id for v: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  print *, "length of unlimited dimension are : ",nrecords
  dims(1)=1
  dims(2)=1
  dims(3)=nrecords  
  cnts(1)=ncell
  cnts(2)=num
  cnts(3)=1

    ierr = nf90_put_var(nc_fid,varid,temp(:,1:num),start=dims,count=cnts)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'3 error getting variable: ',trim(varname)
!       print *, temp,dims ! should be mark
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF
     !--Close file
     ierr = nf90_close(nc_fid)


END SUBROUTINE ncd_write_v10


!################################


SUBROUTINE ncd_write_t(INFILE,nnode,nlayer,temp)
USE mod_ncd
IMPLICIT NONE
  INTEGER, INTENT(IN)                   :: nnode,nlayer
  CHARACTER(LEN=120), INTENT(IN)        :: INFILE
  CHARACTER(LEN=100)                  :: varname
  CHARACTER(LEN=100)                  :: TSTRING,NETCDF_CONVENTION
  INTEGER                             :: IERR,VARID
  INTEGER                             :: TIME_DID,TIME_VID,NLAG_DID,IND_VID,LAB_VID
  INTEGER                             :: X_VID,Y_VID,Z_VID,U_VID,V_VID,W_VID
  INTEGER                             :: DYNMTIME(1)
  INTEGER                             :: STAT1D(1),DYNM1D(2)
  REAL(SP), DIMENSION(nnode,nlayer)   :: TEMP
      integer :: RecordDimID
      integer :: nRecords
      integer :: dims(3)
      character(len = nf90_max_name) :: RecordDimName

  varname='temp'

  ierr = nf90_open(trim(infile),nf90_write,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR OPENING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     stop
  END IF

  !write time to file
!  CALL PUTDVAR(NC_FID,LEN_TRIM('time'),'time',1,TEMP,NT)
   ierr = nf90_inq_varid(nc_fid,trim(varname),varid)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error getting variable id: ',trim(varname)
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF
  ierr = nf90_inquire(nc_fid, unlimiteddimid = recorddimid)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire unlimited id for t: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  ierr = nf90_inquire_dimension(nc_fid, recorddimid, &
                                      name = recorddimname, len = nrecords)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire length of unlimited id for t: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  print *, "length of unlimited dimension are : ",nrecords
  dims(1)=1
  dims(2)=1
  dims(3)=nrecords
    ierr = nf90_put_var(nc_fid,varid,temp,start=dims)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'3 error getting variable: ',trim(varname)
!       print *, temp,dims ! should be mark
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF
     !--Close file
     ierr = nf90_close(nc_fid)
END SUBROUTINE ncd_write_t


SUBROUTINE ncd_write_s(INFILE,nnode,nlayer,temp)
USE mod_ncd
IMPLICIT NONE
  INTEGER, INTENT(IN)                   ::nnode,nlayer
  CHARACTER(LEN=120), INTENT(IN)        :: INFILE
  CHARACTER(LEN=100)                  :: varname
  CHARACTER(LEN=100)                  :: TSTRING,NETCDF_CONVENTION
  INTEGER                             :: IERR,VARID
  INTEGER                             :: TIME_DID,TIME_VID,NLAG_DID,IND_VID,LAB_VID
  INTEGER                             :: X_VID,Y_VID,Z_VID,U_VID,V_VID,W_VID
  INTEGER                             :: DYNMTIME(1)
  INTEGER                             :: STAT1D(1),DYNM1D(2)
  REAL(SP), DIMENSION(nnode,nlayer)   :: TEMP
      integer :: RecordDimID
      integer :: nRecords
      integer :: dims(3)
      character(len = nf90_max_name) :: RecordDimName

  varname='salinity'

  ierr = nf90_open(trim(infile),nf90_write,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR OPENING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     stop
  END IF

  !write time to file
!  CALL PUTDVAR(NC_FID,LEN_TRIM('time'),'time',1,TEMP,NT)
   ierr = nf90_inq_varid(nc_fid,trim(varname),varid)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error getting variable id: ',trim(varname)
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF

  ierr = nf90_inquire(nc_fid, unlimiteddimid = recorddimid)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire unlimited id for s: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  ierr = nf90_inquire_dimension(nc_fid, recorddimid, &
                                      name = recorddimname, len = nrecords)
  IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'error inquire length of unlimited id for s: '
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
  END IF
  print *, "length of unlimited dimension are : ",nrecords
  dims(1)=1
  dims(2)=1
  dims(3)=nrecords



    ierr = nf90_put_var(nc_fid,varid,temp,start=dims)
    IF(ierr /=nf90_noerr)THEN
       WRITE(*,*)'3 error getting variable: ',trim(varname)
!       print *, temp,dims ! should be mark
       WRITE(*,*)trim(nf90_strerror(ierr))
       stop
    END IF
     !--Close file
     ierr = nf90_close(nc_fid)
END SUBROUTINE ncd_write_s













!=======================================================================================
SUBROUTINE ncd_write(INFILE,NPTS,TIME,LABEL,INDOMAIN,XP,YP,ZP,UP,VP,WP,EP,HP,INWATER,NT)

  USE mod_ncd
  IMPLICIT NONE
  !----------------------------------------------------------------------------! 
  INTEGER, INTENT(IN)                   :: NPTS, NT
  INTEGER, DIMENSION(NPTS),INTENT(IN)   :: LABEL,INDOMAIN,INWATER   ! 07/07 ADDED INWATER
  REAL(SP), DIMENSION(NPTS),INTENT(IN)  :: XP,YP,ZP,UP,VP,WP,EP,HP  ! 04/07 ADDED EP & HP
  REAL(SP), INTENT(IN)                  :: TIME
  CHARACTER(LEN=100), INTENT(IN)        :: INFILE
  !----------------------------------------------------------------------------! 
  CHARACTER(LEN=100)                  :: TSTRING,NETCDF_CONVENTION
  INTEGER                             :: IERR
  INTEGER                             :: TIME_DID,TIME_VID,NLAG_DID,IND_VID,LAB_VID
  INTEGER                             :: X_VID,Y_VID,Z_VID,U_VID,V_VID,W_VID
  INTEGER                             :: DYNMTIME(1)
  INTEGER                             :: STAT1D(1),DYNM1D(2)
  REAL(SP), ALLOCATABLE, DIMENSION(:) :: TEMP
  !----------------------------------------------------------------------------!

  IF(nt == 1)THEN

     !--Create File --> nc_fid (modified by JHC 4/07 to allow for input of file into MATLAB)
     ierr = nf90_create(trim(infile) ,nf90_clobber, nc_fid)
     IF(ierr /=nf90_noerr)THEN
        WRITE(*,*)'error creating',trim(infile)
        WRITE(*,*)trim(nf90_strerror(ierr))
        stop
     END IF

     !--Netcdf Convention String
     netcdf_convention = 'CF-1.0'

     !Global Attributes     
     tstring = "FVCOM Offline Lagrangian Particle Data"
     ierr = nf90_put_att(nc_fid,nf90_global,"title",tstring)
     ierr = nf90_put_att(nc_fid,nf90_global,"institution","SMAST") 
     ierr = nf90_put_att(nc_fid,nf90_global,"source"     ,"OFFLINE_FVCOM") 
     ierr = nf90_put_att(nc_fid,nf90_global,"modeler"    ,"PHIL MCCRACKEN") 
     ierr = nf90_put_att(nc_fid,nf90_global,"Conventions",trim(netcdf_convention))

     !Dimensioning
     ierr = nf90_def_dim(nc_fid,"nlag" ,npts,nlag_did)
     ierr = nf90_def_dim(nc_fid,"time" ,nf90_unlimited,time_did)

     dynmtime     = (/time_did/)                 !!Time
     dynm1d       = (/nlag_did,time_did/)        !!Dynamic 1d var 
     stat1d       = (/nlag_did/)                 !!Static 1d var 
     !--Variable Definitions

     !!====time===============================================!
     ierr = nf90_def_var(nc_fid,"time",nf90_float,dynmtime,time_vid)
     ierr = nf90_put_att(nc_fid,time_vid,"long_name","time")
     ierr = nf90_put_att(nc_fid,time_vid,"units","seconds")

     !!=====label==============================================!
     ierr = nf90_def_var(nc_fid,"label",nf90_float,stat1d,lab_vid)
     ierr = nf90_put_att(nc_fid,lab_vid,"long_name","particle label")
     ierr = nf90_put_att(nc_fid,lab_vid,"units","")

     !!=====indomain=================================================!
     ierr = nf90_def_var(nc_fid,"indomain",nf90_float,dynm1d,ind_vid)
     ierr = nf90_put_att(nc_fid,ind_vid,"long_name","particle indomain (1)")
     ierr = nf90_put_att(nc_fid,ind_vid,"units","")
     !!=====x=================================================!
     ierr = nf90_def_var(nc_fid,"x",nf90_float,dynm1d,x_vid)
     ierr = nf90_put_att(nc_fid,x_vid,"long_name","particle x position")
     ierr = nf90_put_att(nc_fid,x_vid,"units","m")

     !!=====y=================================================!
     ierr = nf90_def_var(nc_fid,"y",nf90_float,dynm1d,y_vid)
     ierr = nf90_put_att(nc_fid,y_vid,"long_name","particle y position")
     ierr = nf90_put_att(nc_fid,y_vid,"units","m")

     !!=====z=================================================!
     ierr = nf90_def_var(nc_fid,"z",nf90_float,dynm1d,z_vid)
     ierr = nf90_put_att(nc_fid,z_vid,"long_name","particle z position")
     ierr = nf90_put_att(nc_fid,z_vid,"units","m")

     !!=====u=================================================!
     ierr = nf90_def_var(nc_fid,"u",nf90_float,dynm1d,u_vid)
     ierr = nf90_put_att(nc_fid,u_vid,"long_name","particle u velocity")
     ierr = nf90_put_att(nc_fid,u_vid,"units","cm/s")

     !!=====v=================================================!
     ierr = nf90_def_var(nc_fid,"v",nf90_float,dynm1d,v_vid)
     ierr = nf90_put_att(nc_fid,v_vid,"long_name","particle v velocity")
     ierr = nf90_put_att(nc_fid,v_vid,"units","cm/s")

     !!=====w=================================================!
     ierr = nf90_def_var(nc_fid,"omega",nf90_float,dynm1d,w_vid)
     ierr = nf90_put_att(nc_fid,w_vid,"long_name","particle w velocity")
     ierr = nf90_put_att(nc_fid,w_vid,"units","mm/s")

 !   ADDED BY JCH 04/07 TO OUTPUT SURFACE ELEVATION AND BOTTOM DEPTH
     !!=====elev=============================================!
     ierr = nf90_def_var(nc_fid,"elev",nf90_float,dynm1d,w_vid)
     ierr = nf90_put_att(nc_fid,w_vid,"long_name","surface elevation above particle")
     ierr = nf90_put_att(nc_fid,w_vid,"units","m")
      !!=====depth============================================!
     ierr = nf90_def_var(nc_fid,"depth",nf90_float,dynm1d,w_vid)
     ierr = nf90_put_att(nc_fid,w_vid,"long_name","bottom depth at particle")
     ierr = nf90_put_att(nc_fid,w_vid,"units","m")

 !   ADDED BY JCH 07/07 TO OUTPUT INWATER   
     ierr = nf90_def_var(nc_fid,"inwater",nf90_float,dynm1d,ind_vid)
     ierr = nf90_put_att(nc_fid,ind_vid,"long_name","particle inwater (1)")
     ierr = nf90_put_att(nc_fid,ind_vid,"units","")
  


 !   END ADDITION
     !--End definition section
     ierr = nf90_enddef(nc_fid)
    
     !--Write Particle label
     ALLOCATE(temp(npts))
     temp(:) = float(label(:))
     CALL putsvar(nc_fid,len_trim('label'),'label',npts,temp)
     DEALLOCATE(temp)

     !--Close file
     ierr = nf90_close(nc_fid)

  ENDIF
  
  ierr = nf90_open(trim(infile),nf90_write,nc_fid)
  IF(ierr /=nf90_noerr)THEN
     WRITE(*,*)'ERROR OPENING ',trim(infile)
     WRITE(*,*)trim(nf90_strerror(ierr))
     stop
  END IF

  !write time to file
  ALLOCATE(temp(1))
  temp(1)=time
  CALL putdvar(nc_fid,len_trim('time'),'time',1,temp,nt)
  DEALLOCATE(temp)

  !write indomain to file
  ALLOCATE(temp(npts))
  temp(:) = float(indomain(:))
  CALL putdvar(nc_fid,len_trim('indomain'),'indomain',npts,temp,nt)
  DEALLOCATE(temp)

  !write position and velocity to file
  CALL putdvar(nc_fid,len_trim('x'),'x',npts,xp,nt)
  CALL putdvar(nc_fid,len_trim('y'),'y',npts,yp,nt)
  CALL putdvar(nc_fid,len_trim('z'),'z',npts,zp,nt)
  CALL putdvar(nc_fid,len_trim('u'),'u',npts,up,nt)
  CALL putdvar(nc_fid,len_trim('v'),'v',npts,vp,nt)
  CALL putdvar(nc_fid,len_trim('omega'),'omega',npts,wp,nt)

  ! ADDED BY JHC 04/07 TO WRITE SURFACE ELEVATION AND BOTTOM DEPTH TO FILE
  CALL putdvar(nc_fid,len_trim('elev'),'elev',npts,ep,nt)
  CALL putdvar(nc_fid,len_trim('depth'),'depth',npts,hp,nt)
  
  !   ADDED BY JCH 07/07 TO OUTPUT INWATER  
  ALLOCATE(temp(npts))
  temp(:) = float(inwater(:))
  CALL putdvar(nc_fid,len_trim('inwater'),'inwater',npts,temp,nt)
  DEALLOCATE(temp)

  ! END ADDITION

  !close file
  ierr = nf90_close(nc_fid)
!  write(*,'(i4)') IERR
END SUBROUTINE ncd_write