PROGRAM MESO2GRADS

!	This is a program to read Oklahoma Mesonet
!	data and to write the information to
!	a binary station file suitable for use
!	with GrADS.
!	Christopher M. Godfrey		21 June 2004

!	Modified to read and write soil variables
!	Christopher M. Godfrey		30 August 2005

!	NOTE: GrADS (specifically the stnmap utility) is
!	very finicky about the binary data format.  This
!	program works, but only with these compilation options
!	on NSSLcluster1: "pgf90 -Mfreeform -pc 32 -byteswapio".
!	The executable is named meso2grads.exe.  Complete the
!	GrADS setup with "stnmap [-v] -i [descriptor].ctl".

!	This does NOT include FITTSTOWN!!!

	IMPLICIT NONE
	INTEGER, PARAMETER :: nstations=126,maxstations=500	!Number of stations in geomeso.tbl
	CHARACTER :: stid*8,station(nstations)*4,rstid*4,stidarr(maxstations)*4
	REAL :: lat(nstations),lon(nstations)
	REAL :: tair,wspd,wvec,wdsd,wssd,wmax,rain,pres,ta9m,ws2m
	INTEGER :: elev(nstations)
	INTEGER :: relh,wdir,srad
	INTEGER :: flag,noflag
	INTEGER :: match,index,soil
	INTEGER :: i,m,nlev
	DOUBLE PRECISION :: time
	CHARACTER :: date*12,datafile*52,binfile*45
	REAL :: wc05(maxstations),wc25(maxstations),wc60(maxstations),wc75(maxstations)
	REAL :: fw05(maxstations),fw25(maxstations),fw60(maxstations),fw75(maxstations)
	REAL :: ts05(maxstations),ts10(maxstations),ts30(maxstations)
	REAL :: mp05(maxstations),mp25(maxstations),mp60(maxstations),mp75(maxstations)

	!Initialize variables
	stid='XXXX    '
	flag=0
	time=0.0
	date='::::::::::::'

	datafile='/dave/cgodfrey/mesonet/mesodata/d05/YYYYMMDDTTTT.d05'
	binfile='/dave/cgodfrey/mesonet/grads/YYYYMMDDTTTT.dat'
	call arguments(date)

	!Read soil temperature and moisture data for each station
	if(date(11:12).eq."00".or.date(11:12).eq."15".or.date(11:12).eq."30".or.date(11:12).eq."45")then
	  CALL RDMESO(stidarr,wc05,wc25,wc60,wc75,ts05,ts10,ts30,fw05,fw25,fw60,fw75,mp05,mp25,mp60,mp75,date)
	  !Set funny obs to missing
	  do i=1,nstations
	    if(ts05(i).lt.-900) ts05(i)=-999.00
	    if(ts10(i).lt.-900) ts10(i)=-999.00
	    if(ts30(i).lt.-900) ts30(i)=-999.00
	  enddo
	else
	  do i=1,nstations
	    wc05(i)=-999.00
	    wc25(i)=-999.00
	    wc60(i)=-999.00
	    wc75(i)=-999.00
	    ts05(i)=-999.00
	    ts10(i)=-999.00
	    ts30(i)=-999.00
	    fw05(i)=-999.00
	    fw25(i)=-999.00
	    fw60(i)=-999.00
	    fw75(i)=-999.00
	    mp05(i)=-999.00
	    mp25(i)=-999.00
	    mp60(i)=-999.00
	    mp75(i)=-999.00
	  enddo
	endif

	datafile(37:48)=date
	binfile(30:41)=date
	open(unit=1,file=datafile,status='old')   !Data file (*.d05)
	open(unit=2,file='/dave/cgodfrey/mesonet/plotutil/geomeso.tbl',status='old') 	    !Lat/lon table
	open(unit=3,file=binfile,form='unformatted')!GrADS binary file
	!open(unit=4,file='test.out',status='unknown')

	!Read latitude/longitude pairs for each station
	read(2,"42(/)")		!Header
	do 10 i=1,nstations
	  read(2,100)station(i),lat(i),lon(i),elev(i)
	  !write(*,*)'test',station(i),lat(i),lon(i),elev(i)
