mod_station_timeseries.f90

Go to the documentation of this file.

MODULE mod_station_timeseries

  USE mod_time
  USE netcdf
  USE vars_wave
  IMPLICIT NONE
  SAVE
  
  LOGICAL out_station_timeseries_on
  CHARACTER(LEN=80) station_file
  CHARACTER(LEN=80) location_type
  LOGICAL out_elevation
  LOGICAL out_velocity_3d
  LOGICAL out_velocity_2d
  LOGICAL out_wind_velocity
  LOGICAL out_salt_temp
  CHARACTER(LEN=80) out_interval
  
  namelist /nml_station_timeseries/        &
       & out_station_timeseries_on,        &
       & station_file,                     &
       & location_type,                    &
       & out_elevation,                    &
       & out_velocity_3d,                  &
       & out_velocity_2d,                  &
       & out_wind_velocity,                &
       & out_salt_temp,                    &
       & out_interval
       
  INTEGER nsta       
  CHARACTER(LEN=20), ALLOCATABLE :: name_sta(:)
  REAL(sp), ALLOCATABLE :: lat_sta(:),lon_sta(:),h_sta(:)
  INTEGER, ALLOCATABLE  :: node_sta(:),element_sta(:),idummy(:)
  INTEGER, ALLOCATABLE  :: ntvegl(:),nbvegl(:,:)
  
  TYPE(time) :: interval_time_series, time_series
  INTEGER    :: kd_start
  TYPE(time) :: kdd,kdd1
  
!--Control Variables----------------------------------------------!
   integer,private :: out_cnt            !!counts number of outputs
   integer,private :: stck_cnt           !!counts number of outputs in each file
   character(len=120),private :: cdfname !!netcdf file name
!--NetCDF IDs----------------------------------------------------!

   !--NetCDF File 
   integer,private :: nc_ofid

   !--Dimensions
   integer,private :: station_did,clen_did
   integer,private :: siglay_did,siglev_did
   integer,private :: time_did

   !--Grid Variables
   integer,private :: x_s_vid,y_s_vid,lat_s_vid,lon_s_vid
   integer,private :: siglay_vid,siglev_vid

   !--Flow Variables 
   integer,private :: time_s_vid
   integer,private :: iint_vid
   integer,private :: u_s_vid
   integer,private :: v_s_vid
   integer,private :: ww_s_vid
   integer,private :: s1_s_vid
   integer,private :: t1_s_vid
   integer,private :: el_s_vid
   integer,private :: h_s_vid
   integer,private :: ua_s_vid
   integer,private :: va_s_vid
   integer,private :: uuwind_s_vid
   integer,private :: vvwind_s_vid
   integer,private :: atmpres_s_vid
   integer,private :: name_s_vid
   

   !--Info Variables
   character(len=120),public :: netcdf_timestring 


   INTERFACE putvar
      MODULE PROCEDURE putvar1d_int
      MODULE PROCEDURE putvar1d_real
      MODULE PROCEDURE putvar2d_int
      MODULE PROCEDURE putvar2d_real
   END INTERFACE

  CONTAINS
!====================================================================================
   SUBROUTINE station_name_list_initialize
   
   IMPLICIT NONE
   
   out_station_timeseries_on = .false.
   station_file         = "'none'"
   location_type             = "'node' or 'cell'"
   out_elevation             = .false.
   out_velocity_3d           = .false.
   out_velocity_2d           = .false.
   out_wind_velocity         = .false.
   out_salt_temp             = .false.
   out_interval              = "A length of time: 'seconds= ','days= ', or 'cycles= '"
   
   
   RETURN
   END SUBROUTINE station_name_list_initialize
!----------------------------------------------------------------------------------
   SUBROUTINE station_name_list_print
   USE control, ONLY : ipt
   
   IMPLICIT NONE
   
   WRITE(unit=ipt,nml=nml_station_timeseries)
   
   RETURN
   END SUBROUTINE station_name_list_print
