diff --git a/postproc/comparevalww3board.py b/postproc/comparevalww3board.py
index 15ff3537a..79b777de8 100755
--- a/postproc/comparevalww3board.py
+++ b/postproc/comparevalww3board.py
@@ -155,11 +155,11 @@
 	# WAVEWATCH III --------------------------
 	if it==0:
 		fname=wrun1+"/ww3gefs."+cdate+"_spec.nc"
-		spm_time,lixo,lixo,lixo,spm_freq,spm_freq1,spm_freq2,spm_dfreq,spm_pspec1,spm_dmspec1,spm_dire,spm_dspec1 = wread.spec_ww3(fname,stname)
+		spm_time,lixo,lixo,lixo,spm_freq,spm_freq1,spm_freq2,spm_dfreq,spm_pspec1,spm_dmspec1,spm_dire,spm_dspec1,lixo,lixo = wread.spec_ww3(fname,stname)
 		del fname
 	else:
 		fname=wrun1+"/ww3gefs."+cdate+"_spec.nc"
-		aspm_time,lixo,lixo,lixo,lixo,lixo,lixo,lixo,aspm_pspec1,aspm_dmspec1,lixo,aspm_dspec1 = wread.spec_ww3(fname,stname)
+		aspm_time,lixo,lixo,lixo,lixo,lixo,lixo,lixo,aspm_pspec1,aspm_dmspec1,lixo,aspm_dspec1,lixo,lixo = wread.spec_ww3(fname,stname)
 		del fname
 		spm_time = np.append(spm_time,aspm_time)
 		spm_pspec1 = np.append(spm_pspec1,aspm_pspec1,axis=0)
@@ -168,11 +168,11 @@
 
 	if it==0:
 		fname=wrun2+"/ww3gefs."+cdate+"_spec.nc"
-		lixo,lixo,lixo,lixo,lixo,lixo,lixo,lixo,spm_pspec2,spm_dmspec2,lixo,spm_dspec2 = wread.spec_ww3(fname,stname)
+		lixo,lixo,lixo,lixo,lixo,lixo,lixo,lixo,spm_pspec2,spm_dmspec2,lixo,spm_dspec2,lixo,lixo = wread.spec_ww3(fname,stname)
 		del fname
 	else:
 		fname=wrun2+"/ww3gefs."+cdate+"_spec.nc"
-		lixo,lixo,lixo,lixo,lixo,lixo,lixo,lixo,aspm_pspec2,aspm_dmspec2,lixo,aspm_dspec2 = wread.spec_ww3(fname,stname)
+		lixo,lixo,lixo,lixo,lixo,lixo,lixo,lixo,aspm_pspec2,aspm_dmspec2,lixo,aspm_dspec2,lixo,lixo = wread.spec_ww3(fname,stname)
 		del fname
 		spm_pspec2 = np.append(spm_pspec2,aspm_pspec2,axis=0)
 		spm_dmspec2 = np.append(spm_dmspec2,aspm_dmspec2,axis=0)
diff --git a/postproc/wread.py b/postproc/wread.py
index 4a330640b..c96e687b9 100755
--- a/postproc/wread.py
+++ b/postproc/wread.py
@@ -354,6 +354,6 @@ def spec_ww3(*args):
 				d1sp[t,il]=np.float(aux)
 				del a,b,aux
 