10	continue
	close(2)

	!Read data file header
15	read(1,"///")

	do 20 m=1,1000
	  !We're only writing one time level; extra flags not necessary
	  read(1,200,end=99)rstid,relh,tair,wspd,wvec,wdir,wdsd,&
	    wssd,wmax,rain,pres,srad,ta9m,ws2m
	  write(*,*)rstid,relh,tair,wspd,wvec,wdir,wdsd,&
	    wssd,wmax,rain,pres,srad,ta9m,ws2m

	  !Convert station ID to 8 characters
	  stid(1:4)=rstid

	  !Convert missing/non-QA data to -999
	  if(relh.lt.-99) relh=-999
	  if(tair.lt.-99) tair=-999.0
	  if(wspd.lt.-99) wspd=-999.0
	  if(wvec.lt.-99) wvec=-999.0
	  if(wdir.lt.-99) wdir=-999
	  if(wdsd.lt.-99) wdsd=-999.0
	  if(wssd.lt.-99) wssd=-999.0
	  if(wmax.lt.-99) wmax=-999.0
	  if(rain.lt.-99) rain=-999.0
	  if(pres.lt.-99) pres=-999.0
	  if(srad.lt.-99) srad=-999
	  if(ta9m.lt.-99) ta9m=-999.0
	  if(ws2m.lt.-99) ws2m=-999.0

	  !Search for appropriate lat/lon coords for station ID
	  match=0
	  index=0
	  do while (match.eq.0)
	    index=index+1
	    if(station(index).eq.stid(1:4)) match=1
	    if(index.gt.nstations)then
	      write(*,*)"ERROR: STATION DOES NOT EXIST ",stid
	    endif
	  enddo

	  !Search for correct soil variables to match station ID
	  match=0
	  soil=0
	  if(date(11:12).eq."00".or.date(11:12).eq."15".or.date(11:12).eq."30".or.date(11:12).eq."45")then
	    do while (match.eq.0)
	      soil=soil+1
	      if(stidarr(soil).eq.stid(1:4)) match=1
	      if(soil.gt.maxstations)then
	        write(*,*)"ERROR: SOIL STATION DOES NOT EXIST ",stid
	      endif
	    enddo
	  else
	    soil=1
	  endif

	  !Error check
	  if(stid(1:4).ne.station(index))then
	    write(*,*)"ERROR: STATION MISMATCH ",stid,station(index)," ",index,m
	    call abort()
	  endif
	  if(date(11:12).eq."00".or.date(11:12).eq."15".or.date(11:12).eq."30".or.date(11:12).eq."45")then
	    if(stid(1:4).ne.stidarr(soil))then
	      write(*,*)"ERROR: SOIL STATION MISMATCH ",stid,stidarr(soil)," ",soil,m
	      call abort()
	    endif
	  endif

	  !Write data
	  time=0.0
          nlev=1
          noflag=1
          write(3) stid,real(lat(index)),real(lon(index)),real(time),nlev,noflag
          write(3) real(relh),tair,wspd,wvec,real(wdir),wdsd,&
	    wssd,wmax,rain,pres,real(srad),ta9m,ws2m,wc05(soil),wc25(soil),&
	    wc60(soil),wc75(soil),ts05(soil),ts10(soil),ts30(soil),fw05(soil),&
	    fw25(soil),fw60(soil),fw75(soil),mp05(soil),mp25(soil),mp60(soil),mp75(soil)
          !write (4,900) stid,lat(index),lon(index),time,nlev,noflag
!900	  format(a8,2x,f7.4,2x,f9.4,2x,f3.1,2x,i1,2x,i1)
          !write (4,*)relh,tair,wspd,wvec,wdir,wdsd,&
	    !wssd,wmax,rain,pres,srad,ta9m,ws2m
	  !write(4,*) real(relh)

20	continue

	!At end of file write last time group terminator. 
