subroutine filter(tb19v, tb19h, tb22v, tb37v,          &
                     tb37h, tb85v, tb85h, ipass, iprecip, xlat )
   implicit none

   real    ,  intent (in)  :: tb19v, tb19h, tb22v, tb37v, &
                              tb37h, tb85v, tb85h, xlat
   real                     :: tb(7)
   logical ,  intent (inout):: ipass
   integer ,  intent (inout):: iprecip

   integer                 :: ipass1, ipass2
   real                    :: sst, dsst
   real                    :: theta,tc,wv,u,alw19,alw37,alw85

   !  write (unit=stdout,fmt=*) 'tbdata',tb19v, tb19h, tb22v, tb37v,tb37h, tb85v, tb85h
   ! tc : cloud top temperature (C)

   theta = 53.1
   tc = 10.0
   sst = 23.0
   dsst = 1000.0

   ipass1 = 0
   ipass2 = 0

   tb(1) = tb19v 
   tb(2) = tb19h
   tb(3) = tb22v
   tb(4) = tb37v
   tb(5) = tb37h
   tb(6) = tb85v
   tb(7) = tb85h

   ! upper and lower boundary 19h (90,280), others (100,280)

   if (tb(1) .lt. 280 .and. tb(1) .gt. 100. .and. &
       tb(2) .lt. 280 .and. tb(2) .gt.  90. .and. &
       tb(3) .lt. 280 .and. tb(3) .gt. 100. .and. &
       tb(4) .lt. 280 .and. tb(4) .gt. 100. .and. &
       tb(5) .lt. 280 .and. tb(5) .gt. 100. .and. &
       tb(6) .lt. 280 .and. tb(6) .gt. 100. .and. &
       tb(7) .lt. 280 .and. tb(7) .gt. 100. .and. &

       !    horizontal always much less than vertical

       tb(1) .gt. tb(2) .and. &
       tb(4) .gt. tb(5) .and. &
       tb(6) .gt. tb(7) .and. &

       ! T19V always at least 4 K less than  T22V, except in heavy rain

       tb(3)-tb(1) .gt. 4.) then
      ipass1 = 1
   end if

   ! second check

   if (ipass1 .eq. 1) then

      call pettyalgos(tb,theta,tc,sst,dsst,          &
                          wv,u,alw19,alw37,alw85,iprecip, xlat)

         if (wv    .gt. -999. .and. u     .gt. -999. .and. &
             alw19 .gt. -999. .and. alw37 .gt. -999. .and. &
             alw85 .gt. -999. ) then
            ipass2 = 1
            ! write(unit=stdout,*)a1,a2,lat,long,theta,tb(1),tb(2),tb(3), &
            !    tb(4),tb(5),tb(6),tb(7)
         end if

      ! third check

      ! if (ipass2 .eq. 1) then
      !    if (wv    .ge. 5.    .and. wv .le. 70 .and.     &
      !        u     .ge. -5    .and. u  .le. 30 .and.     &
      !        alw19 .ge. -0.01 .and. alw19 .le. 0.5 .and. &
      !        alw37 .ge. -0.01 .and. alw37 .le. 0.5 .and. &
      !        alw85 .ge. -0.01 .and. alw85 .le. 0.5) then
      !       write(unit=stdout,*)a1,a2,lat,long,theta,tb(1),tb(2),tb(3), &
      !          tb(4),tb(5),tb(6),tb(7)
      !       ipass = .true.
      !    end if
      ! end if
      ! if (iprecip .eq. 1) then
      !    ipass = .true.
      !    write(unit=stdout,*) 'precipitation',lat,long,wv,u,alw19,alw37,alw85,iprecip
      ! else
      !    ipass = .true.
      ! end if
   end if

   ! write(unit=stdout,*) 'pass=',ipass1,ipass2,iprecip,ipass