-	return mdate,mtime,lat,lon,freq,freq1,freq2,dfreq,pwst,d1sp,dire,dspec
+	return mdate,mtime,lat,lon,freq,freq1,freq2,dfreq,pwst,d1sp,dire,dspec,wnds,wndd
 
 
diff --git a/validation/modelBuoy_collocation.py b/validation/modelBuoy_collocation.py
new file mode 100755
index 000000000..17a7c2224
--- /dev/null
+++ b/validation/modelBuoy_collocation.py
@@ -0,0 +1,372 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+modelBuoy_collocation.py
+
+VERSION AND LAST UPDATE:
+ v1.0  05/11/2022
+ v1.1  05/18/2022
+
+PURPOSE:
+ Collocation/pairing ww3 point output results with wave buoys.
+ Matchups of ww3 results and buoy data are generated, for the same
+  points (lat/lon) and time.
+
+USAGE:
+ This code is designed for ww3 hindcast or general simulations, not 
+  forecast with consecutive cycles (overlapped time) or ensemble forecasts.
+ Multiple files can be entered and appended, creating a continuous 
+  time array. Enter a list at ww3list.txt
+ WW3 netcdf results for point output tables (tab) is utilized.
+ It uses two public buoy databases, NDBC and Copernicus,
+  which (at least one) must have been previously downloaded. See
+  get_buoydata_copernicus.sh and retrieve_ndbc_nc.py
+ Edit ndbcp and copernp to set paths.
+ Users should confirm the buoys' names at the "f=nc.Dataset" lines below.
+ Python code can be run directly, without input arguments.
+ Additionally, users and enter two extra arguments, gridIndo and CycloneMap,
+  generated by prepGridMask.py and procyclmap.py, where the information
+  for the buoy's position (nearest grid point) will be extracted and
+  included in the output netcdf file.
+
+OUTPUT:
+ netcdf file ww3buoy_collocation_*.nc containing the matchups of buoy 
+  and ww3 data, for the stations (lat/lon) where both data sources 
+  are available.
+
+DEPENDENCIES:
+ See dependencies.py and the imports below.
+ ww3 table results in netcdf format (list of files ww3list.txt)
+ NDBC buoy data (see retrieve_ndbc_nc.py)
+ Copernicus buoy data (see get_buoydata_copernicus.sh)
+
+AUTHOR and DATE:
+ 05/11/2022: Ricardo M. Campos, first version.
+ 05/18/2022: Ricardo M. Campos, inclusion of gridInfo mask information,
+  with water depth, distance to coast, ocean names, forecast areas, and
+  cyclone information.
+
+PERSON OF CONTACT:
+ Ricardo M Campos: ricardo.campos@noaa.gov
+
+"""
+
+import numpy as np
+from matplotlib.mlab import *
+from pylab import *
+import xarray as xr
+import netCDF4 as nc
+import time
+from time import strptime
+from calendar import timegm
+import warnings; warnings.filterwarnings("ignore")
+# netcdf format
+fnetcdf="NETCDF4"
+
+# Paths
+# NDBC buoys
+# ndbcp="/data/buoys/NDBC/wparam"
+ndbcp="/work/noaa/marine/ricardo.campos/data/buoys/NDBC/ncformat/wparam"
+# Copernicus buoys
+# copernp="/data/buoys/Copernicus/wtimeseries"
+copernp="/work/noaa/marine/ricardo.campos/data/buoys/Copernicus/wtimeseries"
+# import os; os.system("ls /modelpath/*tab.nc > ww3list.txt &")
+wlist=np.loadtxt('ww3list.txt',dtype=str)
+if size(wlist)==1:
+	wlist=[wlist]
+
+# Options of including grid and cyclone information
+gridinfo=np.int(0); cyclonemap=np.int(0)
+if len(sys.argv) >= 2 :
+	gridinfo=np.str(sys.argv[1])
+if len(sys.argv) >= 3:
+	cyclonemap=np.str(sys.argv[2])
+if len(sys.argv) > 3:
+	sys.exit(' Too many inputs')
+
+# READ DATA
+print(" ")
+if gridinfo!=0:
+	f=nc.Dataset(gridinfo)
+	mlat=np.array(f.variables['latitude'][:]); mlon=np.array(f.variables['longitude'][:])
+	mask=f.variables['mask'][:,:]; distcoast=f.variables['distcoast'][:,:]; depth=f.variables['depth'][:,:]
+	oni=f.variables['GlobalOceansSeas'][:,:]; hsmz=f.variables['HighSeasMarineZones'][:,:]
+	ocnames=f.variables['names_GlobalOceansSeas'][:]; hsmznames=f.variables['names_HighSeasMarineZones'][:] 
+	f.close(); del f
+	print(" GridInfo Ok. "+gridinfo)
+
+	if cyclonemap!=0:
+		fcy=nc.MFDataset(cyclonemap, aggdim='time')
+		clat=np.array(fcy.variables['lat'][:]); clon=np.array(fcy.variables['lon'][:])
+		cmap=fcy.variables['cmap']; ctime=fcy.variables['time'][:]
+		cinfo=np.str(fcy.info); cinfo=np.array(np.str(cinfo).split(':')[1].split(';'))
+		if np.array_equal(clat,mlat)==True & np.array_equal(clon,mlon)==True: 
+			print(" CycloneMap Ok. "+cyclonemap)
+		else:
+			sys.exit(' Error: Cyclone grid and Mask grid are different.')
+
+# Model
+for t in range(0,size(wlist)):
+	try:
+		f=nc.Dataset(np.str(wlist[t]))
+	except:
+		print(" Cannot open "+wlist[t])
+	else:
+		# list of station/buoy names
+		if t==0:
+			auxstationname=f.variables['station_name'][:,:]; stname=[]
+			for i in range(0,auxstationname.shape[0]):
+				stname=np.append(stname,"".join(np.array(auxstationname[i,:]).astype('str')))
+
+		ahs = np.array(f.variables['hs'][:,:]).T
+		atm = np.array(f.variables['tr'][:,:]).T
+		adm = np.array(f.variables['th1m'][:,:]).T
+		at = np.array(f.variables['time'][:]*24*3600 + timegm( strptime(np.str(f.variables['time'].units).split(' ')[2][0:4]+'01010000', '%Y%m%d%H%M') )).astype('double')
+		f.close(); del f
+		if t==0:
+			mhs=np.copy(ahs)
+			mtm=np.copy(atm)
+			mdm=np.copy(adm)
+			mtime=np.copy(at)
+		else:
+			mhs=np.append(mhs,ahs,axis=1)
+			mtm=np.append(mtm,atm,axis=1)
+			mdm=np.append(mdm,adm,axis=1)
+			mtime=np.append(mtime,at)
+
+		del ahs,atm,adm,at
+
+print(" Read WW3 data OK. Start building the matchups model/buoy ...")
+
+# Buoys
+bhs=np.zeros((size(stname),size(mtime)),'f')*np.nan
+btm=np.zeros((size(stname),size(mtime)),'f')*np.nan
+bdm=np.zeros((size(stname),size(mtime)),'f')*np.nan
+lat=np.zeros(size(stname),'f')*np.nan; lon=np.zeros(size(stname),'f')*np.nan
+# help reading NDBC buoys, divided by year
+yrange=np.array(np.arange(time.gmtime(mtime.min())[0],time.gmtime(mtime.min())[0]+1,1)).astype('int')
+# loop buoys
+for b in range(0,size(stname)):
+
+	ahs=[]
+	try:
+
+		ahs=[];atm=[];adm=[];atime=[]
+		for y in yrange:
+
+			f=nc.Dataset(ndbcp+"/"+stname[b]+"h"+repr(y)+".nc")
+			ahs = np.append(ahs,f.variables['wave_height'][:,0,0])
+			atm = np.append(atm,f.variables['average_wpd'][:,0,0])
+			adm = np.append(adm,f.variables['mean_wave_dir'][:,0,0])
+			atime = np.append(atime,np.array(f.variables['time'][:]).astype('double'))
+			lat[b] = f.variables['latitude'][:]; lon[b] = f.variables['longitude'][:]
+			f.close(); del f
+
+	except:
+		try:
+			f=nc.Dataset(copernp+"/GL_TS_MO_"+stname[b]+".nc")
+			ahs = np.nanmean(f.variables['VHM0'][:,:],axis=1)
+			atm = np.nanmean(f.variables['VTM02'][:,:],axis=1)
+			adm = np.nanmean(f.variables['VMDR'][:,:],axis=1)
+			atime = np.array(f.variables['TIME'][:]*24*3600 + timegm( strptime('195001010000', '%Y%m%d%H%M') )).astype('double')
+			lat[b] = np.nanmean(f.variables['LATITUDE'][:]); lon[b] = np.nanmean(f.variables['LONGITUDE'][:])
+			f.close(); del f
+		except:
+			ahs=[]
+
+	else:
+		if size(ahs>0):
+			c=0
+			for t in range(0,size(mtime)):
+				indt=np.where(np.abs(atime-mtime[t])<1800.)
+				if size(indt)>0:
+					if np.any(ahs[indt[0]].mask==False):
+						bhs[b,t] = np.nanmean(ahs[indt[0]][ahs[indt[0]].mask==False])
+						c=c+1
+					if np.any(atm[indt[0]].mask==False):
+						btm[b,t] = np.nanmean(atm[indt[0]][atm[indt[0]].mask==False])
+					if np.any(adm[indt[0]].mask==False):
+						bdm[b,t] = np.nanmean(adm[indt[0]][adm[indt[0]].mask==False])
+
+					del indt
+
+			# print("counted "+repr(c)+" at "+stname[b])
+
+	print("done "+stname[b])
+	del ahs
+
+# Simple quality-control (range)
+ind=np.where((bhs>30.)|(bhs<0.0))
+if size(ind)>0:
+	bhs[ind]=np.nan; del ind
+
+ind=np.where((btm>40.)|(btm<0.0))
+if size(ind)>0:
+	btm[ind]=np.nan; del ind
+
+ind=np.where((bdm>360.)|(bdm<-180.))
+if size(ind)>0:
+	bdm[ind]=np.nan; del ind
+
+ind=np.where((mhs>30.)|(mhs<0.0))
+if size(ind)>0:
+	mhs[ind]=np.nan; del ind
+
+ind=np.where((mtm>40.)|(mtm<0.0))
+if size(ind)>0:
+	mtm[ind]=np.nan; del ind
+
+ind=np.where((mdm>360.)|(mdm<-180.))
+if size(ind)>0:
+	mdm[ind]=np.nan; del ind
+
+# Clean data. Select matchups only when model and buoy are available.
+ind=np.where( (np.isnan(lat)==False) & (np.isnan(lon)==False) & (np.isnan(np.nanmean(mhs,axis=1))==False) & (np.isnan(np.nanmean(bhs,axis=1))==False) )
+if size(ind)>0:
+	stname=np.array(stname[ind[0]])
+	lat=np.array(lat[ind[0]])
+	lon=np.array(lon[ind[0]])
+	mhs=np.array(mhs[ind[0],:])
+	mtm=np.array(mtm[ind[0],:])
+	mdm=np.array(mdm[ind[0],:])
+	bhs=np.array(bhs[ind[0],:])
+	btm=np.array(btm[ind[0],:])
+	bdm=np.array(bdm[ind[0],:])
+else:
+	sys.exit(' Error: No matchups Model/Buoy available.')
+
+print(" Matchups model/buoy complete. Total of "+repr(size(ind))+" points avaliable."); del ind
+
+# Processing grid and/or cyclone information
+if gridinfo!=0:
+	alon=np.copy(lon); alon[alon<0]=alon[alon<0]+360.
+	indgplat=[]; indgplon=[]
+	for i in range(0,lat.shape[0]):
+		# indexes nearest point.
+		indgplat = np.append(indgplat,np.where( abs(mlat-lat[i])==abs(mlat-lat[i]).min() )[0][0])
+		indgplon = np.append(indgplon,np.where( abs(mlon-alon[i])==abs(mlon-alon[i]).min() )[0][0])
+
+	indgplat=np.array(indgplat).astype('int'); indgplon=np.array(indgplon).astype('int')
+	pdistcoast=np.zeros(lat.shape[0],'f')*np.nan
+	pdepth=np.zeros(lat.shape[0],'f')*np.nan
+	poni=np.zeros(lat.shape[0],'f')*np.nan
+	phsmz=np.zeros(lat.shape[0],'f')*np.nan
+	for i in range(0,lat.shape[0]):
+		pdistcoast[i]=distcoast[indgplat[i],indgplon[i]]
+		pdepth[i]=depth[indgplat[i],indgplon[i]]
+		poni[i]=oni[indgplat[i],indgplon[i]]
+		phsmz[i]=hsmz[indgplat[i],indgplon[i]]
+
+	print(" Included Grid Information.")
+
+	# Excluding shallow water points too close to the coast (mask information not accurate)
+	ind=np.where( (np.isnan(pdistcoast)==False) & (np.isnan(pdepth)==False) )
+	if size(ind)>0:
+		stname=np.array(stname[ind[0]])
+		lat=np.array(lat[ind[0]])
+		lon=np.array(lon[ind[0]])
+		mhs=np.array(mhs[ind[0],:])
+		mtm=np.array(mtm[ind[0],:])
+		mdm=np.array(mdm[ind[0],:])
+		bhs=np.array(bhs[ind[0],:])
+		btm=np.array(btm[ind[0],:])
+		bdm=np.array(bdm[ind[0],:])
+		pdistcoast=np.array(pdistcoast[ind[0]])
+		pdepth=np.array(pdepth[ind[0]])
+		poni=np.array(poni[ind[0]])
+		phsmz=np.array(phsmz[ind[0]])
+	else:
+		sys.exit(' Error: No matchups Model/Buoy available after using grid mask.')
+
+	del ind
+
+	if cyclonemap!=0:
+		fcmap=np.zeros((lat.shape[0],mtime.shape[0]),'f')*np.nan
+		for t in range(0,size(mtime)):
+			# search for cyclone time index and cyclone map
+			indt=np.where(np.abs(ctime-mtime[t])<5400.)
+			if size(indt)>0:
+				for i in range(0,lat.shape[0]):
+					fcmap[i,t] = np.array(cmap[indt[0][0],indgplat[i],indgplon[i]])
+
+				del indt
+			else:
+				print(' Warning: No cyclone information for this time step: '+repr(t))
+
+			print(' Done cyclone analysis at time-step: '+repr(t))
+
+		ind=np.where(fcmap<0)
+		if size(ind)>0:
+			fcmap[ind]=np.nan
+
+		print(" Included Cyclone Information.")
+
+# Save netcdf output file 
+lon[lon>180.]=lon[lon>180.]-360.
+initime=np.str(time.gmtime(mtime.min())[0])+np.str(time.gmtime(mtime.min())[1]).zfill(2)+np.str(time.gmtime(mtime.min())[2]).zfill(2)+np.str(time.gmtime(mtime.min())[3]).zfill(2)
+fintime=np.str(time.gmtime(mtime.max())[0])+np.str(time.gmtime(mtime.max())[1]).zfill(2)+np.str(time.gmtime(mtime.max())[2]).zfill(2)+np.str(time.gmtime(mtime.max())[3]).zfill(2)
+ncfile = nc.Dataset('WW3.Buoy_'+initime+'to'+fintime+'.nc', "w", format=fnetcdf)
+ncfile.history="Matchups of WAVEWATCHIII point output (table) and NDBC and Copernicus Buoys. Total of "+repr(bhs[bhs>0.].shape[0])+" observations or pairs model/observation."
+# create  dimensions. 2 Dimensions
+ncfile.createDimension('buoypoints', bhs.shape[0] )
+ncfile.createDimension('time', bhs.shape[1] )
+if gridinfo!=0:
+	ncfile.createDimension('GlobalOceansSeas', ocnames.shape[0] )
+	ncfile.createDimension('HighSeasMarineZones', hsmznames.shape[0] )
+if cyclonemap!=0:
+	ncfile.createDimension('cycloneinfo', cinfo.shape[0] )
+	vcinfo = ncfile.createVariable('cycloneinfo',dtype('a25'),('cycloneinfo'))
+
+# create variables.
+vt = ncfile.createVariable('time',np.dtype('float64').char,('time'))
+vstname = ncfile.createVariable('buoyID',dtype('a25'),('buoypoints'))
+vlat = ncfile.createVariable('latitude',np.dtype('float32').char,('buoypoints'))
+vlon = ncfile.createVariable('longitude',np.dtype('float32').char,('buoypoints'))
+#
+vmhs = ncfile.createVariable('model_hs',np.dtype('float32').char,('buoypoints','time'))
+vmtm = ncfile.createVariable('model_tm',np.dtype('float32').char,('buoypoints','time'))
+vmdm = ncfile.createVariable('model_dm',np.dtype('float32').char,('buoypoints','time'))
+vbhs = ncfile.createVariable('obs_hs',np.dtype('float32').char,('buoypoints','time'))
+vbtm = ncfile.createVariable('obs_tm',np.dtype('float32').char,('buoypoints','time'))
+vbdm = ncfile.createVariable('obs_dm',np.dtype('float32').char,('buoypoints','time'))
+if gridinfo!=0:
+	vpdistcoast = ncfile.createVariable('distcoast',np.dtype('float32').char,('buoypoints')) 
+	vpdepth = ncfile.createVariable('depth',np.dtype('float32').char,('buoypoints')) 
+	vponi = ncfile.createVariable('GlobalOceansSeas',np.dtype('float32').char,('buoypoints'))
+	vocnames = ncfile.createVariable('names_GlobalOceansSeas',dtype('a25'),('GlobalOceansSeas'))
+	vphsmz = ncfile.createVariable('HighSeasMarineZones',np.dtype('float32').char,('buoypoints')) 
+	vhsmznames = ncfile.createVariable('names_HighSeasMarineZones',dtype('a25'),('HighSeasMarineZones'))
+if cyclonemap!=0:
+	vcmap = ncfile.createVariable('cyclone',np.dtype('float32').char,('buoypoints','time'))
+
+# Assign units
+vlat.units = 'degrees_north' ; vlon.units = 'degrees_east'
+vt.units = 'seconds since 1970-01-01T00:00:00+00:00'
+vmhs.units='m'; vbhs.units='m'
+vmtm.units='s'; vbtm.units='s'
+vmdm.units='degrees'; vbdm.units='degrees'
+if gridinfo!=0:
+	vpdepth.units='m'; vpdistcoast.units='km'
+
+# Allocate Data
+vt[:]=mtime[:]; vstname[:]=stname[:]
+vlat[:] = lat[:]; vlon[:] = lon[:]
+vmhs[:,:]=mhs[:,:]
+vmtm[:,:]=mtm[:,:]
+vmdm[:,:]=mdm[:,:]
+vbhs[:,:]=bhs[:,:]
+vbtm[:,:]=btm[:,:]
+vbdm[:,:]=bdm[:,:]
+if gridinfo!=0:
+	vpdistcoast[:]=pdistcoast[:]
+	vpdepth[:]=pdepth[:]
+	vponi[:]=poni[:]; vocnames[:] = ocnames[:]
+	vphsmz[:]=phsmz[:]; vhsmznames[:] = hsmznames[:]
+if cyclonemap!=0:
+	vcmap[:,:]=fcmap[:,:]; vcinfo[:] = cinfo[:]
+#
+ncfile.close()
+print(' ')
+print('Done. Netcdf ok. New file saved: WW3.Buoy_'+initime+'to'+fintime+'.nc')
+
diff --git a/validation/modelBuoy_collocation_GEFSforecast.py b/validation/modelBuoy_collocation_GEFSforecast.py
index 5b197d285..6dcd7c40e 100755
--- a/validation/modelBuoy_collocation_GEFSforecast.py
+++ b/validation/modelBuoy_collocation_GEFSforecast.py
@@ -356,10 +356,10 @@
 vmhs = ncfile.createVariable('model_hs',np.dtype('float32').char,('ensmembers','buoypoints','fcycletime','forecastleadtime'))
 vmtm = ncfile.createVariable('model_tm',np.dtype('float32').char,('ensmembers','buoypoints','fcycletime','forecastleadtime'))
 vmdm = ncfile.createVariable('model_dm',np.dtype('float32').char,('ensmembers','buoypoints','fcycletime','forecastleadtime'))
-vbwnd = ncfile.createVariable('buoy_wsp',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
-vbhs = ncfile.createVariable('buoy_hs',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
-vbtm = ncfile.createVariable('buoy_tm',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
-vbdm = ncfile.createVariable('buoy_dm',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
+vbwnd = ncfile.createVariable('obs_wsp',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
+vbhs = ncfile.createVariable('obs_hs',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
+vbtm = ncfile.createVariable('obs_tm',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
+vbdm = ncfile.createVariable('obs_dm',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
 vpdfc = ncfile.createVariable('distcoast',np.dtype('float32').char,('buoypoints'))
 vpbat = ncfile.createVariable('depth',np.dtype('float32').char,('buoypoints'))
 vpmask = ncfile.createVariable('mask',np.dtype('float32').char,('buoypoints'))
diff --git a/validation/modelBuoy_collocation_forecast.py b/validation/modelBuoy_collocation_forecast.py
index 541a6d612..ce3d4da0d 100755
--- a/validation/modelBuoy_collocation_forecast.py
+++ b/validation/modelBuoy_collocation_forecast.py
@@ -206,9 +206,9 @@
 vmhs = ncfile.createVariable('model_hs',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
 vmtm = ncfile.createVariable('model_tm',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
 vmdm = ncfile.createVariable('model_dm',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
-vbhs = ncfile.createVariable('buoy_hs',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
-vbtm = ncfile.createVariable('buoy_tm',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
-vbdm = ncfile.createVariable('buoy_dm',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
+vbhs = ncfile.createVariable('obs_hs',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
+vbtm = ncfile.createVariable('obs_tm',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
+vbdm = ncfile.createVariable('obs_dm',np.dtype('float32').char,('buoypoints','fcycletime','forecastleadtime'))
 # Assign units
 vlat.units = 'degrees_north' ; vlon.units = 'degrees_east'
 vt.units = 'seconds since 1970-01-01T00:00:00+00:00'
diff --git a/validation/modelBuoy_collocation_hindcast.py b/validation/modelBuoy_collocation_hindcast.py
deleted file mode 100755
index 6cd5c52dc..000000000
--- a/validation/modelBuoy_collocation_hindcast.py
+++ /dev/null
@@ -1,221 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-"""
-modelBuoy_collocation_hindcast.py
-
-VERSION AND LAST UPDATE:
- v1.0  04/04/2022
-
-PURPOSE:
- Collocation/pairing ww3 point output results with wave buoys.
- Matchups of ww3 results and buoy data are generated, for the same
-  points (lat/lon) and time.
-
-USAGE:
- This code is designed for ww3 hindcast or general simulations, not 
-  forecast with consecutive cycles or ensemble forecasts.
- Multiple files are entered and appended, creating a continuous 
-  time array. Enter a list at ww3list.txt
- Ww3 netcdf results for point output tables (tab) is utilized.
- It uses two public buoy databases, NDBC and Copernicus,
-  which (at least one) must have been previously downloaded. See
-  get_buoydata_copernicus.sh and retrieve_ndbc_nc.py
- Edit mpath, ndbcp, and copernp to set paths.
- Users have to confirm the buoys' names at the "f=nc.Dataset" lines below.
- Python code can be run directly, without input arguments.
-
-OUTPUT:
- netcdf file ww3buoy_collocation_*.nc containing the matchups of buoy 
-  and ww3 data, for the stations (lat/lon) where both data sources 
-  are available.
-
-DEPENDENCIES:
- See dependencies.py and the imports below.
- ww3 table results in netcdf format (list of files ww3list.txt)
- NDBC buoy data (see retrieve_ndbc_nc.py)
- Copernicus buoy data (see get_buoydata_copernicus.sh)
-
-AUTHOR and DATE:
- 04/04/2022: Ricardo M. Campos, first version.
-
-PERSON OF CONTACT:
- Ricardo M Campos: ricardo.campos@noaa.gov
-
-"""
-
-import numpy as np
-from matplotlib.mlab import *
-from pylab import *
-import xarray as xr
-import netCDF4 as nc
-import time
-from time import strptime
-from calendar import timegm
-import warnings; warnings.filterwarnings("ignore")
-# netcdf format
-fnetcdf="NETCDF4"
-
-# Paths
-# NDBC buoys
-ndbcp="/data/buoys/NDBC/wparam"
-# Copernicus buoys
-copernp="/data/buoys/Copernicus/wtimeseries"
-# import os; os.system("ls /modelpath/*tab.nc > ww3list.txt &")
-wlist=list(np.loadtxt('ww3list.txt',dtype=str))
-
-# Read Data
-# Model
-for t in range(0,size(wlist)):
-	try:
-		ds = xr.open_dataset(wlist[t]); f=nc.Dataset(wlist[t])
-	except:
-		print(" Cannot open "+wlist[t])
-	else:
-		# list of station/buoy names
-		if t==0:
-			auxstationname=ds['station_name'].values[:,:]; stname=[]
-			for i in range(0,auxstationname.shape[0]):
-				stname=np.append(stname,"".join(np.array(auxstationname[i,:]).astype('str')))
-
-		ahs = np.array(ds['hs'].values[:,:]).T
-		atm = np.array(ds['tr'].values[:,:]).T
-		adm = np.array(ds['th1m'].values[:,:]).T
-		at = np.array(f.variables['time'][:]*24*3600 + timegm( strptime(np.str(f.variables['time'].units).split(' ')[2][0:4]+'01010000', '%Y%m%d%H%M') )).astype('double')
-		ds.close();f.close(); del ds,f
-		if t==0:
-			mhs=np.copy(ahs)
-			mtm=np.copy(atm)
-			mdm=np.copy(adm)
-			mtime=np.copy(at)
-		else:
-			mhs=np.append(mhs,ahs,axis=1)
-			mtm=np.append(mtm,atm,axis=1)
-			mdm=np.append(mdm,adm,axis=1)
-			mtime=np.append(mtime,at)
-
-		del ahs,atm,adm,at
-
-# Buoys
-bhs=np.zeros((size(stname),size(mtime)),'f')*np.nan
-btm=np.zeros((size(stname),size(mtime)),'f')*np.nan
-bdm=np.zeros((size(stname),size(mtime)),'f')*np.nan
-lat=np.zeros(size(stname),'f')*np.nan; lon=np.zeros(size(stname),'f')*np.nan
-# help reading NDBC buoys, divided by year
-yrange=np.array(np.arange(time.gmtime(mtime.min())[0],time.gmtime(mtime.min())[0]+1,1)).astype('int')
-# loop buoys
-for b in range(0,size(stname)):
-
-	ahs=[]
-	try:
-
-		ahs=[];atm=[];adm=[];atime=[]
-		for y in yrange:
-
-			f=nc.Dataset(ndbcp+"/"+stname[b]+"h"+repr(y)+".nc")
-			ahs = np.append(ahs,f.variables['wave_height'][:,0,0])
-			atm = np.append(atm,f.variables['average_wpd'][:,0,0])
-			adm = np.append(adm,f.variables['mean_wave_dir'][:,0,0])
-			atime = np.append(atime,np.array(f.variables['time'][:]).astype('double'))
-			lat[b] = f.variables['latitude'][:]; lon[b] = f.variables['longitude'][:]
-			f.close(); del f
-
-	except:
-		try:
-			f=nc.Dataset(copernp+"/GL_TS_MO_"+stname[b]+".nc")
-			ahs = np.nanmean(f.variables['VHM0'][:,:],axis=1)
-			atm = np.nanmean(f.variables['VTM02'][:,:],axis=1)
-			adm = np.nanmean(f.variables['VMDR'][:,:],axis=1)
-			atime = np.array(f.variables['TIME'][:]*24*3600 + timegm( strptime('195001010000', '%Y%m%d%H%M') )).astype('double')
-			lat[b] = np.nanmean(f.variables['LATITUDE'][:]); lon[b] = np.nanmean(f.variables['LONGITUDE'][:])
-			f.close(); del f
-		except:
-			ahs=[]
-
-	else:
-		if size(ahs>0):
-			c=0
-			for t in range(0,size(mtime)):
-				indt=np.where(np.abs(atime-mtime[t])<1800.)
-				if size(indt)>0:
-					if np.any(ahs[indt[0]].mask==False):
-						bhs[b,t] = np.nanmean(ahs[indt[0]][ahs[indt[0]].mask==False])
-						c=c+1
-					if np.any(atm[indt[0]].mask==False):
-						btm[b,t] = np.nanmean(atm[indt[0]][atm[indt[0]].mask==False])
-					if np.any(adm[indt[0]].mask==False):
-						bdm[b,t] = np.nanmean(adm[indt[0]][adm[indt[0]].mask==False])
-
-					del indt
-
-			# print("counted "+repr(c)+" at "+stname[b])
-
-	print("done "+stname[b])
-	del ahs
-
-# Simple quality-control (range)
-ind=np.where((bhs>30.)|(bhs<0.0))
-if size(ind)>0:
-	bhs[ind]=np.nan; del ind
-
-ind=np.where((btm>40.)|(btm<0.0))
-if size(ind)>0:
-	btm[ind]=np.nan; del ind
-
-ind=np.where((bdm>360.)|(bdm<-180.))
-if size(ind)>0:
-	bdm[ind]=np.nan; del ind
-
-ind=np.where((mhs>30.)|(mhs<0.0))
-if size(ind)>0:
-	mhs[ind]=np.nan; del ind
-
-ind=np.where((mtm>40.)|(mtm<0.0))
-if size(ind)>0:
-	mtm[ind]=np.nan; del ind
-
-ind=np.where((mdm>360.)|(mdm<-180.))
-if size(ind)>0:
-	mdm[ind]=np.nan; del ind
-
-# Save netcdf output file 
-lon[lon>180.]=lon[lon>180.]-360.
-inidate=np.str(time.gmtime(mtime.min())[0])+np.str(time.gmtime(mtime.min())[1]).zfill(2)+np.str(time.gmtime(mtime.min())[2]).zfill(2)+np.str(time.gmtime(mtime.min())[3]).zfill(2)
-findate=np.str(time.gmtime(mtime.max())[0])+np.str(time.gmtime(mtime.max())[1]).zfill(2)+np.str(time.gmtime(mtime.max())[2]).zfill(2)+np.str(time.gmtime(mtime.max())[3]).zfill(2)
-ncfile = nc.Dataset("ww3buoy_collocation_"+inidate+"to"+findate+".nc", "w", format=fnetcdf) 
-ncfile.history="Collocation of WW3 point output (table) and NDBC/Copernicus Buoys. Total of "+repr(bhs[bhs>0.].shape[0])+" observations or pairs model/observation."
-# create  dimensions. 2 Dimensions
-ncfile.createDimension('buoypoints', bhs.shape[0] )
-ncfile.createDimension('time', bhs.shape[1] )
-# create variables.
-vt = ncfile.createVariable('time',np.dtype('float64').char,('time'))
-vstname = ncfile.createVariable('buoyID',dtype('a25'),('buoypoints'))
-vlat = ncfile.createVariable('latitude',np.dtype('float32').char,('buoypoints'))
-vlon = ncfile.createVariable('longitude',np.dtype('float32').char,('buoypoints'))
-#
-vmhs = ncfile.createVariable('model_hs',np.dtype('float32').char,('buoypoints','time'))
-vmtm = ncfile.createVariable('model_tm',np.dtype('float32').char,('buoypoints','time'))
-vmdm = ncfile.createVariable('model_dm',np.dtype('float32').char,('buoypoints','time'))
-vbhs = ncfile.createVariable('buoy_hs',np.dtype('float32').char,('buoypoints','time'))
-vbtm = ncfile.createVariable('buoy_tm',np.dtype('float32').char,('buoypoints','time'))
-vbdm = ncfile.createVariable('buoy_dm',np.dtype('float32').char,('buoypoints','time'))
-# Assign units
-vlat.units = 'degrees_north' ; vlon.units = 'degrees_east'
-vt.units = 'seconds since 1970-01-01T00:00:00+00:00'
-vmhs.units='m'; vbhs.units='m'
-vmtm.units='s'; vbtm.units='s'
-vmdm.units='degrees'; vbdm.units='degrees'
-# Allocate Data
-vt[:]=mtime[:]; vstname[:]=stname[:]
-vlat[:] = lat[:]; vlon[:] = lon[:]
-vmhs[:,:]=mhs[:,:]
-vmtm[:,:]=mtm[:,:]
-vmdm[:,:]=mdm[:,:]
-vbhs[:,:]=bhs[:,:]
-vbtm[:,:]=btm[:,:]
-vbdm[:,:]=bdm[:,:]
-#
-ncfile.close()
-print(' ')
-print('netcdf ok ')
-
diff --git a/validation/modelSat_collocation.py b/validation/modelSat_collocation.py
index 89e5945f0..17aba2984 100755
--- a/validation/modelSat_collocation.py
+++ b/validation/modelSat_collocation.py
@@ -26,6 +26,7 @@
    years in separated files, you can enter CycloneMap_*.nc and all the
    files/years in that directory will be properly read.
  - Altimeter Gridded netcdf files (generated by gridSatGlobal_Altimeter.py)
+   that must be written into a list named satlist.txt. Each line contains
    that must be written into a list named satlist.txt. Each line containing
    one satellite mission.
  - WAVEWATCHIII netcdf field outputs, also saved in a list ww3list.txt
@@ -219,36 +220,38 @@
 	ncfile = nc.Dataset('WW3.Altimeter_'+initime+'to'+fintime+'.nc', "w", format=fnetcdf)
 	ncfile.history="Matchups of WAVEWATCHIII and AODN Altimeter data. Total of "+repr(size(ind))+" observations or pairs model/observation."
 	# create  dimensions. 2 Dimensions
+	ncfile.createDimension('index',ftime.shape[0])
 	ncfile.createDimension('time',ftime.shape[0])
 	ncfile.createDimension('satellite', sdname.shape[0] )
 	ncfile.createDimension('GlobalOceansSeas', ocnames.shape[0] )
 	ncfile.createDimension('HighSeasMarineZones', hsmznames.shape[0] )
 	ncfile.createDimension('cycloneinfo', cinfo.shape[0] )
 	# create variables.
-	vt = ncfile.createVariable('time',np.dtype('float64').char,('time'))
-	vmonth = ncfile.createVariable('month',np.dtype('int16').char,('time'))	
-	vlat = ncfile.createVariable('latitude',np.dtype('float32').char,('time'))
-	vlon = ncfile.createVariable('longitude',np.dtype('float32').char,('time'))
-	vdistcoast = ncfile.createVariable('distcoast',np.dtype('float32').char,('time'))
-	vdepth = ncfile.createVariable('depth',np.dtype('float32').char,('time'))
-	voni = ncfile.createVariable('GlobalOceansSeas',np.dtype('int16').char,('time'))
+	vt = ncfile.createVariable('time',np.dtype('float64').char,('index'))
+	vmonth = ncfile.createVariable('month',np.dtype('int16').char,('index'))	
+	vlat = ncfile.createVariable('latitude',np.dtype('float32').char,('index'))
+	vlon = ncfile.createVariable('longitude',np.dtype('float32').char,('index'))
+	vdistcoast = ncfile.createVariable('distcoast',np.dtype('float32').char,('index'))
+	vdepth = ncfile.createVariable('depth',np.dtype('float32').char,('index'))
+	voni = ncfile.createVariable('GlobalOceansSeas',np.dtype('int16').char,('index'))
 	vocnames = ncfile.createVariable('names_GlobalOceansSeas',dtype('a25'),('GlobalOceansSeas'))
-	vhsmz = ncfile.createVariable('HighSeasMarineZones',np.dtype('int16').char,('time'))
+	vhsmz = ncfile.createVariable('HighSeasMarineZones',np.dtype('int16').char,('index'))
 	vhsmznames = ncfile.createVariable('names_HighSeasMarineZones',dtype('a25'),('HighSeasMarineZones'))
-	vcmap = ncfile.createVariable('cyclone',np.dtype('int16').char,('time'))
+	vcmap = ncfile.createVariable('cyclone',np.dtype('int16').char,('index'))
 	vcinfo = ncfile.createVariable('cycloneinfo',dtype('a25'),('cycloneinfo'))
-	vsid = ncfile.createVariable('satelliteID',np.dtype('int16').char,('time'))
+	vsid = ncfile.createVariable('satelliteID',np.dtype('int16').char,('index'))
 	vsdname = ncfile.createVariable('names_satellite',dtype('a25'),('HighSeasMarineZones'))
 	# results
-	vwhs = ncfile.createVariable('hs_model',np.dtype('float32').char,('time'))
-	vwwnd = ncfile.createVariable('wnd_model',np.dtype('float32').char,('time'))
-	vshs = ncfile.createVariable('hs_obsv',np.dtype('float32').char,('time'))
-	vswnd = ncfile.createVariable('wnd_obsv',np.dtype('float32').char,('time'))
+	vwhs = ncfile.createVariable('model_hs',np.dtype('float32').char,('index'))
+	vwwnd = ncfile.createVariable('model_wnd',np.dtype('float32').char,('index'))
+	vshs = ncfile.createVariable('obs_hs',np.dtype('float32').char,('index'))
+	vswnd = ncfile.createVariable('obs_wnd',np.dtype('float32').char,('index'))
 	# Assign units
 	vlat.units = 'degrees_north' ; vlon.units = 'degrees_east'
 	vt.units = 'seconds since 1970-01-01T00:00:00+00:00'
 	vwhs.units='m'; vshs.units='m'
 	vwwnd.units='m/s'; vswnd.units='m/s'
+	vdepth.units='m'; vdistcoast.units='km'
 	# Allocate Data
 	vt[:]=ftime[:]; vmonth[:]=fmonth[:]
 	vlat[:] = flat[:]; vlon[:] = flon[:]
diff --git a/validation/procyclmap.py b/validation/procyclmap.py
index 99746985f..e5dc3d7e9 100755
--- a/validation/procyclmap.py
+++ b/validation/procyclmap.py
@@ -70,7 +70,7 @@
 # READ mask
 f=nc.Dataset('gridInfo.nc')
 latm=f.variables['latitude'][:]; lonm=f.variables['longitude'][:]
-maskm=f.variables['mask'][:,:]; mdist=f.variables['distcoast'][:,:]
+maskm=f.variables['mask'][:,:];
 f.close(); del f
 # -------------