C Copyright(c) 1997, Space Science and Engineering Center, UW-Madison
C Refer to "McIDAS Software Acquisition and Distribution Policies"
C in the file  mcidas/data/license.txt

C *** $Id: mcdewpt.f,v 1.8 1998/11/30 15:12:04 billl Rel $ ***

*$ Name:
*$	mcdewpt - Computes dewpoint temperature given t and rh.
*$
*$ Interface:
*$	integer function
*$	mcdewpt(double precision t, double precision rh,
*$         double precision td)
*$
*$ Input:
*$	t  - temperature (Kelvin)
*$	rh - relative humidity
*$
*$ Input and Output:
*$	none
*$
*$ Output:
*$	td - dewpoint temperature (Kelvin)
*$
*$ Return values:
*$        1     - computed td is greater than temperature
*$        0   	- all is okay
*$       -3     - input t is not reasonable.  See remarks
*$       -4     - input rh is not reasonable.  See remarks
*$
*$ Remarks:
*$	This subroutine computes dewpoint temperature using the Clausius-
*$        Clapeyron equation:  RH = exp((L/Rv)*(1/t - 1/td)) where L is the
*$        latent heat of vaporization and Rv is the gas constant for water
*$        vapor:  L = 2500800. - (t - 273.15)*2397.5 and Rv = 461.5.  To
*$        compute dewpoint, then:
*$
*$        Td = 1./[ 1/T - (Rv ln (RH))/Lv ]
*$
*$      If the RH is reported as less than 1%, it is set to 1%.  Since the
*$        TD computation involved the log of RH, as it goes to zero, the TD
*$        goes to zero.  This is not observed too often.
*$
*$      This program requires that input values are reasonable.  RH cannot
*$        exceed 110 (110%) or be less than 0, T must be between 150 and 400K.
*$
*$ Categories:
*$	met/science

        integer function
     &	  mcdewpt(t,rh,td)

	implicit none
	include 'hex80.inc'

C-----symbolic constants and shared data

	double precision	t	! temperature, assumed Kelvin
	double precision	rh	! relative humidity, 0<=RH<=1.1
	double precision	td	! dewpoint temperature

C-----external functions

C-----Local variables....

 	integer		 rc		! return code
	double precision gasconst	! moist gas constant
	double precision latheat	! latent heat of vaporization
	double precision relhum 	! relative humidity

        data gasconst/4.615D2/
	data latheat/2500800.D0/

	rc = 0

	if (t .lt. 150.D0 .or. t .gt. 400.D0) then
	   rc = -3
	   goto 2000
	endif

	if ( rh .gt. 110.D0 .or. rh .lt. 0.D0) then
	   rc = -4
	   goto 2000
	endif

	relhum = rh/100.D0

C--- constrain the lower bound of RH to be 1% -- otherwise the
C	log of RH gets a little problematic...

	if (relhum .lt. 0.01D0)  relhum = 0.01D0

C---now compute td

	td = 1.D0/(1.D0/t - (gasconst*dlog(relhum))/(latheat -
     &			(t-273.15D0)*2397.5D0) )

	if (td .gt. t) then
	   rc = 1
	endif

 2000   mcdewpt = rc
	return
        end

*$ Name:
*$	mcdewptgrd - Compute dewpoint temperature given t, rh, and p.
*$
*$ Interface:
*$	integer function
*$	mcdewptgrd(integer gridt(*) integer gridrh(*)
*$         integer gridtd(*) headt(*) integer headrh(*)
*$	   integer headtd(*) )
*$
*$ Input:
*$	headt  - gridheader of temperature
*$	headrh - gridheader of relative humidity
*$	gridt  - grid of temperature
*$	gridrh - grid of relative humidity
*$
*$ Input and Output:
*$	none
*$
*$ Output:
*$	headtd - gridheader of dewpoint
*$	gridtd - dewpoint grid
*$
*$ Return values:
*$	 0   	- all is okay
*$	 -1	- t and rh grids are of different sizes/times
*$	 -2	- t and rh grids have different projections
*$       -3     - RH computation failed for all points
*$
*$ Remarks:
*$	Computes a grid of dewpoint temperature.  This function sets up the
*$        grid header.  It calls mcdewpt to actually compute dewpoint.
*$      If the RH computation fails for a point, then the RH is set to missing.
*$        If the computation fails for all points, then an error status is
*$        returned.  The computation will fail if the RH is everywhere >110%
*$        or < 0%.
*$
*$ Categories:
*$	grid

        integer function
     &	  mcdewptgrd(gridt,gridrh,gridtd,headt,headrh,headtd)

	implicit none
	include 'gridparm.inc'
	include 'hex80.inc'