99	nlev=0
	noflag=1
	write(3) stid,lat(index),lon(index),real(time),nlev,noflag
	!write(4,*)stid,lat(index),lon(index),time,nlev,noflag

	!NOTE: The GrADS control (descriptor) file is created by
	!mesoplot.cgi on webtest when a CGI plot is submitted.

100	format(5x,a4,55x,f7.4,1x,f9.4,1x,i4)
200	format(1x,a4,17x,i4,3(1x,f6.1),2x,i4,3(1x,f6.1),1x,f7.2,&
	  2x,f7.2,2x,i4,2(1x,f6.1))

	STOP
	END

	SUBROUTINE RDMESO(stidarr,wc05,wc25,wc60,wc75,ts05,ts10,ts30,fw05,fw25,fw60,fw75,mp05,mp25,mp60,mp75,datetime)
!	Modified from the 8/25/2005 version of RDMESO in SOIL/mesointerp.f
!
!	************************************************************
!	MESONET CONSTANTS (see geomeso.coef/geomeso.csv for details)
!	stnm	Station number
!	stid	Station identifier
!	name	Station name
!	city	Nearest incorporated town
!	rang	Range from town to station
!	cdir	Compass direction from town to station
!	cnty	County
!	nlat	North latitude
!	elon	East longitude
!	elev	Elevation in meters
!	qloc	Quality of geographic location
!	cdiv	Oklahoma climate division
!	clas	Station class (STANDARD/OASIS/SUPER)
!	subc	Station subclass
!	wcr05	5 cm residual water content (m^3 water/m^3 soil)
!	wcs05	5 cm saturated water content (m^3 water/m^3 soil)
!	a05	5 cm alpha constant (empirical)
!	n05	5 cm n constant (empirical)
!	wcr25	25 cm residual water content (m^3 water/m^3 soil)
!	wcs25	25 cm saturated water content (m^3 water/m^3 soil)
!	a25	25 cm alpha constant (empirical)
!	n25	25 cm n constant (empirical)
!	wcr60	60 cm residual water content (m^3 water/m^3 soil)
!	wcs60	60 cm saturated water content (m^3 water/m^3 soil)
!	a60	60 cm alpha constant (empirical)
!	n60	60 cm n constant (empirical)
!	wcr75	75 cm residual water content (m^3 water/m^3 soil)
!	wcs75	75 cm saturated water content (m^3 water/m^3 soil)
!	a75	75 cm alpha constant (empirical)
!	n75	75 cm n constant (empirical)
!	datc	Date commissioned
!	datd	Date de-commissioned
!
!	MESONET DATA
!	tr05	DeltaT reference temperature at a depth of 5 cm
!	tr25	DeltaT reference temperature at a depth of 25 cm
!	tr60	DeltaT reference temperature at a depth of 60 cm
!	tr75	DeltaT reference temperature at a depth of 75 cm
!	ts05	Sod soil temperature at a depth of 5 cm (degrees C)
!	ts10	Sod soil temperature at a depth of 10 cm (degrees C)
!	ts30	Sod soil temperature at a depth of 30 cm (degrees C)
!	error	Bias added to all DeltaT reference temperatures
!
!	DERIVED VARIABLES
!	mp05	Soil water potential/matric potential (bars) at 5 cm
!	mp25	Soil water potential/matric potential (bars) at 25 cm
!	mp60	Soil water potential/matric potential (bars) at 60 cm
!	mp75	Soil water potential/matric potential (bars) at 75 cm
!	wc05	Volumetric soil water content (m^3/m^3) at 5 cm
!	wc25	Volumetric soil water content (m^3/m^3) at 25 cm
!	wc60	Volumetric soil water content (m^3/m^3) at 60 cm
!	wc75	Volumetric soil water content (m^3/m^3) at 75 cm
!	fw05	Fractional water content at 5 cm
!	fw25	Fractional water content at 25 cm
!	fw60	Fractional water content at 60 cm
!	fw75	Fractional water content at 75 cm
!
!	REFERENCES:
!	Oklahoma Climatological Survey, 2002: Estimates of soil moisture
!		from the Oklahoma Mesonet.  Version 2.1, October 2002.
!		(How do I format this one?)
!	************************************************************

	IMPLICIT NONE
	DOUBLE PRECISION :: dist
	INTEGER :: geomesounit,moistunit,tempunit,maxstations,i,n
	INTEGER :: stnm,elev,qloc,cdiv,subc,datc,datd
	INTEGER :: nstations,mismatch,noread,override,test,skip
	PARAMETER (geomesounit=1)	!Coordinates and constants file
	PARAMETER (moistunit=2)		!Mesonet deltaT reference data file
	PARAMETER (tempunit=3)		!Mesonet soil temperature data file
	PARAMETER (maxstations=500)
	REAL :: rang,nlat,elon
	REAL :: dwcr05,dwcs05,da05,dn05,dwcr25,dwcs25
	REAL :: da25,dn25,dwcr60,dwcs60,da60,dn60,dwcr75,dwcs75,da75,dn75
	REAL :: wcr05(maxstations),wcs05(maxstations),a05(maxstations)
	REAL :: n05(maxstations),wcr25(maxstations),wcs25(maxstations)
	REAL :: a25(maxstations),n25(maxstations),wcr60(maxstations)
	REAL :: wcs60(maxstations),a60(maxstations),n60(maxstations)
	REAL :: wcr75(maxstations),wcs75(maxstations),a75(maxstations),n75(maxstations)
	DOUBLE PRECISION :: lat(maxstations),lon(maxstations)
	REAL :: tr05(maxstations),tr25(maxstations),tr60(maxstations),tr75(maxstations)
	REAL :: a,c
	REAL :: mp05(maxstations),mp25(maxstations),mp60(maxstations),mp75(maxstations)
	REAL :: wc05(maxstations),wc25(maxstations),wc60(maxstations),wc75(maxstations)
	REAL :: fw05(maxstations),fw25(maxstations),fw60(maxstations),fw75(maxstations)
	REAL :: ts05(maxstations),ts10(maxstations),ts30(maxstations)
	REAL :: countwc05,meanwc05,sumstdev,stdevwc05
	REAL :: countwc25,meanwc25,stdevwc25
	REAL :: countwc60,meanwc60,stdevwc60
	REAL :: countwc75,meanwc75,stdevwc75
	REAL :: ts05count,ts05mean,ts05stdev
	REAL :: ts10count,ts10mean,ts10stdev
	REAL :: ts30count,ts30mean,ts30stdev
	REAL :: error
	CHARACTER(len=4) :: stid,tstid,stidarr(maxstations),stid2(maxstations),time
	CHARACTER(len=13) :: name,city
	CHARACTER(len=12) :: cnty,mesodate,datecheck
	CHARACTER(len=3) :: cdir
	CHARACTER(len=8) :: clas,date
	CHARACTER(len=*) :: datetime
	CHARACTER(len=52) :: moistfile,tempfile
	CHARACTER(len=1) :: domain,rundir

	!Open Oklahoma Mesonet site description file
	open(unit=geomesounit,file='/dave/cgodfrey/mesonet/plotutil/geomeso.coef',&
	  status='old',action='read')

	!Read site locations and constants
	read(geomesounit,'74(/)')
	nstations=0
	do 10 i=1,1000
	  read(geomesounit,*,end=11)stnm,stid,name,city,rang,cdir,cnty,nlat,elon,elev,qloc,&
	    cdiv,clas,subc,dwcr05,dwcs05,da05,dn05,dwcr25,dwcs25,da25,dn25,dwcr60,dwcs60,&
	    da60,dn60,dwcr75,dwcs75,da75,dn75,datc,datd
	  !Remove stations that were de-commissioned before March 2004
	  if(datd.eq.20991231.or.datd.gt.20040300)then
	    nstations=nstations+1
	    lat(nstations)=nlat
	    lon(nstations)=elon	!yes, elon is negative already
	    stidarr(nstations)=stid
	    wcr05(nstations)=dwcr05
	    wcs05(nstations)=dwcs05
	    a05(nstations)=da05
	    n05(nstations)=dn05
	    wcr25(nstations)=dwcr25
	    wcs25(nstations)=dwcs25
	    a25(nstations)=da25
	    n25(nstations)=dn25
	    wcr60(nstations)=dwcr60
	    wcs60(nstations)=dwcs60
	    a60(nstations)=da60
	    n60(nstations)=dn60
	    wcr75(nstations)=dwcr75
	    wcs75(nstations)=dwcs75
	    a75(nstations)=da75
	    n75(nstations)=dn75
	  endif
