program sample_application

c     This is a sample code demonstrating application of the diurnal
c     warming codes.  The main work is done in the subroutine applyLUT,
c     but this code illustrates how that code would be called with a
c     sample of cruise data.  

c     This test is for data at a single point (or moving track).  The
c     code can be generalized to application on a grid, but separate
c     storage of the integrated wind speed and insolation must be 
c     maintained at each point on the grid.

c     Created by Sandra Castro and Gary Wick
c     University of Colorado/CCAR and NOAA ESRL
c     March 2009

c     External routines
c        applyLUT

	implicit none

	real jdayutc, jdaylocal, utchour, lat, lon
	real cws, csol, ngwd, ngsol
	real lasthr, hr, dth
	real qsolint, wsint
	real dtinst, dtintq, dtintuq

	integer iyr, doy
	integer swap
	integer i

c     The LUT are provided for both big and little endian architecture
c     and the application code is capable of handling both (see comments
c     in applyLUT.f)

c     Initialize swap for linux/PC architecture
	swap = 0

c     The code first opens the sample cruise data

	open(21,file='sample_cruisedata.dat',form='formatted')

c     Fields in the sample data are as follows:
c     1: year
c     2: fractional julian day in utc time (eg. 50.38)
c     3: fractional julian day in local time
c     4: hour in utc time
c     5: latitude
c     6: longitude
c     7: wind speed from ship data (m/s)
c     8: insolation from ship data (W/m2)
c     9: corresponding wind speed from NOGAPS data
c    10: corresponding insolation from NOGAPS data

c     Note that the data are provided every 3 hours as might be expected
c     from NOGAPS data.  A key factor is that the data begin at 6 am local
c     time.  This is important for initiation of the integrated wind speed
c     and insolation values, where the tables are designed for integration
c     from this time.  Application of the code should only begin with data
c     at 6 am.  For other start times, the data should either be skipped
c     until the following 6 am, or the integrated quantities from the
c     previous 6 am must be estimated.

c     We initialize the hour corresponding to the time of what the last 
c     data prior to the available data would have been.  This is part of
c     setting up proper integration of the wind speed and insolation from
c     6 am local time.

c     Since the data are 3 hourly and the data begin at 6 am, we assume
c     the last data would have been at 3 am.

	lasthr = 3.0

c     The data have 23 lines so we enter an explicit loop through the
c     23 time steps

	do 5 i=1,23

c     Read in the data line

	   read(21,*) iyr, jdayutc, jdaylocal, utchour, lat, lon, cws,
     A                csol, ngwd, ngsol

c     Extract the local day of year and hour from jdaylocal

	   doy = int(jdaylocal)
	   hr = (jdaylocal - doy) * 24.0

c    We now compare the current hour with the previous hour to see if we
c    have passed 6 am local time.  If so, the integrated wind speed and
c    insolation must be reinitialized to 0

	   if (hr .ge. 6.0 .and. lasthr .lt. 6.0) then

c    Initialize integrated wind speed and insolation to zero
	      qsolint = 0.0
	      wsint = 0.0

c     We also reset the last hour to 6 am so that integration can take
c     place from 6 am

	      lasthr = 6.0

	   endif

c     Prior to calling the code to apply the LUT, compute the time step
c     (in hours) from the last data point for use in integrating the wind
c     speed and insolation.

	   dth = hr - lasthr

c     Must watch for day change where hr becomes < lasthr

	   if (dth .lt. 0.0) then
	      dth = dth + 24.0
	   endif

c     Now all quantities are available to call the applyLUT routine.  Note
c     that integration of the wind speed and insolation are done in the
c     routine and the updated integrated quantities are returned back to this
c     driver code

	   call applyLUT(jdayutc, doy, hr, cws, csol, dth, lat, lon,
     A                   swap, wsint, qsolint, dtinst, dtintq, 
     B                   dtintuq)

c     Print out the extracted warming values
c     dtinst:  Warming from instantaneous table
c     dtintq:  Warming from instantaneous wind speed and integrated solar
c              table
c     dtintuq: Warming from integrated wind speed and insolation LUT

c	   print*, cws, wsint, csol, qsolint
	   print*, hr, dtinst, dtintq, dtintuq

c     Before moving to the next time step, replace lasthr with the current
c     hour

	   lasthr = hr

   5	continue

	close(21)

	stop
	end