end subroutine filter

   !***************************************************************
   !                    NOTICE
   ! 
   ! This package of subroutines has been developed by Grant W. Petty
   ! (Purdue University) for the retrieval of column water vapor, column
   ! cloud liquid water, and surface wind speed from SSM/I brightness
   ! temperatures.  The algorithm versions contained herein were developed
   ! in part using data sets provided by NASDA/EORC, Japan and are
   ! submitted to NASDA for evaluation as part of ADEOS-2/AMSR algorithm
   ! development activities.  Because of the need to meet a submission
   ! deadline, they have been subject only to limited testing, and it is
   ! possible that undetected bugs or problems remain in the code that
   ! follows.
   ! 
   ! Problems or questions should be directed to 
   !       Grant W. Petty
   !       Earth and Atmospheric Sciences Dept.
   !       Purdue University
   !       West Lafayette, in 47907-1397
   !       (765) 494-2544
   !       gpetty@rain.atms.purdue.edu
   ! 
   ! Version:   July 13, 1997
   ! Copyright 1997, Grant W. Petty 
   !
   ! Permission is granted to all scientific users to utilize
   ! and/or test this set of algorithms without restriction, provided only
   ! that
   !       (1) they are not used for commercial Purposes without the
   !          author's permission
   ! 
   !       (2) the author is given due credit for their development
   ! 
   !       (3) any modifications to the algorithms by 3d parties must be
   !          clearly identified as such and distinguished from the author's
   !          original code
   ! 
   !********************************************************************
   !      GENERAL COMMENTS ON USAGE
   ! 
   ! The algorithms that follow can act on single-pixel values of brightness
   ! temperature and derived variables, such as water vapor, wind speed,
   ! etc.  However, the way they are normally used at Purdue University is
   ! as follows:
   ! 
   !  1    Retrieve Water Vapor and Surface Wind Speed on a pixel-by-pixel basis
   ! 
   !  2    Apply mild spatial smoothing to derived wind speed and water vapor
   !       fields
   ! 
   !  3    Used filtered values of wind speed and water vapor to compute
   !       cloud-free polarization differences
   ! 
   !  4    Compute liquid water from unsmoothed, full-resolution brightness
   !       temperatures, normalized by smoothed cloud-free polarizations.
   ! 
   ! The smoothing step described above is optional and should not have a
   ! strong effect on the results outside of precipitation.  Normally the
   ! smoothing step is necessary only for interpolating water vapor and
   ! wind speed into areas of precipitation, where normalized polarizations
   ! at 19, 37 and 85 GHz are needed for precipitation retrievals.
   ! 
   ! 
   !********************************************************************
   ! This subroutine accepts SSM/I  brightness temperatures and returns
   ! all over-ocean variables using the algorithms of Petty (1997, unpublished)
   ! It, and the algorithms it calls, are generally intended for use only outside
   ! of precipitation, except for the water vapor algorithm, which is believed to
   !  tolerate significant precipitation without serious errors.
   !
   ! Inputs:   tb(7) = 19v,t19h,t22v,t37v,t37h,t85v,t85h
   !                --- SSM/I brightness temperatures (K) 
   !           theta -- sensor viewing angle
   !           tc  --- assumed cloud temperature in degrees Celsius
   !           Note: If TC not known, set to a reasonable 
   !                 average value, such as 8.0 degrees
   !
   !           sst   -- sea surface temperature (deg. C)
   !           dsst   -- uncertainty in the above SST (deg. C)
   !                     (if dsst > 2.8 C, then an wind speed algorithm is
   !                      used that doesn't require SST)
   !
   ! Output:   wv     ---- column water vapor (kg/m**2) from 19.35 GHz
   !           u      ---- surface wind speed (m/sec)
   !           alw19  ---- column cloud liquid water (kg/m**2) from 19.35 GHz
   !           alw37  ---- column cloud liquid water (kg/m**2) from 37 GHz
   !           alw85  ---- column cloud liquid water (kg/m**2) from 85.5 GHz
   !           iprecip --- set to 1 if precipitation flag invoked, else 0
   ! 
   ! 
   ! NOTE: It is not expected that liquid water values from different SSM/I
   ! frequencies will always closely agree, on account of significant
   ! differences in spatial resolution, sensitivity, and dynamic range.
   ! For best comparisons, average results to a common spatial resolution
   ! first.  19 GHz channels have the least sensitivity to thin clouds; 85
   ! GHz channels are most severely affected by precipitation and/or strong
   ! inhomogeneity in cloud fields.
   ! 
   ! Also, liquid water values in excess of 0.5 kg/m**2 are generally
   ! flagged as precipitating, and a value of MISSinG is returned.  There
   ! is no known method for separating cloud water from precipitation
   ! water.  Furthermore, any attempt to retrieve total liquid water path
   ! when precipitation is present will require a considerably more
   ! sophisticated algorithm than the ones provided here.
   !******************************************************************

subroutine pettyalgos(tb,theta,tc,sst,dsst, &
   wv,u,alw19,alw37,alw85,iprecip,xlat)

   implicit none

   real tb(7),theta,tc,wv,u,alw19,alw37,alw85,xlat, sst,dsst
   real, parameter :: BAD = -1000.0
   integer iprecip

   ! get water vapor
   !     write (unit=stdout,fmt=*) 'tb=',tb
   call pettyvapor(tb,wv)
   !     write (unit=stdout,fmt=*) 'wv=',wv

   ! get surface wind speed
   call pettywind(tb,theta,sst,dsst,u)

   ! get column cloud liquid water, using 3 different algorithms
   call pettylwp(tb,theta,1,tc,alw19,iprecip,xlat)
   call pettylwp(tb,theta,2,tc,alw37,iprecip,xlat)
   call pettylwp(tb,theta,3,tc,alw85,iprecip,xlat)
   if (alw19 .eq. BAD .or. alw19 .eq. BAD .or. alw85 .eq. BAD) iprecip=1
   !     if (iprecip .eq. 1) write (unit=stdout,fmt=*) 'iprecip B=',iprecip,'alw=',alw19,alw37,alw85

end subroutine pettyalgos


   !********************************************************************
   ! This subroutine accepts SSM/I brightness temperatures and returns
   ! cloud liquid water in kg/m**2.  
   ! 
   ! Inputs:   tb(7) = 19v,t19h,t22v,t37v,t37h,t85v,t85h
   !                --- SSM/I brightness temperatures (K) 
   !           theta -- sensor viewing angle
   !           ifreq -- frequency to use for computing cloud liquid water
   !                    (1 = 19.35 GHz  2 = 37 GHz  3 = 85.5 GHz)
   !           tc  --- assumed cloud temperature in degrees Celsius
   !           Note: If not known, set to a reasonable value, such as 8.0 degrees
   !
   ! Output:   alw  ---- column water vapor (mm)
   !           iprecip --- set to 1 if precipitation flag invoked, else 0
 
subroutine pettylwp(tb,theta,ifreq,tc,alw,iprecip,xlat)

   implicit none
   real tb(7),theta,alw,tc,xlat
   integer ifreq,iprecip
   real, parameter :: BAD = -1000.0
   real t19v,t19h,t22v,t37v,t37h,t85v,t85h,MISSinG,pclr
   real :: p,s85,v0,wv,u
   integer :: ich
 
   iprecip = 0
   MISSinG = -8888.
   pclr=0.
   if (ifreq.lt.1 .or. ifreq.gt.3) then
      pause "BAD VALUE OF ifREQ in PettyLWP"
      alw = BAD
      return
   end if

   ! initialize brightness temperatures

   t19v = tb(1)
   t19h = tb(2)
   t22v = tb(3)
   t37v = tb(4)
   t37h = tb(5)
   t85v = tb(6)
   t85h = tb(7)

   ! estimate the clear-sky polarization difference for the scene

   call clearpol(tb,theta,ifreq,pclr,iprecip)

   ! if value of Tc is invalid, then set to reasonable intermediate value
   if (tc .lt. -20.0 .or. tc .gt. 40.0) then
      tc = 8.0
   end if

   ! calculate normalized polarization
   if (ifreq .eq. 1) then
      p = (t19v-t19h)/pclr
   else if (ifreq .eq. 2) then
      p = (t37v-t37h)/pclr
   else if (ifreq .eq. 3) then
      p = (t85v-t85h)/pclr

      ! at 85.5 GHz, there could be errors due to ice, so check S85
      call PettyVapor(tb,wv)
      if (wv .eq. BAD) then
         alw = BAD
         return
      end if
      call PettyWind(tb,theta,-99.9,1000.0,u)
      if (u .eq. BAD) then
         u = 5.0
         alw = BAD
      end if
     ich = 6 
     call tb_clear(ich,u,wv,v0)
     s85 = p*v0 + (1.0-p)*t_kelvin - t85v 
     if (s85 .gt. 10.0 .or. (t85v-t85h) .lt. 5.0) p = MISSinG 
   end if

   ! convert normalized polarization to LWP
   if (p .ne. MISSinG .and. p .lt. 1.4 .and. p .gt. 0.1) then
      call P2LWP(ifreq,p,tc,theta,alw,xlat)
   else
      alw = BAD
   end if

end subroutine PettyLWP

   !*********************************************************
   ! This routine computes liquid water path from normalized polarization
   ! at 19.35, 37, or 85.5 GHz  (indicated by ifREQ)

subroutine p2lwp(ifreq,p,tc,theta,alw,xlat) 

   implicit none        
   integer ifreq
   real p,alw,xlat,tc,theta
   real alpha(3),threshhold
   data alpha/2.08,2.05,1.78/
   real, parameter :: BAD = -1000.0
   real :: ext, costheta
   integer :: iprecip

   ! check for extreme value due to precipitation
   if (p .lt. 0.01) then
      alw = BAD
      iprecip = 1
      return
   end if

   !  get liquid water mass extinction coefficient
   ext = alw_ext(ifreq,tc)

   ! compute colum cloud liquid water using Beer's Law
   costheta = cos(theta*pi/180.0)
   alw = (-costheta/(alpha(ifreq)*ext))*log(p)

   ! check for precipitation
   threshhold=0.5
   !     write (unit=stdout,fmt=*) 'xlat=',xlat
   if (xlat .gt. 20) threshhold=0.3
   if (xlat .gt. 30) threshhold=0.1
   !     if (xlat .gt. 40) threshhold=0.25
   !     if (alw .gt. 0.5) then
   if (alw .gt. threshhold) then
      alw = BAD
      iprecip = 1
   end if

end subroutine p2lwp

   !***********************************************************************
   ! This subroutine accepts SSM/I brightness temperatures from a scene and returns
   ! the predicted cloud-free polarization difference for that scene
   ! 
   ! Inputs:   tb(7) = 19v,t19h,t22v,t37v,t37h,t85v,t85h
   !                --- SSM/I brightness temperatures (K) 
   !           theta -- sensor viewing angle
   !           ifreq -- frequency to use for computing cloud liquid water
   !                    (1 = 19.35 GHz  2 = 37 GHz  3 = 85.5 GHz)
   !
   ! Output:   pclr  ---- cloud free polarization difference (K)
   !           iprecip --- set to 1 if precipitation flag invoked, else 0

subroutine ClearPol(tb,theta,ifreq,pclr,iprecip)

   implicit none
   real tb(7),theta,pclr
   integer ifreq,iprecip
   real, parameter :: BAD = -1000.0
   real :: wv,u
   integer :: ich

   ! check for valid ifREQ
   if (ifreq.lt.1 .or. ifreq.gt.3) then
      pause 'BAD VALUE OF ifREQ in ClearPol'
      pclr = BAD
      return
   end if

   iprecip = 0

   ! get water vapor estimate.  If this value is returned as 'bad', then
   ! either the brightness temperatures are not valid for over ocean or
   ! else there is unusually heavy precipitation

   call PettyVapor(tb,wv)
   if (wv .lt. 0.0) then
      pclr = BAD
      return
   end if

   ! get wind speed estimate.  If this value is returned as 'bad', assume
   ! that it is due to precipitation, and substitute a reasonable 'average'
   ! wind speed.

   call PettyWind(tb,theta,-99.9,1000.0,u)
   if (u .eq. BAD) then
      iprecip = 1
      u = 5.0
   end if

   ! now that we finally have estimates of water vapor and wind speed,
   ! compute the clear-sky polarization difference at the desired frequency

   ich = ifreq + 7
   call tb_clear(ich,u,wv,pclr)

         end subroutine ClearPol

   !**********************************************************
   ! The following function returns the microwave mass extinction coefficient
   ! of cloud water at 19.35, 37.0, or 85.5 GHz.
   !
   !     inputs:
   !           ifreq -- frequency to use for computing cloud liquid water
   !                    (1 = 19.35 GHz  2 = 37 GHz  3 = 85.5 GHz)
   !           tc  --- cloud temperature in degrees Celsius
   !
   !      returned value:  Mass extinction coefficient (m**2/kg)
   !

function alw_ext(ifreq,tc)

   implicit none

   integer, intent(in) :: ifreq
   real, intent(in) :: tc

   real :: alw_ext
   real :: tc2,tc3

   tc2 = tc*tc
   tc3 = tc2*tc

   if (ifreq .eq. 1) then
      alw_ext = exp(-2.55-2.98e-2*tc-6.81e-4*tc2+3.35e-6*tc3)
   else if (ifreq .eq. 2) then
      alw_ext = exp(-1.35-0.0234*tc-1.22e-4*tc2+5.48e-6*tc3) 
   else if (ifreq .eq. 3) then
      alw_ext = exp(-0.0713-8.00e-3*tc-1.62e-4*tc2+1.18e-6*tc3)
   end if

end function alw_ext

   !*****************************************************************************
   ! The following routine returns approximate clear-sky SSM/I brightness temperature 
   ! over the ocean as a function of column water vapor and surface wind speed.
   ! This implementation is a bivariate polynomial fit to model computed brightness 
   ! temperatures (see Petty (1992,1994)), where the model was in turn carefully calibrated
   ! empirically against 200,000 cloud-free SSM/I pixels.
   ! This routine assumes a viewing angle not too different from 53.4 degrees, and
   ! it assumes "typical" values for other ocean and atmospheric parameters.
   ! 
   !      inputs:
   !         ich = channel no. (1 = 19V, 2 = 19H, ..., 7 = 85V, 
   !                            8 = 19V-19H, 9 = 37V-37H, 10 = 85V-85H)
   !         u = estimated surface wind speed (m/sec)
   !         wv = estimated column water vapor (kg/m**2)
   !
   !      output:
   !         tb = brightness temperature

subroutine tb_clear(ich,u,wv,tb)
   implicit none
   integer, intent(in)  :: ich
   real,    intent(in)  :: u
   real,    intent(in)  :: wv
   real,    intent(out) :: tb

   real :: wvp,up,x,x2,x3,y,y2,y3,xy,x2y,xy2

   real a(10,10)
   data a/ &
    0.1984E+03, 0.1674E+02, 0.1998E+00,-0.3868E+00,-0.1818E+00, &
    0.7065E-01,-0.8148E+00, 0.6116E-01, 0.1400E-01, 0.0000E+00, &
    0.1334E+03, 0.2587E+02, 0.3819E+01,-0.7762E-01, 0.1564E+00, &
    0.4968E-01,-0.1232E+01,-0.6779E-01,-0.3083E-01, 0.0000E+00, &
    0.2301E+03, 0.2865E+02, 0.5589E+00,-0.4668E+01,-0.7274E-01, &
    0.1043E-01,-0.6731E+00, 0.8141E-02,-0.1081E-01, 0.0000E+00, &
    0.2126E+03, 0.1103E+02,-0.2488E+00, 0.8614E+00,-0.1298E+00, &
    0.4531E-01,-0.8315E+00, 0.5104E-01, 0.1785E-01, 0.0000E+00, &
    0.1500E+03, 0.1864E+02, 0.3950E+01, 0.2095E+01,-0.1743E+00, &
    0.1170E+00,-0.1345E+01,-0.8584E-01, 0.8143E-02, 0.0000E+00, &
    0.2589E+03, 0.1706E+02,-0.7729E+00,-0.2905E+01, 0.2821E+00, &
    0.3375E-02,-0.4366E+00,-0.1150E+00, 0.1881E-01, 0.0000E+00, &
    0.2265E+03, 0.3363E+02, 0.2050E+01,-0.5701E+01,-0.4185E+00, &
    0.8181E-01,-0.4101E+00,-0.1836E+00, 0.2666E-01, 0.0000E+00, &
    0.6512E+02,-0.9112E+01,-0.3430E+01,-0.4870E+00, 0.3890E-01, &
    0.1031E+00, 0.4424E+00, 0.7264E-01, 0.3514E-02, 0.0000E+00, &
    0.6242E+02,-0.7553E+01,-0.4117E+01,-0.1154E+01, 0.2788E+00, &
    0.1618E+00, 0.5605E+00, 0.6024E-01,-0.1659E-01, 0.0000E+00, &
    0.3255E+02,-0.1698E+02,-0.3049E+01, 0.2481E+01, 0.9171E+00, &
    0.1086E+00, 0.2261E+00, 0.4896E-01,-0.3257E-01, 0.0000E+00/

   wvp = (wv-30.0)/20.0
   up = (u-6.0)/3.0

   x = wvp
   x2 = x*x
   x3 = x2*x
   y = up
   y2 = y*y
   y3 = y2*y
   xy = x*y
   x2y = x2*y
   xy2 = x*y2

   tb = a(1,ich) + a(2,ich)*x + a(3,ich)*y + a(4,ich)*x2  &
        + a(5,ich)*xy + a(6,ich)*y2 + a(7,ich)*x3         &
        + a(8,ich)*x2y + a(9,ich)*xy2 + a(10,ich)*y3

end subroutine tb_clear

   !**************************************************************
   ! This subroutine accepts SSM/I brightness temperatures and returns
   ! column water vapor in mm (or kg per m**2).  The algorithm is described
   ! in the document by Petty (1997) accompanying this FORTRAN module.
   ! 
   ! Inputs:   tb(5) = 19v,t19h,t22v,t37v,t37h
   !                --- SSM/I brightness temperatures (K) 
   !  
   ! Output:   wv  --- column water vapor (mm)
   ! 
   ! 

subroutine PettyVapor(tb,wv)
   implicit none
   real tb(*),wv

   real t19v,t19h,t22v,t37v,t37h
   real tbold(5),wvold
   save tbold,wvold

   logical flag
   real, parameter :: BAD = -1000.0
   integer :: i
   real :: del, t1,t2,t3,wv2,wv3

   ! check for same TBs used in previous call in order to avoid
   ! redundant calculations

   flag = .false. 
   do i=1,5
      if (tbold(i) .ne. tb(i)) flag = .true.
      tbold(i) = tb(i)
   end do
   if (.not. flag) then
      wv = wvold
      return
   end if

   ! initialize brightness temperatures

   t19v = tb(1)
   t19h = tb(2)
   t22v = tb(3)
   t37v = tb(4)
   t37h = tb(5)

   ! check for land, ice, and heavy precipitation

   del = t22v-t19v
   if (del .lt. 4.0) then
      wv = BAD
      wvold = wv
      return
   end if

   ! special check for sea ice

   t1 = 260.0 - 2.0*del
   t2 = 270.0 - (120.0/15.0)*del
   t3 = 130 + (20./15.)*del

   if (t19h .lt. t1 .and. t19h .lt. t2 .and. t19h .gt. t3) then
      wv = BAD
      wvold = wv
      return
   end if

   ! check for abnormally warm brightness temperatures
   if (t22v .gt. 285.0 .or.  &
        t19h .gt. 285.0 .or. &
        t37h .gt. 285.0) then
      wv = BAD
      wvold = wv
      return
   end if

   ! compute water vapor from first-stage regression algorithm

   wv = 0.1670E+03  &
        + 0.4037E+01*log(295-t19h) &
        - 0.5322E+02*log(295-t22v) &
        + 0.1296E+02*log(295-t37h)

   ! apply polynomial correction to eliminate biases at high
   ! and low end of range

   wv2 = wv*wv
   wv3 = wv*wv2
   !      wv = -0.2079E+01 + 0.1366E+01*wv+  &
   !           -0.1504E-01*wv2+0.1639E-03*wv3
   wv = -2.079 + 1.366*wv-0.01504*wv2+0.0001639*wv3
   wvold = wv

end subroutine PettyVapor

   !**********************************************************************
   ! This subroutine accepts SSM/I brightness temperatures and returns
   ! surface wind speed in  m/sec.  The algorithm is described
   ! in the document accompanying this FORTRAN module.
   ! 
   ! Inputs:   t19v,t19h,t22v,t37v,t37h  
   !                --- SSM/I brightness temperatures (K) 
   !           theta -- sensor viewing angle
   !           sst   -- sea surface temperature (deg. C)
   !           dsst   -- uncertainty in the above SST (deg. C)
   !                     (if dsst > 2.8 C, then an algorithm is
   !                      used that doesn't require SST)
   !
   ! Output:   u  ---- column water vapor (mm)
   ! 
   ! 

subroutine PettyWind(tb,theta,sst,dsst,u)
   implicit none
   real tb(5),theta,sst,dsst,u

   logical flag
   real t19v,t19h,t22v,t37v,t37h
   real, parameter :: BAD = -1000.0

   real wv, alw37,errwv,p37clr,errlw,erru,p37,wv2,wv3
   integer iaold,i,ia
   real tbold(5),uold
   save tbold,uold,iaold

   ! choose wind algorithm to use, depending on how well SST is
   ! known

   if (abs(dsst).gt.2.8.or.sst.lt.-4.0.or.sst.gt.35.0) then
      ia = 2
   else
      ia = 1
   end if

   ! check for same TBs used in previous call in order to avoid
   ! redundant calculations

   flag = .false. 

   if (ia .ne. iaold)  flag = .true.
   iaold = ia

   do i=1,5
      if (tbold(i) .ne. tb(i)) flag = .true.
      tbold(i) = tb(i)
   end do
   if (.not. flag) then
      u = uold
      return
   end if

   ! initialize brightness temperatures

   t19v = tb(1)
   t19h = tb(2)
   t22v = tb(3)
   t37v = tb(4)
   t37h = tb(5)

   ! get estimate of column water vapor for use in 
   ! cross talk corrections

   call PettyVapor(tb,wv)

   ! check for exceptionally high or bad water vapor values
   if (wv .gt. 68.0 .or. wv .lt. 0.0) then
      u = BAD
      uold = u
      return
   end if


   call ualg(ia,theta,sst,tb,u)

   ! calculate 37 GHz liquid water using simple algorithm
   p37clr = 77.68 -.1782*wv -1.1546*u +  &
        0.001838*wv*u -.003160*wv*wv + 0.01127*u*u
   p37 = (t37v-t37h)/p37clr
   alw37 = -(1.0/0.8614)*log(p37)

   ! check for cloud water exceeding threshold

   if (alw37 .gt. 0.5) then
      u = BAD
      uold = u
      return
   end if

   ! select appropriate corrections
   if (ia .eq. 1) then

      ! water vapor cross-talk correction
      wv2 = wv*wv
      wv3 = wv*wv2
      errwv = -0.7173 + 0.1160*wv -0.4238E-02*wv2 + 0.4372E-04*wv3
      u = u - errwv

      ! cloud liquid water cross-talk correction
      if (alw37 .lt.0.08) then
         errlw =  -6.3889*alw37 + 0.36111
      else
         errlw = 0.83333*alw37 - 0.21667
      end if
      u = u - errlw

      ! check for unphysical values
      if (u .lt. -2.0) then
         u = BAD
         uold = u
         return
      end if

      ! apply correction for non-linearity
      if (u .lt. 2.5) then
         erru  =  -1.267 + 0.7142*u -0.8395E-01*u*u
      else if (u .gt. 10.0 ) then
         erru  =  3.243 -0.3478*u + 0.3679E-02*u*u
      else
         erru = 0.0
      end if
      u = u - erru

   else if (ia .eq. 2) then

      ! water vapor cross-talk correction
      wv2 = wv*wv
      wv3 = wv*wv2
      errwv = 0.2945 + 0.01270*wv -0.001654*wv2+ 0.2596E-04*wv3
      u = u - errwv

      ! cloud liquid water cross-talk correction
      if (alw37 .lt.0.0) then
         errlw =  4.0*alw37 + 0.4
      else if (alw37 .ge. 0.0 .and. alw37 .lt. 0.08) then
         errlw = -8.125*alw37 + 0.4
      else
         errlw = -0.25
      end if
      u = u - errlw

      ! check for unphysical values
      if (u .lt. -2.0) then
         u = BAD
         uold = u
         return
      end if

      ! apply correction for non-linearity
      if (u .lt. 2.5) then
         erru = -1.219+ 0.9327*u -0.1852*u*u
      else if (u .gt. 10.0 ) then
         erru = 3.360-0.3918*u + 0.8121E-02*u*u
      else
         erru = 0.0
      end if
      u = u - erru
   end if
   uold = u

end subroutine PettyWind

subroutine ualg(ia,theta,sst,tb,u)

   !***********************************************************
   !  This subroutine implements the actual first stage wind speed algorithm.
   !  When IA=1, an algorithm is employed that makes use of SST.
   !  When IA=2, the SST information is not used.

   implicit none

   integer ia,i
   real theta,sst,tb(5),u

   real coeff1(8)
   real coeff2(8)
   data coeff1/ 0.1862E+03,0.9951E-01,0.0,  &
        -0.1829E-01, -0.1438E+01,0.7029,    &
        0.2329E+01,0.1687/
   data coeff2/0.1719E+03,0.2827, 0.0, -0.2549E-01, &
        -0.1473E+01, 0.6425, 0.2045E+01, 0.0/

   ! if theta is out of range, then substitute default value
   if (theta .lt. 50.0 .or. theta .gt. 55.0) theta = 53.1

   if (ia .eq. 1) then
      u = coeff1(1)
      do i=1,5
         u = u + coeff1(i+1)*tb(i)
      end do
      u = u + coeff1(7)*(theta-53.0)
      u = u + coeff1(8)*sst
   else if (ia .eq. 2) then
      u = coeff2(1)
      do i=1,5
         u = u + coeff2(i+1)*tb(i)
      end do
      u = u + coeff2(7)*(theta-53.0)
   else
      write (0,*) "Invalid algorithm ID"
   end if

end subroutine ualg