10	continue
11	close(geomesounit)

	!*********
	!SOIL DATA
	!*********
	date=datetime(1:8)
	time=datetime(9:12)
	moistfile='/dave/cgodfrey/mesonet/mesodata/d30/YYYYMMDDTTTT.d30'
	tempfile='/dave/cgodfrey/mesonet/mesodata/d15/YYYYMMDDTTTT.d15'
	moistfile(37:44)=date
	moistfile(45:48)=time
	tempfile(37:44)=date
	tempfile(45:48)=time
	datecheck(1:8)=date
	datecheck(9:12)=time

	!To ease the programming burden, change the time for soil moisture if the requested time
	! is not on the half hour, read the data as usual, then set the data to missing at the
	! end of the subroutine. CMG 8/30/2005
	skip=0
	if(time(3:4).ne."00".and.time(3:4).ne."30")then
	  moistfile(45:48)="0000"
	  skip=1
	endif

	open(moistunit,file=moistfile,status='old',action='read')
	read(moistunit,'3(/)')	!Skip header
	open(tempunit,file=tempfile,status='old',action='read')
	read(tempunit,'3(/)')	!Skip header

	!Presumably, both files can be read at once, since for the same
	!date and time, the same stations should be present in both the
	!soil moisture file and the soil temperature file.  A check for
	!anything different is hard-coded below.

	!Read DeltaT reference temperatures and soil temperatures
	print*, 'READING OKLAHOMA MESONET DATA AT ',time,'Z FOR ',date
	mismatch=0
	noread=0
	do 20 i=1,1000
	  if(noread.eq.0)then

	    !SOIL MOISTURE
	    read(moistunit,100,end=21)stid,mesodate,tr05(i-mismatch),&
	      tr25(i-mismatch),tr60(i-mismatch),tr75(i-mismatch)
	    !Error check
	    if(mesodate.ne.datecheck.and.skip.eq.0)then
	      write(*,*)"***ERROR (81): SOIL MOISTURE FILE DATE DOES NOT MATCH FILE NAME"
	      STOP 81
	    endif

	    !SOIL TEMPERATURE (Sod temperature only)