C-----symbolic constants and shared data

	integer gridt(*)		! grid data of u-component velocity
	integer gridrh(*)		! grid data of v-component velocity
	integer gridtd(*)		! grid data of vorticity
	integer headt(*)		! header for u-component velocity
	integer headrh(*)		! header for v-component velocity
	integer headtd(*)		! header for vorticity grid

C-----external functions

	integer lit
	integer mcdewpt

C-----Local variables....

        character*80     errstring
        integer		 bad		! RH at this point is bad
        integer		 badpoint	! number of points with bad rh
        integer		 i		! loop bound
 	integer		 iret		! return code
        integer		 miss		! missing value code
 	integer		 rc		! return code
	integer		 rhscale	! scaling for relative humidity grid
	integer		 tscale		! scaling for temperature grid
        double precision temp		! temperature
	double precision rh		! relative humidity
	double precision td		! dewpoint

	data miss/HEX80/

        badpoint = 0
        bad = 0
	rc = 0

C-----Check grids to see if they are the same time/space
C    Check gridheaders 1-6:  total size, nrows, ncols

        do 100 i = 1,6

           if (   headt(i) .ne. headrh(i)  )  then

		rc = -1
		goto 2000
		endif

 100       continue


C----Check gridheaders 34-40:  Projection parameters

        do 110 i = 34,40

           if (   headt(i) .ne. headrh(i)  )  then

		rc = -2
		goto 2000
		endif

 110       continue

C----Copy gridheader into dewpt gridheader

         do 120 i = 1,64

           headtd(i) = headt(i)
	   if (i .ge. 53)  headtd(i) = 0

 120       continue

C-----change name in header

	headtd(9) = lit('K   ')
	headtd(7) = lit('TD  ')
	call movcw('Dewpoint temperature',headtd(53))
	tscale = 10.**headt(8)
        rhscale = 10.**headrh(8)

        do 200 i = 1,headt(1)

C---check if Temp or RH is missing

	     if (gridt(i) .eq. miss  .or.
     &			gridrh(i) .eq. miss) then
		gridtd(i) = miss

	     else

C---check to make sure the temperature is in kelvin;  convert
C	if not

	        temp = gridt(i)/tscale

                if (headt(9) .eq. lit('C   ') )  then
		    temp = temp + 273.15D0
		elseif (headt(9) .eq. lit('F   ') ) then
		    temp = (5.D0/9.D0)*(temp - 32.D0) + 273.15D0
		endif

		rh = gridrh(i)/rhscale

C -- mcdewpt assumes that the relative humidity is stored as a number between
C	0 and 110 -- 100% RH is saturated, but we'll allow for a little
C	supersaturation
C
		if (rh.gt.110.D0 .or. rh .lt. 0.D0) then
		   badpoint = badpoint + 1
		   gridtd(i) = miss
		else

C---now compute td -- but no need to if the RH is bad

		   iret = mcdewpt(temp,rh,td)
		   if (iret .lt. 0)  then

		      badpoint = badpoint + 1
		      gridtd(i) = miss

		   else

C --- if the computed dewpoint if much greater than
C	the temperature, return a missing value and
C	increment the number of bad points

		      if (td .gt. (temp+0.05D0)) then
		         badpoint = badpoint + 1
			 gridtd(i) = miss

C --- if the computed dewpoint is slightly greater than
C	the temperature, return the temperature -- this
C	is allowing for roundoff error.  Note that it
C	also reduces supersaturated conditions to
C	simply saturated.

		      elseif (td .gt. temp) then
			 gridtd(i) = nint(temp*tscale)
		      else

C --- Finally, we reach the 'normal' case...the dewpoint
C	is computed, and it's less than the temperature!

			 gridtd(i) = nint(td*tscale)
		      endif

		   endif

	        endif

	     endif

 200       continue

	if (badpoint .eq. headt(1)) then
	   rc = -3
	   goto 2000
        endif

 2000   continue

	mcdewptgrd = rc
	return
        end