!----------------------------------------------------------------------------------   
   SUBROUTINE station_name_list_read
   USE control, ONLY : casename,nmlunit
   USE mod_utils
   USE mod_set_time, ONLY : get_output_file_interval
   
   IMPLICIT NONE
   
   INTEGER :: IOS, I
   CHARACTER(LEN=120) :: FNAME
   CHARACTER(LEN=160) :: PATHNFILE
   
   IF(dbg_set(dbg_sbr)) &
         & WRITE(ipt,*) "Subroutine Begins: Read_Station_Name_List;"

   ios = 0

   fname = "./"//trim(casename)//"_run.nml"

   CALL fopen(nmlunit,trim(fname),'cfr')

   !READ NAME LIST FILE
    rewind(nmlunit)

   ! Read IO Information
   READ(unit=nmlunit, nml=nml_station_timeseries,iostat=ios)
   if(ios .NE. 0 ) Then
     if(dbg_set(dbg_log)) write(unit=ipt,nml=nml_station_timeseries)
     Call fatal_error("Can Not Read NameList NML_STATION_TIMESERIES from file: "//trim(fname))
   end if
   CLOSE(nmlunit)

!   CALL GET_OUTPUT_FILE_INTERVAL(TRIM(OUT_INTERVAL),INTERVAL_TIME_SERIES)

   RETURN
   END SUBROUTINE station_name_list_read
!----------------------------------------------------------------------------------   
   SUBROUTINE read_station_file
   USE control, ONLY : msr,input_dir,output_dir,casename,use_real_world_time,date_reference
   USE lims, ONLY : mgl,ngl,myid,msrid
   USE all_vars, ONLY : serial, par, nprocs, h, one_third, nvg,start_date
   USE mod_utils, ONLY : fatal_error
   USE mod_utils, ONLY : fopen
   USE mod_obcs, ONLY : gday1
   
   IMPLICIT NONE
   CHARACTER(LEN=120) :: FNAME
   CHARACTER(LEN=3)   :: NAC
   INTEGER            :: IZAJ_MAX,IOS,IDUMMY1,I
   LOGICAL            :: FEXIST
   
   REAL(SP), ALLOCATABLE :: FTEMP(:),FTEMPC(:)
   INTEGER               :: IDD,IMM,IYY,ICC,IHH,IMI,ISS 
   INTEGER               :: IDD1,IMM1,IYY1,ICC1,IHH1,IMI1,ISS1 
   INTEGER               :: KD_REFERENCE  
   
   out_cnt = 0

   ALLOCATE(ftemp(mgl))  ; ftemp = 0.0_sp
   IF(serial) ftemp = h

   IF(trim(location_type) == "cell" .OR. trim(location_type) == "CELL")THEN
     ALLOCATE(ftempc(ngl)) ;ftempc = 0.0_sp
     DO i = 1,ngl
       ftempc(i) = one_third*(ftemp(nvg(i,1))+ftemp(nvg(i,2))+ftemp(nvg(i,3)))
     END DO
   END IF    
     
   fname = trim(input_dir)//"/"//trim(station_file)
   CALL fopen(127,trim(fname),'cfr')

   IF(msr)THEN
     izaj_max = 0
     READ(127,*,iostat=ios)
     DO WHILE(.true.)
       READ(127,*,iostat=ios)idummy1
       IF(ios < 0)EXIT
       izaj_max = izaj_max + 1
     END DO
     nsta=izaj_max
     print *,'NSTA=',nsta
     ALLOCATE(name_sta(nsta))
     ALLOCATE(lat_sta(nsta))
     ALLOCATE(lon_sta(nsta))
     ALLOCATE(h_sta(nsta))
     ALLOCATE(idummy(nsta))
     IF(trim(location_type) == 'node' .OR. trim(location_type) == 'NODE')THEN
       ALLOCATE(node_sta(nsta))
     ELSE IF(trim(location_type) == 'cell' .OR. trim(location_type) == 'CELL')THEN  
       ALLOCATE(element_sta(nsta))
     ELSE
       CALL fatal_error("LOCATION_TYPE should be either node or cell")
     END IF  
     rewind(127)
 
     READ(127,*)
     DO i=1,nsta
       IF(trim(location_type) == "node" .OR. trim(location_type) == "NODE")THEN
         READ(127,*)idummy(i),lon_sta(i),lat_sta(i),node_sta(i),h_sta(i),name_sta(i)
         WRITE(6,*) idummy(i),lon_sta(i),lat_sta(i),node_sta(i),h_sta(i),ftemp(node_sta(i)),name_sta(i)
       ELSE IF(trim(location_type) == "cell" .OR. trim(location_type) == "CELL")THEN
         READ(127,*)idummy(i),lon_sta(i),lat_sta(i),element_sta(i),h_sta(i),name_sta(i)
         WRITE(6,*) idummy(i),lon_sta(i),lat_sta(i),element_sta(i),h_sta(i),ftempc(element_sta(i)),name_sta(i)
       END IF    
     ENDDO
     CLOSE(127)
     
   END IF

   DEALLOCATE(ftemp)
   IF(trim(location_type) == "cell" .OR. trim(location_type) == "CELL")DEALLOCATE(ftempc)
   
   IF(use_real_world_time)THEN
     READ(start_date(1:2),*) icc
     READ(start_date(3:4),*) iyy
     READ(start_date(6:7),*) imm
     READ(start_date(9:10),*) idd
     READ(start_date(12:13),*) ihh
     READ(start_date(15:16),*) imi
     READ(start_date(18:19),*) iss
     
     CALL gday1(idd,imm,iyy,icc,kd_start)

     IF(date_reference /= 'default')THEN
       READ(date_reference(1:2),*) icc1
       READ(date_reference(3:4),*) iyy1
       READ(date_reference(6:7),*) imm1
       READ(date_reference(9:10),*) idd1
       READ(date_reference(12:13),*) ihh1
       READ(date_reference(15:16),*) imi1
       READ(date_reference(18:19),*) iss1

       CALL gday1(idd1,imm1,iyy1,icc1,kd_reference)
     ELSE
       CALL gday1(17,11,58,18,kd_reference)
     END IF  

     kd_start = kd_start - kd_reference
     kdd%MJD = kd_start      
     kdd%MuSOD = ihh*3600.+imi*60.+iss
   END IF
   
   RETURN
   END SUBROUTINE read_station_file
!----------------------------------------------------------------------------------
   SUBROUTINE out_station_timeseries 
   
   USE all_vars
   USE netcdf
   use mod_nctools,only : handle_ncerr
   IMPLICIT NONE
   REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: UTMP,VTMP,T1TMP,S1TMP
   REAL(SP), ALLOCATABLE, DIMENSION(:)   :: UATMP,VATMP,UUWINDTMP,VVWINDTMP,ELTMP 
   REAL(SP) :: XTT,YTT
   INTEGER :: I1,I2,J 

   INTEGER  ::IZAJ_MAX,IZAJ_MIN,IZAJ,IZAJ_MAX_s,IZAJ_MIN_s,IZAJ_MAX_t,IZAJ_MIN_t
   INTEGER  ::kZAJ_MAX_s,kZAJ_MIN_s,kZAJ_MAX_t,kZAJ_MIN_t
   INTEGER  ::k,ierr
   REAL(SP) :: THOUR,THOUR1
   REAL(SP) :: STATMP(NSTA),STATMP1(NSTA),STATMPT(NSTA,KB),STATMPS(NSTA,KB)
   integer :: dims(1)
   real(sp), allocatable :: ftemp(:)
   integer :: VARID
   REAL(SP) :: KDD_TMP

!------------------------------------------------------------------------------!
!  WRITE TO FILES (SERIAL EXECUTION)                                           !
!------------------------------------------------------------------------------!
   IF(time_series > inttime) RETURN
   
   time_series = inttime + interval_time_series
   thour = dti*float(iint-istart+1)/3600.0_sp
   thour1 = dti*float(iint)/3600.0_sp
!  WRITE ELEVATION, SALINITY, TEMPERATURE, VELOCITY AND WIND DATA
   out_cnt = out_cnt + 1
   stck_cnt = stck_cnt + 1 
   if(out_cnt == 1) call write_netcdf_setup

   dims(1) = stck_cnt
   
!--Open File
   if(msr)then
     ierr = nf90_open(cdfname,nf90_write,nc_ofid)
     if(ierr /= nf90_noerr)then
       call handle_ncerr(ierr,"file open error")
     end if

!--Dump Time/IINT to File
     ierr    = nf90_put_var(nc_ofid,iint_vid,iint,start=dims)
     if(ierr /= nf90_noerr)then
       call handle_ncerr(ierr,"error writing variable to netcdf")
     end if

     IF(use_real_world_time)THEN
       kdd1%MJD = kdd%MJD + int((kdd%MuSOD/3600.+thour)/24.0)
       kdd1%MuSOD = kdd%MuSOD + thour * 3600 - int((kdd%MuSOD/3600.+thour)/24.0) * 24 * 3600
       kdd_tmp = kdd1%MJD + kdd1%MuSOD/86400.0

       ierr    = nf90_put_var(nc_ofid,time_s_vid,kdd_tmp,start=dims)
       if(ierr /= nf90_noerr)then
         call handle_ncerr(ierr,"error writing variable to netcdf")
       end if

     ELSE
       ierr    = nf90_put_var(nc_ofid,time_s_vid,thour1*3600.,start=dims)
       if(ierr /= nf90_noerr)then
         call handle_ncerr(ierr,"error writing variable to netcdf")
       end if
     END IF
   end if

!--Write Variables to File
   if(msr) write(ipt,*)'dumping to netcdf file: ',trim(cdfname),stck_cnt

     IF(out_elevation)THEN
       i1 = lbound(el,1) ; i2 = ubound(el,1)
       call putvar(i1,i2,m,mgl,1,1,"n",el,nc_ofid,el_s_vid,myid,nprocs&
            &,ipt, stck_cnt) 
     END IF

     IF(out_salt_temp)THEN 
       i1 = lbound(t1,1) ; i2 = ubound(t1,1)
       call putvar(i1,i2,m,mgl,kb,kb-1,"n",t1,nc_ofid,t1_s_vid,myid&
            &,nprocs,ipt, stck_cnt) 
       i1 = lbound(s1,1) ; i2 = ubound(s1,1)
       call putvar(i1,i2,m,mgl,kb,kb-1,"n",s1,nc_ofid,s1_s_vid,myid&
            &,nprocs,ipt, stck_cnt) 
     END IF

     IF(out_velocity_3d)THEN 
       i1 = lbound(u,1) ; i2 = ubound(u,1)
       call putvar(i1,i2,n,ngl,kb,kb-1,"e",u,nc_ofid,u_s_vid,myid&
            &,nprocs,ipt, stck_cnt) 
       i1 = lbound(v,1) ; i2 = ubound(v,1)
       call putvar(i1,i2,n,ngl,kb,kb-1,"e",v,nc_ofid,v_s_vid,myid&
            &,nprocs,ipt, stck_cnt) 
       i1 = lbound(ww,1) ; i2 = ubound(ww,1)
       call putvar(i1,i2,n,ngl,kb,kb-1,"e",ww,nc_ofid,ww_s_vid,myid&
            &,nprocs,ipt, stck_cnt) 
     END IF 

     IF(out_velocity_2d)THEN 
       allocate(ftemp(n))
       ftemp =ua(1:n)
       i1 = lbound(ftemp,1) ; i2 = ubound(ftemp,1)
       call putvar(i1,i2,n,ngl,1,1,"e",ftemp,nc_ofid,ua_s_vid,myid&
            &,nprocs,ipt, stck_cnt)
       deallocate(ftemp)
       allocate(ftemp(n))
       ftemp =va(1:n)
       i1 = lbound(ftemp,1) ; i2 = ubound(ftemp,1)
       call putvar(i1,i2,n,ngl,1,1,"e",ftemp,nc_ofid,va_s_vid,myid&
            &,nprocs,ipt, stck_cnt)
       deallocate(ftemp)
     END IF 

     IF(out_wind_velocity)THEN 
       allocate(ftemp(n))
       ftemp =uuwind(1:n)
       i1 = lbound(ftemp,1) ; i2 = ubound(ftemp,1)
       call putvar(i1,i2,n,ngl,1,1,"e",ftemp,nc_ofid,uuwind_s_vid,myid&
            &,nprocs,ipt, stck_cnt)
       deallocate(ftemp)
       allocate(ftemp(n))
       ftemp =vvwind(1:n)
       i1 = lbound(ftemp,1) ; i2 = ubound(ftemp,1)
       call putvar(i1,i2,n,ngl,1,1,"e",ftemp,nc_ofid,vvwind_s_vid,myid&
            &,nprocs,ipt, stck_cnt)
       deallocate(ftemp)
     END IF 



  ierr = nf90_close(nc_ofid)

   RETURN
   END SUBROUTINE out_station_timeseries
!==============================================================================|

!==============================================================================|
!  Write NetCDF Header and Static Variables                                    |
!==============================================================================|
   SUBROUTINE write_netcdf_setup 

   use all_vars
   use mod_clock
   use mod_nctools

   use netcdf
   use mod_types
   use mod_utils
   implicit none
   integer, dimension(3) :: dynm2de_lay,dynm2dn_lay
   integer, dimension(2) :: dynm2ds
   integer, dimension(1) :: stat2ds
   integer, dimension(2) :: stat2ds_lev,stat2ds_lay 
   integer, dimension(2) :: stat2ds_char
   integer, dimension(1) :: dynmtime
   character(len=100)    :: netcdf_convention
   character(len=100)    :: timestamp ,temp
   integer               :: i,j,ierr,i1,i2

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

!==============================================================================|
!  Set up Constants and Initialize Counters                                    |
!==============================================================================|
netcdf_timestring = 'seconds after 00:00:00'
!--Initialize Stack Count
   stck_cnt = 1

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

!--Time Stamp for History
   call get_timestamp(temp)
   timestamp = 'model started at: '//trim(temp)


!==============================================================================|
!  OPEN FILE AND DEFINE VARIABLES                                              |
!==============================================================================|
   if(msr)then

    cdfname = trim(output_dir)//trim(casename)//'_station_timeseries.nc'

!--Create File 
    ierr = nf90_create(path=cdfname,cmode=nf90_clobber,ncid=nc_ofid)
    if(ierr /= nf90_noerr)then
     call handle_ncerr(ierr,"file creation error")
    end if

!--Description of File Contents
    ierr = nf90_put_att(nc_ofid,nf90_global,"title"      ,trim(case_title))
    ierr = nf90_put_att(nc_ofid,nf90_global,"institution",trim(institution))
    ierr = nf90_put_att(nc_ofid,nf90_global,"source"     ,trim(fvcom_version))
    ierr = nf90_put_att(nc_ofid,nf90_global,"history"    ,trim(timestamp))
    ierr = nf90_put_att(nc_ofid,nf90_global,"references" ,trim(fvcom_website))
    ierr = nf90_put_att(nc_ofid,nf90_global,"Conventions",trim(netcdf_convention))

!--Define Fixed Model Dimensions 
    ierr = nf90_def_dim(nc_ofid,"siglay" ,kbm1   ,siglay_did )
    ierr = nf90_def_dim(nc_ofid,"siglev" ,kb     ,siglev_did )
    ierr = nf90_def_dim(nc_ofid,"station"   ,nsta      ,station_did   )
    ierr = nf90_def_dim(nc_ofid,"clen", 20, clen_did )

!--Define Unlimited Model Dimension
    ierr = nf90_def_dim(nc_ofid,"time"   ,nf90_unlimited,time_did)

!--Set Up Data Dimensioning - Static Vars
    stat2ds      = (/station_did/)            !!2d station vars
    stat2ds_lev  = (/station_did,siglev_did/)
    stat2ds_lay  = (/station_did,siglay_did/)
    stat2ds_char = (/clen_did,station_did/)

!--Set Up Data Dimensioning - Dynamic Vars 
    dynm2ds      = (/station_did,time_did/)            !!2d station vars
    dynm2de_lay  = (/station_did,siglay_did,time_did/) 
    dynm2dn_lay  = (/station_did,siglay_did,time_did/) 
    dynmtime     = (/time_did/)   

!--Define Station Name Variables and Attributes

    !!====Station Name (NAME_STA)  ===================!
    ierr = nf90_def_var(nc_ofid,"name_station",nf90_char,stat2ds_char,name_s_vid)
    ierr = nf90_put_att(nc_ofid,name_s_vid,"long_name","Station Name")   

!--Define Coordinate Variables and Attributes
    !!====X Grid Coordinate at Nodes (VX) (Meters)===========!
    ierr = nf90_def_var(nc_ofid,"x",nf90_float,stat2ds,x_s_vid)
    ierr = nf90_put_att(nc_ofid,x_s_vid,"long_name","station x-coordinate")
    ierr = nf90_put_att(nc_ofid,x_s_vid,"units","meters")

    !!====Y Grid Coordinate at Nodes (VY) (Meters)===========!
    ierr = nf90_def_var(nc_ofid,"y",nf90_float,stat2ds,y_s_vid)
    ierr = nf90_put_att(nc_ofid,y_s_vid,"long_name","station y-coordinate")
    ierr = nf90_put_att(nc_ofid,y_s_vid,"units","meters")

    !!====Longitudinal Coordinate at Nodes (LON) (degrees)===!
    ierr = nf90_def_var(nc_ofid,"lon",nf90_float,stat2ds,lon_s_vid)
    ierr = nf90_put_att(nc_ofid,lon_s_vid,"long_name","Longitude")
    ierr = nf90_put_att(nc_ofid,lon_s_vid,"standard_name","longitude")
    ierr = nf90_put_att(nc_ofid,lon_s_vid,"units","degrees_east")

    !!====Latitudinal  Coordinate at Nodes (LAT) (degrees)===!
    ierr = nf90_def_var(nc_ofid,"lat",nf90_float,stat2ds,lat_s_vid)
    ierr = nf90_put_att(nc_ofid,lat_s_vid,"long_name","Latitude")
    ierr = nf90_put_att(nc_ofid,lat_s_vid,"standard_name","latitude")
    ierr = nf90_put_att(nc_ofid,lat_s_vid,"units","degrees_north")
    ierr = nf90_put_att(nc_ofid,lat_s_vid,"grid","Bathymetry_Mesh")

   !!====Sigma Coordinate for Sigma Layers (ZZ)  (-)========!
   ierr = nf90_def_var(nc_ofid,"siglay",nf90_float,stat2ds_lay,siglay_vid)
   ierr = nf90_put_att(nc_ofid,siglay_vid,"long_name","Sigma Layers")
   ierr = nf90_put_att(nc_ofid,siglay_vid,"standard_name","ocean_sigma/general_coordinate")
   ierr = nf90_put_att(nc_ofid,siglay_vid,"positive","up")
   ierr = nf90_put_att(nc_ofid,siglay_vid,"valid_min","-1")
   ierr = nf90_put_att(nc_ofid,siglay_vid,"valid_max","0")
   ierr = nf90_put_att(nc_ofid,siglay_vid,"formula_terms","siglay:siglay eta:zeta depth:depth")

   !!====Sigma Coordinate for Sigma Levels (Z)   (-)========!
   ierr = nf90_def_var(nc_ofid,"siglev",nf90_float,stat2ds_lev,siglev_vid)
   ierr = nf90_put_att(nc_ofid,siglev_vid,"long_name","Sigma Levels")
   ierr = nf90_put_att(nc_ofid,siglev_vid,"standard_name","ocean_sigma/general_coordinate")
   ierr = nf90_put_att(nc_ofid,siglev_vid,"positive","up")
   ierr = nf90_put_att(nc_ofid,siglev_vid,"valid_min","-1")
   ierr = nf90_put_att(nc_ofid,siglev_vid,"valid_max","0")
   ierr = nf90_put_att(nc_ofid,siglev_vid,"formula_terms","siglev:siglev eta:zeta depth:depth")

!--Define Mesh Relevant Variables and Attributes

    !!====Bathymetry at Nodes (H) (meters)===================!
    ierr = nf90_def_var(nc_ofid,"h",nf90_float,stat2ds,h_s_vid)
    ierr = nf90_put_att(nc_ofid,h_s_vid,"long_name","Bathymetry")   
    ierr = nf90_put_att(nc_ofid,h_s_vid,"units","meters")
    ierr = nf90_put_att(nc_ofid,h_s_vid,"positive","down")
    ierr = nf90_put_att(nc_ofid,h_s_vid,"standard_name","depth")

!--Define Model Time Variables and Attributes    
    IF(use_real_world_time)THEN
      ierr = nf90_def_var(nc_ofid,"time",nf90_float,dynmtime,time_s_vid)
      ierr = nf90_put_att(nc_ofid,time_s_vid,"long_name","time")
      if(date_reference == 'default')then
        ierr = nf90_put_att(nc_ofid,time_s_vid,"units",trim("days since 1858-11-17 00:00:00"))
        ierr = nf90_put_att(nc_ofid,time_s_vid,"format",trim("modified julian day (MJD)"))
      else
        ierr = nf90_put_att(nc_ofid,time_s_vid,"units","days since "//trim(date_reference))
        ierr = nf90_put_att(nc_ofid,time_s_vid,"format",trim("defined reference date"))
      end if
!JQI      ierr = nf90_put_att(nc_ofid,time_s_vid,"calendar","none")
      ierr = nf90_put_att(nc_ofid,time_s_vid,"time_zone","UTC")

    ELSE
      ierr = nf90_def_var(nc_ofid,"time",nf90_float,dynmtime,time_s_vid)
      ierr = nf90_put_att(nc_ofid,time_s_vid,"long_name","time")
      ierr = nf90_put_att(nc_ofid,time_s_vid,"units",trim(netcdf_timestring))
!JQI      ierr = nf90_put_att(nc_ofid,time_s_vid,"calendar","none")
      ierr = nf90_put_att(nc_ofid,time_s_vid,"time_zone","none")
    END IF  

    ierr = nf90_def_var(nc_ofid,"iint",nf90_int,dynmtime,iint_vid)
    ierr = nf90_put_att(nc_ofid,iint_vid,"long_name","internal mode iteration number")

!--Define Time Dependent Flow Variables (selected by user from input file)
    if(out_velocity_3d)then
     ierr = nf90_def_var(nc_ofid,"u",nf90_float,dynm2de_lay,u_s_vid)
     ierr = nf90_put_att(nc_ofid,u_s_vid,"long_name","Eastward Water Velocity")
     ierr = nf90_put_att(nc_ofid,u_s_vid,"standard_name","eastward_sea_water_velocity")
     ierr = nf90_put_att(nc_ofid,u_s_vid,"units","meters s-1")
     ierr = nf90_put_att(nc_ofid,u_s_vid,"type","data")
     ierr = nf90_put_att(nc_ofid,u_s_vid,"coordinates","time siglay station")
       
     ierr = nf90_def_var(nc_ofid,"v",nf90_float,dynm2de_lay,v_s_vid)
     ierr = nf90_put_att(nc_ofid,v_s_vid,"long_name","Northward Water Velocity")
     ierr = nf90_put_att(nc_ofid,u_s_vid,"standard_name","northward_sea_water_velocity")
     ierr = nf90_put_att(nc_ofid,v_s_vid,"units","meters s-1")
     ierr = nf90_put_att(nc_ofid,v_s_vid,"type","data")
     ierr = nf90_put_att(nc_ofid,v_s_vid,"coordinates","time siglay station")

     ierr = nf90_def_var(nc_ofid,"ww",nf90_float,dynm2de_lay,ww_s_vid)
     ierr = nf90_put_att(nc_ofid,ww_s_vid,"long_name","Upward Water Velocity")
     ierr = nf90_put_att(nc_ofid,ww_s_vid,"units","meters s-1")
     ierr = nf90_put_att(nc_ofid,ww_s_vid,"type","data")
    end if

    if(out_velocity_2d)then
     ierr = nf90_def_var(nc_ofid,"ua",nf90_float,dynm2ds,ua_s_vid)
     ierr = nf90_put_att(nc_ofid,ua_s_vid,"long_name","Vertically Averaged x-velocity")
     ierr = nf90_put_att(nc_ofid,ua_s_vid,"units","meters s-1")
     ierr = nf90_put_att(nc_ofid,ua_s_vid,"type","data")
       
     ierr = nf90_def_var(nc_ofid,"va",nf90_float,dynm2ds,va_s_vid)
     ierr = nf90_put_att(nc_ofid,va_s_vid,"long_name","Vertically Averaged y-velocity")
     ierr = nf90_put_att(nc_ofid,va_s_vid,"units","meters s-1")
     ierr = nf90_put_att(nc_ofid,va_s_vid,"type","data")
    end if

    if(out_salt_temp)then
     ierr = nf90_def_var(nc_ofid,"temp",nf90_float,dynm2dn_lay,t1_s_vid)
     ierr = nf90_put_att(nc_ofid,t1_s_vid,"long_name","temperature")
     ierr = nf90_put_att(nc_ofid,t1_s_vid,"standard_name","sea_water_temperature")
     ierr = nf90_put_att(nc_ofid,t1_s_vid,"units","degrees_C")
     ierr = nf90_put_att(nc_ofid,t1_s_vid,"type","data")
     ierr = nf90_put_att(nc_ofid,t1_s_vid,"coordinates","time siglay station")

     ierr = nf90_def_var(nc_ofid,"salinity",nf90_float,dynm2dn_lay,s1_s_vid)
     ierr = nf90_put_att(nc_ofid,s1_s_vid,"long_name","salinity")
     ierr = nf90_put_att(nc_ofid,s1_s_vid,"standard_name","sea_water_salinity")
     ierr = nf90_put_att(nc_ofid,s1_s_vid,"units","1e-3")
     ierr = nf90_put_att(nc_ofid,s1_s_vid,"type","data")
     ierr = nf90_put_att(nc_ofid,s1_s_vid,"coordinates","time siglay station")
    end if

    if(out_elevation)then
     ierr = nf90_def_var(nc_ofid,"zeta",nf90_float,dynm2ds,el_s_vid)
     ierr = nf90_put_att(nc_ofid,el_s_vid,"long_name","Water Surface Elevation")
     ierr = nf90_put_att(nc_ofid,el_s_vid,"units","meters")
     ierr = nf90_put_att(nc_ofid,el_s_vid,"positive","up")
     ierr = nf90_put_att(nc_ofid,el_s_vid,"standard_name","sea_surface_height_above_geoid")
     ierr = nf90_put_att(nc_ofid,el_s_vid,"type","data")
     ierr = nf90_put_att(nc_ofid,el_s_vid,"coordinates","time station")
    end if

    if(out_wind_velocity)then
     ierr = nf90_def_var(nc_ofid,"uwind_speed",nf90_float,dynm2ds,uuwind_s_vid)
     ierr = nf90_put_att(nc_ofid,uuwind_s_vid,"long_name","Eastward wind velocity")
     ierr = nf90_put_att(nc_ofid,uuwind_s_vid,"units","(m/s)")
     ierr = nf90_put_att(nc_ofid,uuwind_s_vid,"standard_name","eastward wind")
     ierr = nf90_put_att(nc_ofid,uuwind_s_vid,"type","data")
     ierr = nf90_put_att(nc_ofid,uuwind_s_vid,"coordinates","time station")

     ierr = nf90_def_var(nc_ofid,"vwind_speed",nf90_float,dynm2ds,vvwind_s_vid)
     ierr = nf90_put_att(nc_ofid,vvwind_s_vid,"long_name","Northward wind velocity")
     ierr = nf90_put_att(nc_ofid,vvwind_s_vid,"units","(m/s)")
     ierr = nf90_put_att(nc_ofid,vvwind_s_vid,"standard_name","northward wind")
     ierr = nf90_put_att(nc_ofid,vvwind_s_vid,"type","data")
     ierr = nf90_put_att(nc_ofid,vvwind_s_vid,"coordinates","time station")
    end if


!--Exit Define Mode
    ierr = nf90_enddef(nc_ofid)
    ierr = nf90_close(nc_ofid)

   end if !(msr)

!==============================================================================|
!  WRITE VARIABLES TO FILE                                                     |
!==============================================================================|
   if(msr)then
     ierr = nf90_open(cdfname,nf90_write,nc_ofid)
     if(ierr /= nf90_noerr)then
       call handle_ncerr(ierr,"file open error")
     end if
   end if
   
   !!====Longitude at Nodes (LON) ==========================!
   i1 = lbound(lon,1) ; i2 = ubound(lon,1)
   call putvar(i1,i2,m,mgl,1,1,"n",lon,nc_ofid,lon_s_vid,myid&
        &,nprocs,ipt, stck_cnt)

   !!====Latitude  at Nodes (LAT) ==========================!
   i1 = lbound(lat,1) ; i2 = ubound(lat,1)
   call putvar(i1,i2,m,mgl,1,1,"n",lat,nc_ofid,lat_s_vid,myid&
        &,nprocs,ipt, stck_cnt) 

   !!====X Grid Coordinate at Nodes (VX)====================!
   i1 = lbound(vx,1) ; i2 = ubound(vx,1)
   call putvar(i1,i2,m,mgl,1,1,"n",vx+vxmin,nc_ofid,x_s_vid,myid,nprocs&
        &,ipt, stck_cnt) 

   !!====Y Grid Coordinate at Nodes (VY)====================!
   i1 = lbound(vy,1) ; i2 = ubound(vy,1)
   call putvar(i1,i2,m,mgl,1,1,"n",vy+vymin,nc_ofid,y_s_vid,myid,nprocs&
        &,ipt, stck_cnt) 

   !!====Bathymetry at Nodes (H)============================!
   i1 = lbound(h,1) ; i2 = ubound(h,1)
   call putvar(i1,i2,m,mgl,1,1,"n",h,nc_ofid,h_s_vid,myid,nprocs,ipt,&
        & stck_cnt) 

   !!====Sigma Layers (zz)==================================!
   i1 = lbound(zz,1) ; i2 = ubound(zz,1)
   call putvar(i1,i2,m,mgl,kb,kb-1,"n",zz,nc_ofid,siglay_vid,myid&
        &,nprocs,ipt, stck_cnt) 

   !!====Sigma Levels (z)==================================!
   i1 = lbound(z,1) ; i2 = ubound(z,1)
   call putvar(i1,i2,m,mgl,kb,kb,"n",z,nc_ofid,siglev_vid,myid,nprocs&
        &,ipt, stck_cnt) 

   !!====Station Name (NAME_STA)============================!
   if(msr)then
    i1 = lbound(name_sta,1) ; i2 = ubound(name_sta,1)
    call putvar_char(i1,i2,name_sta,nc_ofid,name_s_vid,myid,nprocs,ipt,&
        & stck_cnt) 
   end if   

!==============================================================================|
!  close the file                                                              |
!==============================================================================|

   if(msr) ierr = nf90_close(nc_ofid)

   return
   end subroutine write_netcdf_setup
!==============================================================================|
 !==============================================================================|
 !  Collect Data to Global Array and Write to Netcdf File                       |
 !==============================================================================|
 
 SUBROUTINE putvar1d_real(i1,i2,n1,n1gl,kt,k1,map_type,var,nc_fid,vid&
      &,myid,nprocs,ipt,stk)
   
   !------------------------------------------------------------------------------|
   implicit none
   integer, intent(in) :: i1,i2,n1,n1gl,kt,k1,nc_fid,vid,myid,nprocs&
        &,ipt,stk
   character(len=*),intent(in)   :: map_type
   real(sp), dimension(i1:i2) :: var
   
   real(sp), allocatable, dimension(:,:) :: temp

   allocate(temp(i1:i2,kt))
   temp(i1:i2,1)=var

   CALL putvar2d_real(i1,i2,n1,n1gl,kt,k1,map_type,temp,nc_fid,vid&
        &,myid,nprocs,ipt,stk)
   
   deallocate(temp)

 END SUBROUTINE putvar1d_real
   

 subroutine putvar2d_real(i1,i2,n1,n1gl,kt,k1,map_type,var,nc_fid,vid&
      &,myid,nprocs,ipt,stk)
!------------------------------------------------------------------------------|

   use all_vars, only : nvg
   use mod_nctools,only : handle_ncerr
   use mod_types
   use mod_utils, only : pstop
   implicit none
   integer, intent(in) :: i1,i2,n1,n1gl,kt,k1,nc_fid,vid,myid,nprocs&
        &,ipt, stk
   character(len=*),intent(in)   :: map_type
   real(sp), allocatable :: var(:,:)

   real(sp), allocatable, dimension(:,:) :: temp,gtemp
   integer :: ierr,k1m1,i,j
   integer, allocatable :: dims(:)
   

   k1m1 = k1 
   if(k1m1 == 1)then
     allocate(dims(2))
     dims(1) = 1 
     dims(2) = stk
   else
     allocate(dims(3))
     dims(1) = 1 
     dims(2) = 1 
     dims(3) = stk      
   end if
     

   if(map_type(1:1) /= "e" .and. map_type(1:1) /= "n")then
     write(ipt,*)'map_type input to putvar should be "e" OR "n"'
     call pstop
   end if

   if(nprocs==1)then
     allocate(temp(nsta,k1m1))  

     if(map_type(1:1) == 'n' .and. trim(location_type) == 'node')then
      do i=1,nsta
       temp(i,1:k1m1) = var(node_sta(i),1:k1m1)
      end do  
     else if(map_type(1:1) == 'e' .and. trim(location_type) == 'cell')then
      do i=1,nsta
       temp(i,1:k1m1) = var(element_sta(i),1:k1m1)
      end do  
     else if(map_type(1:1) == 'n' .and. trim(location_type) == 'cell')then
      do i=1,nsta
       temp(i,:) = 0.0_sp
       do j =1,3
        temp(i,1:k1m1) = temp(i,1:k1m1) + var(nvg(element_sta(i),j),1:k1m1)
       end do
       temp(i,:) = temp(i,:)/3.0_sp
      end do   
     else if(map_type(1:1) == 'e' .and. trim(location_type) == 'node')then
      do i=1,nsta
       temp(i,:) = 0.0_sp
       do j =1, ntvegl(node_sta(i))
        temp(i,1:k1m1) = temp(i,1:k1m1) + var(nbvegl(node_sta(i),j),1:k1m1)
       end do
       temp(i,:) = temp(i,:)/ntvegl(node_sta(i))
      end do   
     end if
   end if


!   if(myid /= 1) return
   if(myid == 1) then
     ierr = nf90_put_var(nc_fid,vid,temp,start=dims)
     if(ierr /= nf90_noerr)then
       call handle_ncerr(ierr,"error writing variable to netcdf")
     end if
   end if  
   deallocate(dims,temp)

   return
 end subroutine putvar2d_real
!==============================================================================|

 SUBROUTINE putvar1d_int(i1,i2,n1,n1gl,kt,k1,map_type,var,nc_fid,vid&
      &,myid,nprocs,ipt,stk )
   
   !------------------------------------------------------------------------------|
   implicit none
   integer, intent(in) :: i1,i2,n1,n1gl,kt,k1,nc_fid,vid,myid,nprocs&
        &,ipt, stk
   character(len=*),intent(in)   :: map_type
   INTEGER, dimension(i1:i2) :: var
   
   INTEGER, allocatable, dimension(:,:) :: temp
   
   allocate(temp(i1:i2,kt))
   temp(i1:i2,kt)= var
   
   call putvar2d_int(i1,i2,n1,n1gl,kt,k1,map_type,temp,nc_fid,vid&
        &,myid,nprocs,ipt, stk)
   
   deallocate(temp)
   
 END SUBROUTINE putvar1d_int
 
 subroutine putvar2d_int(i1,i2,n1,n1gl,kt,k1,map_type,var,nc_fid,vid&
      &,myid,nprocs,ipt, stk)
   
   !------------------------------------------------------------------------------|
   
   use all_vars, only : nvg
   use mod_nctools,only : handle_ncerr
   use mod_utils, only : pstop
   implicit none
   integer, intent(in) :: i1,i2,n1,n1gl,kt,k1,nc_fid,vid,myid,nprocs&
        &,ipt,stk
   character(len=*),intent(in)   :: map_type
   INTEGER, allocatable :: var(:,:)

   INTEGER, allocatable, dimension(:,:) :: temp,gtemp
   integer :: ierr,k1m1,i,j
   integer, allocatable :: dims(:)
   

   k1m1 = k1 
   if(k1m1 == 1)then
     allocate(dims(2))
     dims(1) = 1 
     dims(2) = stk
   else
     allocate(dims(3))
     dims(1) = 1 
     dims(2) = 1 
     dims(3) = stk 
   end if
     

   if(map_type(1:1) /= "e" .and. map_type(1:1) /= "n")then
     write(ipt,*)'map_type input to putvar should be "e" OR "n"'
     call pstop
   end if

   if(nprocs==1)then
     allocate(temp(nsta,k1m1))

     if(map_type(1:1) == 'n' .and. trim(location_type) == 'node')then
      do i=1,nsta
       temp(i,1:k1m1) = var(node_sta(i),1:k1m1)
      end do  
     else if(map_type(1:1) == 'e' .and. trim(location_type) == 'cell')then
      do i=1,nsta
       temp(i,1:k1m1) = var(element_sta(i),1:k1m1)
      end do  
     else if(map_type(1:1) == 'n' .and. trim(location_type) == 'cell')then
      do i=1,nsta
       temp(i,:) = 0.0_sp
       do j =1,3
        temp(i,1:k1m1) = temp(i,1:k1m1) + var(nvg(element_sta(i),j),1:k1m1)
       end do
       temp(i,:) = temp(i,:)/3.0_sp
      end do   
     else if(map_type(1:1) == 'e' .and. trim(location_type) == 'node')then
      do i=1,nsta
       temp(i,:) = 0.0_sp
       do j =1, ntvegl(node_sta(i))
        temp(i,1:k1m1) = temp(i,1:k1m1) + var(nbvegl(node_sta(i),j),1:k1m1)
       end do
       temp(i,:) = temp(i,:)/ntvegl(node_sta(i))
      end do   
     end if
   end if


!   if(myid /= 1) return
   if(myid == 1) then
     ierr = nf90_put_var(nc_fid,vid,temp,start=dims)
     if(ierr /= nf90_noerr)then
       call handle_ncerr(ierr,"error writing variable to netcdf")
     end if
   end if  
   deallocate(dims,temp)

   return
 end subroutine putvar2d_int
!==============================================================================|
!==============================================================================|

 SUBROUTINE putvar_char(i1,i2,var,nc_fid,vid,myid,nprocs,ipt,stk )
   
   !------------------------------------------------------------------------------|
   use mod_nctools,only : handle_ncerr
   implicit none
   integer, intent(in) :: i1,i2,nc_fid,vid,myid,nprocs,ipt, stk
   CHARACTER(LEN=20), dimension(i1:i2) :: var
   CHARACTER(LEN=20), allocatable,dimension(:,:) :: var_tmp
   integer :: ierr
   integer, allocatable :: dims(:)

   allocate(dims(2))
   dims(1) = 1 
   dims(2) = stk
   
   allocate(var_tmp(i1:i2,1))
   var_tmp(i1:i2,1) = var(i1:i2)

   ierr = nf90_put_var(nc_fid,vid,var_tmp,start=dims)
   if(ierr /= nf90_noerr)then
     call handle_ncerr(ierr,"error writing variable to netcdf")
   end if

   deallocate(dims)

   return
 end subroutine putvar_char
!==============================================================================|
!==============================================================================|
   SUBROUTINE triangle_grid_edge_gl

!==============================================================================|
!             DEFINE NTVEGL, NBVEGL                                          !
!                                                                              !
! ntvegl(1:mgl):           total number of the surrounding triangles               !
!                      connected to the given node                             !
! nbvegl(1:mgl, 1:ntvegl+1): the identification number of surrounding                !
!                      triangles with a common node (counted clockwise)        !
!==============================================================================|
   USE all_vars
   IMPLICIT NONE
   INTEGER I,J,NCNT,MX_NBR_ELEM_GL

!
!----DETERMINE MAX NUMBER OF SURROUNDING ELEMENTS------------------------------!
!
   mx_nbr_elem_gl = 0
   DO i=1,mgl
     ncnt = 0
     DO j=1,ngl
       IF( float(nvg(j,1)-i)*float(nvg(j,2)-i)*float(nvg(j,3)-i) == 0.0_sp) &
         ncnt = ncnt + 1
     END DO
     mx_nbr_elem_gl = max(mx_nbr_elem_gl,ncnt)
   END DO

!
!----ALLOCATE ARRAYS BASED ON MX_NBR_ELEM--------------------------------------!
! 
   ALLOCATE(ntvegl(mgl));                  ntvegl = 0
   ALLOCATE(nbvegl(mgl,mx_nbr_elem_gl+1)); nbvegl = 0
!
!--DETERMINE NUMBER OF SURROUNDING ELEMENTS FOR NODE I = NTVEGL(I)---------------!
!--DETERMINE NBVEGL - INDICES OF NEIGHBORING ELEMENTS OF NODE I------------------!
!
       
   DO i=1,mgl
     ncnt=0
     DO j=1,ngl
       IF (float(nvg(j,1)-i)*float(nvg(j,2)-i)*float(nvg(j,3)-i) == 0.0_sp)THEN
         ncnt = ncnt+1
         nbvegl(i,ncnt)=j
       END IF
     ENDDO
     ntvegl(i)=ncnt
   ENDDO

   RETURN
   END SUBROUTINE triangle_grid_edge_gl
!==============================================================================!
END MODULE mod_station_timeseries