21	    read(tempunit,200,end=22)tstid,mesodate,ts10(i-mismatch),&
	      ts05(i-mismatch),ts30(i-mismatch) !Yes, the order is correct
	    !Error check
	    if(mesodate.ne.datecheck)then
	      write(*,*)"***ERROR (82): SOIL TEMPERATURE FILE DATE DOES NOT MATCH FILE NAME"
	      STOP 82
	    endif
	    if(tstid.ne.stid)then
	      write(*,*)"***ERROR (83): MESONET STATION ID MISMATCH"
	      STOP 83
	    endif

	    !Adjust station ID array if stations have been newly decommissioned
	    stid2(i-mismatch)=stid
	    do 30 n=1,nstations
	      if(stidarr(n).eq.stid)then
	        !Adjust latitude/longitude coordinates to match station IDs
	        lat(i-mismatch)=lat(n)
	        lon(i-mismatch)=lon(n)
	        !Adjust instrumentation coefficients to match station IDs 
	        wcr05(i-mismatch)=wcr05(n)
	        wcs05(i-mismatch)=wcs05(n)
	        a05(i-mismatch)=a05(n)
	        n05(i-mismatch)=n05(n)
	        wcr25(i-mismatch)=wcr25(n)
	        wcs25(i-mismatch)=wcs25(n)
	        a25(i-mismatch)=a25(n)
	        n25(i-mismatch)=n25(n)
	        wcr60(i-mismatch)=wcr60(n)
	        wcs60(i-mismatch)=wcs60(n)
	        a60(i-mismatch)=a60(n)
	        n60(i-mismatch)=n60(n)
	        wcr75(i-mismatch)=wcr75(n)
	        wcs75(i-mismatch)=wcs75(n)
	        a75(i-mismatch)=a75(n)
	        n75(i-mismatch)=n75(n)
	      endif
30	    continue
	  endif

	  !Error check
	  if(i.gt.nstations.and.mismatch.eq.0)then
	    !New stations are always added to the end of the file (otherwise this won't work!)
	    write(*,*)"***ERROR (84): A STATION HAS BEEN ADDED!!! UPDATE GEOMESO.COEF."
	    STOP 84
	  endif
	  if(i.ge.nstations.and.mismatch.gt.0)then
	    write(*,*)"A STATION HAS BEEN DECOMMISSIONED!!! UPDATE GEOMESO.COEF."
	  endif
	  noread=0
	  if(stid.ne.stidarr(i))then
	    write(*,*)"STATION ID MISMATCH: ",stid," ",stidarr(i)
	    write(*,*)"   ATTEMPTING TO CONTINUE..."
	    mismatch=mismatch+1
	    noread=1
	  endif
20	continue
22	close(moistunit)
	close(tempunit)

	!Error check (modified 7/20/2005 by CMG)
	if(i.le.nstations.and.stidarr(nstations).eq."FITT")then
	  write(*,*)"NOTE: MESONET DATA FILE DOES NOT INCLUDE FITTSTOWN."
	  nstations=nstations-1
	elseif(i.le.nstations)then
	  write(*,*)"***ERROR (85): MESONET DATA FILE IS TOO SHORT!!! PLEASE CHECK IT."
	  STOP 85
	endif
	if(i.gt.500)then
	  write(*,*)"***ERROR (86): MAXSTATIONS EXCEEDED.  DATA/GEOMESO FILE MISMATCH"
	  STOP 86
	endif

	!Update the number of valid stations and station IDs if necessary
	nstations=nstations-mismatch
	stidarr=stid2

	a=1.788		!Constant coefficient for all sensors
	c=0.717		!Constant coefficient for all sensors
	do 40 i=1,nstations
	  !Introduce errors into soil moisture measurements for sensitivity studies
	  tr05(i)=tr05(i)+error
	  tr25(i)=tr25(i)+error
	  tr60(i)=tr60(i)+error
	  tr75(i)=tr75(i)+error

	  !Soil water potential (matric potential).  See OCS (2002).
	  if(tr05(i).gt.(-100.0))then
	    mp05(i)=-(c*exp(a*tr05(i)))/100.0
	  else
	    mp05(i)=-999.00
	  endif
	  if(tr25(i).gt.(-100.0))then
	    mp25(i)=-(c*exp(a*tr25(i)))/100.0
	  else
	    mp25(i)=-999.00
	  endif
	  if(tr60(i).gt.(-100.0))then
	    mp60(i)=-(c*exp(a*tr60(i)))/100.0
	  else
	    mp60(i)=-999.00
	  endif
	  if(tr75(i).gt.(-100.0))then
	    mp75(i)=-(c*exp(a*tr75(i)))/100.0
	  else
	    mp75(i)=-999.00
	  endif

	  !Remove suspect wet values only if 229-L measured funny values; ignore when adjusting for uncertainty
	  if(error.ge.-0.0001)then
	    !Values of MP between 0.00 and -0.10 bars are not accurate (OCS 2002)
	    if(mp05(i).gt.(-0.10)) mp05(i)=-999.00
	    if(mp25(i).gt.(-0.10)) mp25(i)=-999.00
	    if(mp60(i).gt.(-0.10)) mp60(i)=-999.00
	    if(mp75(i).gt.(-0.10)) mp75(i)=-999.00
	  endif

	  !Volumetric soil water content (m^3/m^3) (OCS 2002)
	  if(mp05(i).gt.(-100.0))then
	    wc05(i)=wcr05(i)+((wcs05(i)-wcr05(i))/((1.+((a05(i)*&
	      (-mp05(i)))**n05(i)))**(1.-(1./n05(i)))))
	  else
	    wc05(i)=-999.00
	  endif
	  if(mp25(i).gt.(-100.0))then
	    wc25(i)=wcr25(i)+((wcs25(i)-wcr25(i))/((1.+((a25(i)*&
	      (-mp25(i)))**n25(i)))**(1.-(1./n25(i)))))
	  else
	    wc25(i)=-999.00
	  endif
	  if(mp60(i).gt.(-100.0))then
	    wc60(i)=wcr60(i)+((wcs60(i)-wcr60(i))/((1.+((a60(i)*&
	      (-mp60(i)))**n60(i)))**(1.-(1./n60(i)))))
	  else
	    wc60(i)=-999.00
	  endif
	  if(mp75(i).gt.(-100.0))then
	    wc75(i)=wcr75(i)+((wcs75(i)-wcr75(i))/((1.+((a75(i)*&
	      (-mp75(i)))**n75(i)))**(1.-(1./n75(i)))))
	  else
	    wc75(i)=-999.00
	  endif

	  !CMG 8/25/2005
	  !NOTE: By default, if the coefficients in geomeso.coef are -999, the value
	  !of volumetric water content is also -999. Just to make sure everything goes
	  !smoothly, force wc## to be -999 if the coefficients are missing.
	  if(wcr05(i).lt.-900.0) wc05(i)=-999.00
	  if(wcr25(i).lt.-900.0) wc25(i)=-999.00
	  if(wcr60(i).lt.-900.0) wc60(i)=-999.00
	  if(wcr75(i).lt.-900.0) wc75(i)=-999.00

40	continue

	!Remove DRY outliers if QA has not done so already.  This was initially
	!necessary to remove Hobart volumetric water content data at 5 cm, but may
	!be applied to all other stations.  Volumetric water content information
	!is removed if the value is more than 3.5 sample standard deviations
	!below the mean of all stations.  Moist outliers are not affected.

	countwc05=0.0
	meanwc05=0.0
	countwc25=0.0
	meanwc25=0.0
	countwc60=0.0
	meanwc60=0.0
	countwc75=0.0
	meanwc75=0.0
	sumstdev=0.0
	!Calculate mean and standard deviation for wc05
	do i=1,nstations
	  if(wc05(i).ne.-999)then
	    meanwc05=meanwc05+wc05(i)
	    countwc05=countwc05+1.0
	  endif
	enddo
	meanwc05=meanwc05/countwc05
	do i=1,nstations
	  if(wc05(i).ne.-999)then
	    sumstdev=sumstdev+((wc05(i)-meanwc05)**2)
	  endif
	enddo
	stdevwc05=((1.0/(countwc05-1.0))*sumstdev)**(0.5)
	!Calculate mean and standard deviation for wc25
	sumstdev=0.0
	do i=1,nstations
	  if(wc25(i).ne.-999)then
	    meanwc25=meanwc25+wc25(i)
	    countwc25=countwc25+1.0
	  endif
	enddo
	meanwc25=meanwc25/countwc25
	do i=1,nstations
	  if(wc25(i).ne.-999)then
	    sumstdev=sumstdev+((wc25(i)-meanwc25)**2)
	  endif
	enddo
	stdevwc25=((1.0/(countwc25-1.0))*sumstdev)**(0.5)
	!Calculate mean and standard deviation for wc60
	sumstdev=0.0
	do i=1,nstations
	  if(wc60(i).ne.-999)then
	    meanwc60=meanwc60+wc60(i)
	    countwc60=countwc60+1.0
	  endif
	enddo
	meanwc60=meanwc60/countwc60
	do i=1,nstations
	  if(wc60(i).ne.-999)then
	    sumstdev=sumstdev+((wc60(i)-meanwc60)**2)
	  endif
	enddo
	stdevwc60=((1.0/(countwc60-1.0))*sumstdev)**(0.5)
	!Calculate mean and standard deviation for wc75
	sumstdev=0.0
	do i=1,nstations
	  if(wc75(i).ne.-999)then
	    meanwc75=meanwc75+wc75(i)
	    countwc75=countwc75+1.0
	  endif
	enddo
	meanwc75=meanwc75/countwc75
	do i=1,nstations
	  if(wc75(i).ne.-999)then
	    sumstdev=sumstdev+((wc75(i)-meanwc75)**2)
	  endif
	enddo
	stdevwc75=((1.0/(countwc75-1.0))*sumstdev)**(0.5)
	!Remove outlying DRY stations
	do i=1,nstations
	  if((meanwc05-(3.0*stdevwc05).gt.wc05(i)).and.(wc05(i).ne.-999))then
	    write(*,'(7x,"Dry station ",a4," removed from analysis of wc05: "&
	    &,f10.5,".")')stidarr(i),wc05(i)
	    wc05(i)=-999.00
	  endif
	  if((meanwc25-(3.0*stdevwc25).gt.wc25(i)).and.(wc25(i).ne.-999))then
	    write(*,'(7x,"Dry station ",a4," removed from analysis of wc25: "&
	    &,f10.5,".")')stidarr(i),wc25(i)
	    wc25(i)=-999.00
	  endif
	  if((meanwc60-(3.0*stdevwc60).gt.wc60(i)).and.(wc60(i).ne.-999))then
	    write(*,'(7x,"Dry station ",a4," removed from analysis of wc60: "&
	    &,f10.5,".")')stidarr(i),wc60(i)
	    wc60(i)=-999.00
	  endif
	  if((meanwc75-(3.0*stdevwc75).gt.wc75(i)).and.(wc75(i).ne.-999))then
	    write(*,'(7x,"Dry station ",a4," removed from analysis of wc75: "&
	    &,f10.5,".")')stidarr(i),wc75(i)
	    wc75(i)=-999.00
	  endif
	enddo

100	format(1x,a4,2x,a12,4(1x,f7.2))
200	format(1x,a4,2x,a12,1x,f6.1,2(8x,f6.1))

	!Fractional Water Index
	do i=1,nstations
	  if(tr05(i).gt.-900)then
	    fw05(i)=(3.96-tr05(i))/(3.96-1.38)
	  else
	    fw05(i)=-999.00
	  endif
	  if(tr25(i).gt.-900)then
	    fw25(i)=(3.96-tr25(i))/(3.96-1.38)
	  else
	    fw25(i)=-999.00
	  endif
	  if(tr60(i).gt.-900)then
	    fw60(i)=(3.96-tr60(i))/(3.96-1.38)
	  else
	    fw60(i)=-999.00
	  endif
	  if(tr75(i).gt.-900)then
	    fw75(i)=(3.96-tr75(i))/(3.96-1.38)
	  else
	    fw75(i)=-999.00
	  endif
	enddo

	close(geomesounit)
	close(moistunit)
	close(tempunit)

	!Set soil moisture to missing if time is not on the half hour
	if(skip.eq.1)then
	  do i=1,nstations
	    wc05(i)=-999.00
	    wc25(i)=-999.00
	    wc60(i)=-999.00
	    wc75(i)=-999.00
	    fw05(i)=-999.00
	    fw25(i)=-999.00
	    fw60(i)=-999.00
	    fw75(i)=-999.00
	  enddo
	endif

	RETURN
	END

	SUBROUTINE ARGUMENTS(date)
	implicit none
	character(len=*) :: date
	character(len=120) :: harg
	integer :: ierr,i,numarg
	integer, external :: iargc

	numarg=iargc()

	i=1
	call getarg(1,harg)
	date=harg
	return
	end