var/da/da_ssmi/read_ssmi.inc

   1 subroutine filter(tb19v, tb19h, tb22v, tb37v,          &
   2                      tb37h, tb85v, tb85h, ipass, iprecip, xlat )
   3    implicit none
   4
   5    real    ,  intent (in)  :: tb19v, tb19h, tb22v, tb37v, &
   6                               tb37h, tb85v, tb85h, xlat
   7    real                     :: tb(7)
   8    logical ,  intent (inout):: ipass
   9    integer ,  intent (inout):: iprecip
  10
  11    integer                 :: ipass1, ipass2
  12    real                    :: sst, dsst
  13    real                    :: theta,tc,wv,u,alw19,alw37,alw85
  14
  15    !  write (unit=stdout,fmt=*) 'tbdata',tb19v, tb19h, tb22v, tb37v,tb37h, tb85v, tb85h
  16    ! tc : cloud top temperature (C)
  17
  18    theta = 53.1
  19    tc = 10.0
  20    sst = 23.0
  21    dsst = 1000.0
  22
  23    ipass1 = 0
  24    ipass2 = 0
  25
  26    tb(1) = tb19v
  27    tb(2) = tb19h
  28    tb(3) = tb22v
  29    tb(4) = tb37v
  30    tb(5) = tb37h
  31    tb(6) = tb85v
  32    tb(7) = tb85h
  33
  34    ! upper and lower boundary 19h (90,280), others (100,280)
  35
  36    if (tb(1) .lt. 280 .and. tb(1) .gt. 100. .and. &
  37        tb(2) .lt. 280 .and. tb(2) .gt.  90. .and. &
  38        tb(3) .lt. 280 .and. tb(3) .gt. 100. .and. &
  39        tb(4) .lt. 280 .and. tb(4) .gt. 100. .and. &
  40        tb(5) .lt. 280 .and. tb(5) .gt. 100. .and. &
  41        tb(6) .lt. 280 .and. tb(6) .gt. 100. .and. &
  42        tb(7) .lt. 280 .and. tb(7) .gt. 100. .and. &
  43
  44        !    horizontal always much less than vertical
  45
  46        tb(1) .gt. tb(2) .and. &
  47        tb(4) .gt. tb(5) .and. &
  48        tb(6) .gt. tb(7) .and. &
  49
  50        ! T19V always at least 4 K less than  T22V, except in heavy rain
  51
  52        tb(3)-tb(1) .gt. 4.) then
  53       ipass1 = 1
  54    end if
  55
  56    ! second check
  57
  58    if (ipass1 .eq. 1) then
  59
  60       call pettyalgos(tb,theta,tc,sst,dsst,          &
  61                           wv,u,alw19,alw37,alw85,iprecip, xlat)
  62
  63          if (wv    .gt. -999. .and. u     .gt. -999. .and. &
  64              alw19 .gt. -999. .and. alw37 .gt. -999. .and. &
  65              alw85 .gt. -999. ) then
  66             ipass2 = 1
  67             ! write(unit=stdout,*)a1,a2,lat,long,theta,tb(1),tb(2),tb(3), &
  68             !    tb(4),tb(5),tb(6),tb(7)
  69          end if
  70
  71       ! third check
  72
  73       ! if (ipass2 .eq. 1) then
  74       !    if (wv    .ge. 5.    .and. wv .le. 70 .and.     &
  75       !        u     .ge. -5    .and. u  .le. 30 .and.     &
  76       !        alw19 .ge. -0.01 .and. alw19 .le. 0.5 .and. &
  77       !        alw37 .ge. -0.01 .and. alw37 .le. 0.5 .and. &
  78       !        alw85 .ge. -0.01 .and. alw85 .le. 0.5) then
  79       !       write(unit=stdout,*)a1,a2,lat,long,theta,tb(1),tb(2),tb(3), &
  80       !          tb(4),tb(5),tb(6),tb(7)
  81       !       ipass = .true.
  82       !    end if
  83       ! end if
  84       ! if (iprecip .eq. 1) then
  85       !    ipass = .true.
  86       !    write(unit=stdout,*) 'precipitation',lat,long,wv,u,alw19,alw37,alw85,iprecip
  87       ! else
  88       !    ipass = .true.
  89       ! end if
  90    end if
  91
  92    ! write(unit=stdout,*) 'pass=',ipass1,ipass2,iprecip,ipass
  93
  94 end subroutine filter
  95
  96    !***************************************************************
  97    !                    NOTICE
  98    !
  99    ! This package of subroutines has been developed by Grant W. Petty
 100    ! (Purdue University) for the retrieval of column water vapor, column
 101    ! cloud liquid water, and surface wind speed from SSM/I brightness
 102    ! temperatures.  The algorithm versions contained herein were developed
 103    ! in part using data sets provided by NASDA/EORC, Japan and are
 104    ! submitted to NASDA for evaluation as part of ADEOS-2/AMSR algorithm
 105    ! development activities.  Because of the need to meet a submission
 106    ! deadline, they have been subject only to limited testing, and it is
 107    ! possible that undetected bugs or problems remain in the code that
 108    ! follows.
 109    !
 110    ! Problems or questions should be directed to
 111    !       Grant W. Petty
 112    !       Earth and Atmospheric Sciences Dept.
 113    !       Purdue University
 114    !       West Lafayette, in 47907-1397
 115    !       (765) 494-2544
 116    !       gpetty@rain.atms.purdue.edu
 117    !
 118    ! Version:   July 13, 1997
 119    ! Copyright 1997, Grant W. Petty
 120    !
 121    ! Permission is granted to all scientific users to utilize
 122    ! and/or test this set of algorithms without restriction, provided only
 123    ! that
 124    !       (1) they are not used for commercial Purposes without the
 125    !          author's permission
 126    !
 127    !       (2) the author is given due credit for their development
 128    !
 129    !       (3) any modifications to the algorithms by 3d parties must be
 130    !          clearly identified as such and distinguished from the author's
 131    !          original code
 132    !
 133    !********************************************************************
 134    !      GENERAL COMMENTS ON USAGE
 135    !
 136    ! The algorithms that follow can act on single-pixel values of brightness
 137    ! temperature and derived variables, such as water vapor, wind speed,
 138    ! etc.  However, the way they are normally used at Purdue University is
 139    ! as follows:
 140    !
 141    !  1    Retrieve Water Vapor and Surface Wind Speed on a pixel-by-pixel basis
 142    !
 143    !  2    Apply mild spatial smoothing to derived wind speed and water vapor
 144    !       fields
 145    !
 146    !  3    Used filtered values of wind speed and water vapor to compute
 147    !       cloud-free polarization differences
 148    !
 149    !  4    Compute liquid water from unsmoothed, full-resolution brightness
 150    !       temperatures, normalized by smoothed cloud-free polarizations.
 151    !
 152    ! The smoothing step described above is optional and should not have a
 153    ! strong effect on the results outside of precipitation.  Normally the
 154    ! smoothing step is necessary only for interpolating water vapor and
 155    ! wind speed into areas of precipitation, where normalized polarizations
 156    ! at 19, 37 and 85 GHz are needed for precipitation retrievals.
 157    !
 158    !
 159    !********************************************************************
 160    ! This subroutine accepts SSM/I  brightness temperatures and returns
 161    ! all over-ocean variables using the algorithms of Petty (1997, unpublished)
 162    ! It, and the algorithms it calls, are generally intended for use only outside
 163    ! of precipitation, except for the water vapor algorithm, which is believed to
 164    !  tolerate significant precipitation without serious errors.
 165    !
 166    ! Inputs:   tb(7) = 19v,t19h,t22v,t37v,t37h,t85v,t85h
 167    !                --- SSM/I brightness temperatures (K)
 168    !           theta -- sensor viewing angle
 169    !           tc  --- assumed cloud temperature in degrees Celsius
 170    !           Note: If TC not known, set to a reasonable
 171    !                 average value, such as 8.0 degrees
 172    !
 173    !           sst   -- sea surface temperature (deg. C)
 174    !           dsst   -- uncertainty in the above SST (deg. C)
 175    !                     (if dsst > 2.8 C, then an wind speed algorithm is
 176    !                      used that doesn't require SST)
 177    !
 178    ! Output:   wv     ---- column water vapor (kg/m**2) from 19.35 GHz
 179    !           u      ---- surface wind speed (m/sec)
 180    !           alw19  ---- column cloud liquid water (kg/m**2) from 19.35 GHz
 181    !           alw37  ---- column cloud liquid water (kg/m**2) from 37 GHz
 182    !           alw85  ---- column cloud liquid water (kg/m**2) from 85.5 GHz
 183    !           iprecip --- set to 1 if precipitation flag invoked, else 0
 184    !
 185    !
 186    ! NOTE: It is not expected that liquid water values from different SSM/I
 187    ! frequencies will always closely agree, on account of significant
 188    ! differences in spatial resolution, sensitivity, and dynamic range.
 189    ! For best comparisons, average results to a common spatial resolution
 190    ! first.  19 GHz channels have the least sensitivity to thin clouds; 85
 191    ! GHz channels are most severely affected by precipitation and/or strong
 192    ! inhomogeneity in cloud fields.
 193    !
 194    ! Also, liquid water values in excess of 0.5 kg/m**2 are generally
 195    ! flagged as precipitating, and a value of MISSinG is returned.  There
 196    ! is no known method for separating cloud water from precipitation
 197    ! water.  Furthermore, any attempt to retrieve total liquid water path
 198    ! when precipitation is present will require a considerably more
 199    ! sophisticated algorithm than the ones provided here.
 200    !******************************************************************
 201
 202 subroutine pettyalgos(tb,theta,tc,sst,dsst, &
 203    wv,u,alw19,alw37,alw85,iprecip,xlat)
 204
 205    implicit none
 206
 207    real tb(7),theta,tc,wv,u,alw19,alw37,alw85,xlat, sst,dsst
 208    real, parameter :: BAD = -1000.0
 209    integer iprecip
 210
 211    ! get water vapor
 212    !     write (unit=stdout,fmt=*) 'tb=',tb
 213    call pettyvapor(tb,wv)
 214    !     write (unit=stdout,fmt=*) 'wv=',wv
 215
 216    ! get surface wind speed
 217    call pettywind(tb,theta,sst,dsst,u)
 218
 219    ! get column cloud liquid water, using 3 different algorithms
 220    call pettylwp(tb,theta,1,tc,alw19,iprecip,xlat)
 221    call pettylwp(tb,theta,2,tc,alw37,iprecip,xlat)
 222    call pettylwp(tb,theta,3,tc,alw85,iprecip,xlat)
 223    if (alw19 .eq. BAD .or. alw19 .eq. BAD .or. alw85 .eq. BAD) iprecip=1
 224    !     if (iprecip .eq. 1) write (unit=stdout,fmt=*) 'iprecip B=',iprecip,'alw=',alw19,alw37,alw85
 225
 226 end subroutine pettyalgos
 227
 228
 229    !********************************************************************
 230    ! This subroutine accepts SSM/I brightness temperatures and returns
 231    ! cloud liquid water in kg/m**2.
 232    !
 233    ! Inputs:   tb(7) = 19v,t19h,t22v,t37v,t37h,t85v,t85h
 234    !                --- SSM/I brightness temperatures (K)
 235    !           theta -- sensor viewing angle
 236    !           ifreq -- frequency to use for computing cloud liquid water
 237    !                    (1 = 19.35 GHz  2 = 37 GHz  3 = 85.5 GHz)
 238    !           tc  --- assumed cloud temperature in degrees Celsius
 239    !           Note: If not known, set to a reasonable value, such as 8.0 degrees
 240    !
 241    ! Output:   alw  ---- column water vapor (mm)
 242    !           iprecip --- set to 1 if precipitation flag invoked, else 0
 243
 244 subroutine pettylwp(tb,theta,ifreq,tc,alw,iprecip,xlat)
 245
 246    implicit none
 247    real tb(7),theta,alw,tc,xlat
 248    integer ifreq,iprecip
 249    real, parameter :: BAD = -1000.0
 250    real t19v,t19h,t22v,t37v,t37h,t85v,t85h,MISSinG,pclr
 251    real :: p,s85,v0,wv,u
 252    integer :: ich
 253
 254    iprecip = 0
 255    MISSinG = -8888.
 256    pclr=0.
 257    if (ifreq.lt.1 .or. ifreq.gt.3) then
 258       pause "BAD VALUE OF ifREQ in PettyLWP"
 259       alw = BAD
 260       return
 261    end if
 262
 263    ! initialize brightness temperatures
 264
 265    t19v = tb(1)
 266    t19h = tb(2)
 267    t22v = tb(3)
 268    t37v = tb(4)
 269    t37h = tb(5)
 270    t85v = tb(6)
 271    t85h = tb(7)
 272
 273    ! estimate the clear-sky polarization difference for the scene
 274
 275    call clearpol(tb,theta,ifreq,pclr,iprecip)
 276
 277    ! if value of Tc is invalid, then set to reasonable intermediate value
 278    if (tc .lt. -20.0 .or. tc .gt. 40.0) then
 279       tc = 8.0
 280    end if
 281
 282    ! calculate normalized polarization
 283    if (ifreq .eq. 1) then
 284       p = (t19v-t19h)/pclr
 285    else if (ifreq .eq. 2) then
 286       p = (t37v-t37h)/pclr
 287    else if (ifreq .eq. 3) then
 288       p = (t85v-t85h)/pclr
 289
 290       ! at 85.5 GHz, there could be errors due to ice, so check S85
 291       call PettyVapor(tb,wv)
 292       if (wv .eq. BAD) then
 293          alw = BAD
 294          return
 295       end if
 296       call PettyWind(tb,theta,-99.9,1000.0,u)
 297       if (u .eq. BAD) then
 298          u = 5.0
 299          alw = BAD
 300       end if
 301      ich = 6
 302      call tb_clear(ich,u,wv,v0)
 303      s85 = p*v0 + (1.0-p)*t_kelvin - t85v
 304      if (s85 .gt. 10.0 .or. (t85v-t85h) .lt. 5.0) p = MISSinG
 305    end if
 306
 307    ! convert normalized polarization to LWP
 308    if (p .ne. MISSinG .and. p .lt. 1.4 .and. p .gt. 0.1) then
 309       call P2LWP(ifreq,p,tc,theta,alw,xlat)
 310    else
 311       alw = BAD
 312    end if
 313
 314 end subroutine PettyLWP
 315
 316    !*********************************************************
 317    ! This routine computes liquid water path from normalized polarization
 318    ! at 19.35, 37, or 85.5 GHz  (indicated by ifREQ)
 319
 320 subroutine p2lwp(ifreq,p,tc,theta,alw,xlat)
 321
 322    implicit none
 323    integer ifreq
 324    real p,alw,xlat,tc,theta
 325    real alpha(3),threshhold
 326    data alpha/2.08,2.05,1.78/
 327    real, parameter :: BAD = -1000.0
 328    real :: ext, costheta
 329    integer :: iprecip
 330
 331    ! check for extreme value due to precipitation
 332    if (p .lt. 0.01) then
 333       alw = BAD
 334       iprecip = 1
 335       return
 336    end if
 337
 338    !  get liquid water mass extinction coefficient
 339    ext = alw_ext(ifreq,tc)
 340
 341    ! compute colum cloud liquid water using Beer's Law
 342    costheta = cos(theta*pi/180.0)
 343    alw = (-costheta/(alpha(ifreq)*ext))*log(p)
 344
 345    ! check for precipitation
 346    threshhold=0.5
 347    !     write (unit=stdout,fmt=*) 'xlat=',xlat
 348    if (xlat .gt. 20) threshhold=0.3
 349    if (xlat .gt. 30) threshhold=0.1
 350    !     if (xlat .gt. 40) threshhold=0.25
 351    !     if (alw .gt. 0.5) then
 352    if (alw .gt. threshhold) then
 353       alw = BAD
 354       iprecip = 1
 355    end if
 356
 357 end subroutine p2lwp
 358
 359    !***********************************************************************
 360    ! This subroutine accepts SSM/I brightness temperatures from a scene and returns
 361    ! the predicted cloud-free polarization difference for that scene
 362    !
 363    ! Inputs:   tb(7) = 19v,t19h,t22v,t37v,t37h,t85v,t85h
 364    !                --- SSM/I brightness temperatures (K)
 365    !           theta -- sensor viewing angle
 366    !           ifreq -- frequency to use for computing cloud liquid water
 367    !                    (1 = 19.35 GHz  2 = 37 GHz  3 = 85.5 GHz)
 368    !
 369    ! Output:   pclr  ---- cloud free polarization difference (K)
 370    !           iprecip --- set to 1 if precipitation flag invoked, else 0
 371
 372 subroutine ClearPol(tb,theta,ifreq,pclr,iprecip)
 373
 374    implicit none
 375    real tb(7),theta,pclr
 376    integer ifreq,iprecip
 377    real, parameter :: BAD = -1000.0
 378    real :: wv,u
 379    integer :: ich
 380
 381    ! check for valid ifREQ
 382    if (ifreq.lt.1 .or. ifreq.gt.3) then
 383       pause 'BAD VALUE OF ifREQ in ClearPol'
 384       pclr = BAD
 385       return
 386    end if
 387
 388    iprecip = 0
 389
 390    ! get water vapor estimate.  If this value is returned as 'bad', then
 391    ! either the brightness temperatures are not valid for over ocean or
 392    ! else there is unusually heavy precipitation
 393
 394    call PettyVapor(tb,wv)
 395    if (wv .lt. 0.0) then
 396       pclr = BAD
 397       return
 398    end if
 399
 400    ! get wind speed estimate.  If this value is returned as 'bad', assume
 401    ! that it is due to precipitation, and substitute a reasonable 'average'
 402    ! wind speed.
 403
 404    call PettyWind(tb,theta,-99.9,1000.0,u)
 405    if (u .eq. BAD) then
 406       iprecip = 1
 407       u = 5.0
 408    end if
 409
 410    ! now that we finally have estimates of water vapor and wind speed,
 411    ! compute the clear-sky polarization difference at the desired frequency
 412
 413    ich = ifreq + 7
 414    call tb_clear(ich,u,wv,pclr)
 415
 416          end subroutine ClearPol
 417
 418    !**********************************************************
 419    ! The following function returns the microwave mass extinction coefficient
 420    ! of cloud water at 19.35, 37.0, or 85.5 GHz.
 421    !
 422    !     inputs:
 423    !           ifreq -- frequency to use for computing cloud liquid water
 424    !                    (1 = 19.35 GHz  2 = 37 GHz  3 = 85.5 GHz)
 425    !           tc  --- cloud temperature in degrees Celsius
 426    !
 427    !      returned value:  Mass extinction coefficient (m**2/kg)
 428    !
 429
 430 function alw_ext(ifreq,tc)
 431
 432    implicit none
 433
 434    integer, intent(in) :: ifreq
 435    real, intent(in) :: tc
 436
 437    real :: alw_ext
 438    real :: tc2,tc3
 439
 440    tc2 = tc*tc
 441    tc3 = tc2*tc
 442
 443    if (ifreq .eq. 1) then
 444       alw_ext = exp(-2.55-2.98e-2*tc-6.81e-4*tc2+3.35e-6*tc3)
 445    else if (ifreq .eq. 2) then
 446       alw_ext = exp(-1.35-0.0234*tc-1.22e-4*tc2+5.48e-6*tc3)
 447    else if (ifreq .eq. 3) then
 448       alw_ext = exp(-0.0713-8.00e-3*tc-1.62e-4*tc2+1.18e-6*tc3)
 449    end if
 450
 451 end function alw_ext
 452
 453    !*****************************************************************************
 454    ! The following routine returns approximate clear-sky SSM/I brightness temperature
 455    ! over the ocean as a function of column water vapor and surface wind speed.
 456    ! This implementation is a bivariate polynomial fit to model computed brightness
 457    ! temperatures (see Petty (1992,1994)), where the model was in turn carefully calibrated
 458    ! empirically against 200,000 cloud-free SSM/I pixels.
 459    ! This routine assumes a viewing angle not too different from 53.4 degrees, and
 460    ! it assumes "typical" values for other ocean and atmospheric parameters.
 461    !
 462    !      inputs:
 463    !         ich = channel no. (1 = 19V, 2 = 19H, ..., 7 = 85V,
 464    !                            8 = 19V-19H, 9 = 37V-37H, 10 = 85V-85H)
 465    !         u = estimated surface wind speed (m/sec)
 466    !         wv = estimated column water vapor (kg/m**2)
 467    !
 468    !      output:
 469    !         tb = brightness temperature
 470
 471 subroutine tb_clear(ich,u,wv,tb)
 472    implicit none
 473    integer, intent(in)  :: ich
 474    real,    intent(in)  :: u
 475    real,    intent(in)  :: wv
 476    real,    intent(out) :: tb
 477
 478    real :: wvp,up,x,x2,x3,y,y2,y3,xy,x2y,xy2
 479
 480    real a(10,10)
 481    data a/ &
 482     0.1984E+03, 0.1674E+02, 0.1998E+00,-0.3868E+00,-0.1818E+00, &
 483     0.7065E-01,-0.8148E+00, 0.6116E-01, 0.1400E-01, 0.0000E+00, &
 484     0.1334E+03, 0.2587E+02, 0.3819E+01,-0.7762E-01, 0.1564E+00, &
 485     0.4968E-01,-0.1232E+01,-0.6779E-01,-0.3083E-01, 0.0000E+00, &
 486     0.2301E+03, 0.2865E+02, 0.5589E+00,-0.4668E+01,-0.7274E-01, &
 487     0.1043E-01,-0.6731E+00, 0.8141E-02,-0.1081E-01, 0.0000E+00, &
 488     0.2126E+03, 0.1103E+02,-0.2488E+00, 0.8614E+00,-0.1298E+00, &
 489     0.4531E-01,-0.8315E+00, 0.5104E-01, 0.1785E-01, 0.0000E+00, &
 490     0.1500E+03, 0.1864E+02, 0.3950E+01, 0.2095E+01,-0.1743E+00, &
 491     0.1170E+00,-0.1345E+01,-0.8584E-01, 0.8143E-02, 0.0000E+00, &
 492     0.2589E+03, 0.1706E+02,-0.7729E+00,-0.2905E+01, 0.2821E+00, &
 493     0.3375E-02,-0.4366E+00,-0.1150E+00, 0.1881E-01, 0.0000E+00, &
 494     0.2265E+03, 0.3363E+02, 0.2050E+01,-0.5701E+01,-0.4185E+00, &
 495     0.8181E-01,-0.4101E+00,-0.1836E+00, 0.2666E-01, 0.0000E+00, &
 496     0.6512E+02,-0.9112E+01,-0.3430E+01,-0.4870E+00, 0.3890E-01, &
 497     0.1031E+00, 0.4424E+00, 0.7264E-01, 0.3514E-02, 0.0000E+00, &
 498     0.6242E+02,-0.7553E+01,-0.4117E+01,-0.1154E+01, 0.2788E+00, &
 499     0.1618E+00, 0.5605E+00, 0.6024E-01,-0.1659E-01, 0.0000E+00, &
 500     0.3255E+02,-0.1698E+02,-0.3049E+01, 0.2481E+01, 0.9171E+00, &
 501     0.1086E+00, 0.2261E+00, 0.4896E-01,-0.3257E-01, 0.0000E+00/
 502
 503    wvp = (wv-30.0)/20.0
 504    up = (u-6.0)/3.0
 505
 506    x = wvp
 507    x2 = x*x
 508    x3 = x2*x
 509    y = up
 510    y2 = y*y
 511    y3 = y2*y
 512    xy = x*y
 513    x2y = x2*y
 514    xy2 = x*y2
 515
 516    tb = a(1,ich) + a(2,ich)*x + a(3,ich)*y + a(4,ich)*x2  &
 517         + a(5,ich)*xy + a(6,ich)*y2 + a(7,ich)*x3         &
 518         + a(8,ich)*x2y + a(9,ich)*xy2 + a(10,ich)*y3
 519
 520 end subroutine tb_clear
 521
 522    !**************************************************************
 523    ! This subroutine accepts SSM/I brightness temperatures and returns
 524    ! column water vapor in mm (or kg per m**2).  The algorithm is described
 525    ! in the document by Petty (1997) accompanying this FORTRAN module.
 526    !
 527    ! Inputs:   tb(5) = 19v,t19h,t22v,t37v,t37h
 528    !                --- SSM/I brightness temperatures (K)
 529    !
 530    ! Output:   wv  --- column water vapor (mm)
 531    !
 532    !
 533
 534 subroutine PettyVapor(tb,wv)
 535    implicit none
 536    real tb(*),wv
 537
 538    real t19v,t19h,t22v,t37v,t37h
 539    real tbold(5),wvold
 540    save tbold,wvold
 541
 542    logical flag
 543    real, parameter :: BAD = -1000.0
 544    integer :: i
 545    real :: del, t1,t2,t3,wv2,wv3
 546
 547    ! check for same TBs used in previous call in order to avoid
 548    ! redundant calculations
 549
 550    flag = .false.
 551    do i=1,5
 552       if (tbold(i) .ne. tb(i)) flag = .true.
 553       tbold(i) = tb(i)
 554    end do
 555    if (.not. flag) then
 556       wv = wvold
 557       return
 558    end if
 559
 560    ! initialize brightness temperatures
 561
 562    t19v = tb(1)
 563    t19h = tb(2)
 564    t22v = tb(3)
 565    t37v = tb(4)
 566    t37h = tb(5)
 567
 568    ! check for land, ice, and heavy precipitation
 569
 570    del = t22v-t19v
 571    if (del .lt. 4.0) then
 572       wv = BAD
 573       wvold = wv
 574       return
 575    end if
 576
 577    ! special check for sea ice
 578
 579    t1 = 260.0 - 2.0*del
 580    t2 = 270.0 - (120.0/15.0)*del
 581    t3 = 130 + (20./15.)*del
 582
 583    if (t19h .lt. t1 .and. t19h .lt. t2 .and. t19h .gt. t3) then
 584       wv = BAD
 585       wvold = wv
 586       return
 587    end if
 588
 589    ! check for abnormally warm brightness temperatures
 590    if (t22v .gt. 285.0 .or.  &
 591         t19h .gt. 285.0 .or. &
 592         t37h .gt. 285.0) then
 593       wv = BAD
 594       wvold = wv
 595       return
 596    end if
 597
 598    ! compute water vapor from first-stage regression algorithm
 599
 600    wv = 0.1670E+03  &
 601         + 0.4037E+01*log(295-t19h) &
 602         - 0.5322E+02*log(295-t22v) &
 603         + 0.1296E+02*log(295-t37h)
 604
 605    ! apply polynomial correction to eliminate biases at high
 606    ! and low end of range
 607
 608    wv2 = wv*wv
 609    wv3 = wv*wv2
 610    !      wv = -0.2079E+01 + 0.1366E+01*wv+  &
 611    !           -0.1504E-01*wv2+0.1639E-03*wv3
 612    wv = -2.079 + 1.366*wv-0.01504*wv2+0.0001639*wv3
 613    wvold = wv
 614
 615 end subroutine PettyVapor
 616
 617    !**********************************************************************
 618    ! This subroutine accepts SSM/I brightness temperatures and returns
 619    ! surface wind speed in  m/sec.  The algorithm is described
 620    ! in the document accompanying this FORTRAN module.
 621    !
 622    ! Inputs:   t19v,t19h,t22v,t37v,t37h
 623    !                --- SSM/I brightness temperatures (K)
 624    !           theta -- sensor viewing angle
 625    !           sst   -- sea surface temperature (deg. C)
 626    !           dsst   -- uncertainty in the above SST (deg. C)
 627    !                     (if dsst > 2.8 C, then an algorithm is
 628    !                      used that doesn't require SST)
 629    !
 630    ! Output:   u  ---- column water vapor (mm)
 631    !
 632    !
 633
 634 subroutine PettyWind(tb,theta,sst,dsst,u)
 635    implicit none
 636    real tb(5),theta,sst,dsst,u
 637
 638    logical flag
 639    real t19v,t19h,t22v,t37v,t37h
 640    real, parameter :: BAD = -1000.0
 641
 642    real wv, alw37,errwv,p37clr,errlw,erru,p37,wv2,wv3
 643    integer iaold,i,ia
 644    real tbold(5),uold
 645    save tbold,uold,iaold
 646
 647    ! choose wind algorithm to use, depending on how well SST is
 648    ! known
 649
 650    if (abs(dsst).gt.2.8.or.sst.lt.-4.0.or.sst.gt.35.0) then
 651       ia = 2
 652    else
 653       ia = 1
 654    end if
 655
 656    ! check for same TBs used in previous call in order to avoid
 657    ! redundant calculations
 658
 659    flag = .false.
 660
 661    if (ia .ne. iaold)  flag = .true.
 662    iaold = ia
 663
 664    do i=1,5
 665       if (tbold(i) .ne. tb(i)) flag = .true.
 666       tbold(i) = tb(i)
 667    end do
 668    if (.not. flag) then
 669       u = uold
 670       return
 671    end if
 672
 673    ! initialize brightness temperatures
 674
 675    t19v = tb(1)
 676    t19h = tb(2)
 677    t22v = tb(3)
 678    t37v = tb(4)
 679    t37h = tb(5)
 680
 681    ! get estimate of column water vapor for use in
 682    ! cross talk corrections
 683
 684    call PettyVapor(tb,wv)
 685
 686    ! check for exceptionally high or bad water vapor values
 687    if (wv .gt. 68.0 .or. wv .lt. 0.0) then
 688       u = BAD
 689       uold = u
 690       return
 691    end if
 692
 693
 694    call ualg(ia,theta,sst,tb,u)
 695
 696    ! calculate 37 GHz liquid water using simple algorithm
 697    p37clr = 77.68 -.1782*wv -1.1546*u +  &
 698         0.001838*wv*u -.003160*wv*wv + 0.01127*u*u
 699    p37 = (t37v-t37h)/p37clr
 700    alw37 = -(1.0/0.8614)*log(p37)
 701
 702    ! check for cloud water exceeding threshold
 703
 704    if (alw37 .gt. 0.5) then
 705       u = BAD
 706       uold = u
 707       return
 708    end if
 709
 710    ! select appropriate corrections
 711    if (ia .eq. 1) then
 712
 713       ! water vapor cross-talk correction
 714       wv2 = wv*wv
 715       wv3 = wv*wv2
 716       errwv = -0.7173 + 0.1160*wv -0.4238E-02*wv2 + 0.4372E-04*wv3
 717       u = u - errwv
 718
 719       ! cloud liquid water cross-talk correction
 720       if (alw37 .lt.0.08) then
 721          errlw =  -6.3889*alw37 + 0.36111
 722       else
 723          errlw = 0.83333*alw37 - 0.21667
 724       end if
 725       u = u - errlw
 726
 727       ! check for unphysical values
 728       if (u .lt. -2.0) then
 729          u = BAD
 730          uold = u
 731          return
 732       end if
 733
 734       ! apply correction for non-linearity
 735       if (u .lt. 2.5) then
 736          erru  =  -1.267 + 0.7142*u -0.8395E-01*u*u
 737       else if (u .gt. 10.0 ) then
 738          erru  =  3.243 -0.3478*u + 0.3679E-02*u*u
 739       else
 740          erru = 0.0
 741       end if
 742       u = u - erru
 743
 744    else if (ia .eq. 2) then
 745
 746       ! water vapor cross-talk correction
 747       wv2 = wv*wv
 748       wv3 = wv*wv2
 749       errwv = 0.2945 + 0.01270*wv -0.001654*wv2+ 0.2596E-04*wv3
 750       u = u - errwv
 751
 752       ! cloud liquid water cross-talk correction
 753       if (alw37 .lt.0.0) then
 754          errlw =  4.0*alw37 + 0.4
 755       else if (alw37 .ge. 0.0 .and. alw37 .lt. 0.08) then
 756          errlw = -8.125*alw37 + 0.4
 757       else
 758          errlw = -0.25
 759       end if
 760       u = u - errlw
 761
 762       ! check for unphysical values
 763       if (u .lt. -2.0) then
 764          u = BAD
 765          uold = u
 766          return
 767       end if
 768
 769       ! apply correction for non-linearity
 770       if (u .lt. 2.5) then
 771          erru = -1.219+ 0.9327*u -0.1852*u*u
 772       else if (u .gt. 10.0 ) then
 773          erru = 3.360-0.3918*u + 0.8121E-02*u*u
 774       else
 775          erru = 0.0
 776       end if
 777       u = u - erru
 778    end if
 779    uold = u
 780
 781 end subroutine PettyWind
 782
 783 subroutine ualg(ia,theta,sst,tb,u)
 784
 785    !***********************************************************
 786    !  This subroutine implements the actual first stage wind speed algorithm.
 787    !  When IA=1, an algorithm is employed that makes use of SST.
 788    !  When IA=2, the SST information is not used.
 789
 790    implicit none
 791
 792    integer ia,i
 793    real theta,sst,tb(5),u
 794
 795    real coeff1(8)
 796    real coeff2(8)
 797    data coeff1/ 0.1862E+03,0.9951E-01,0.0,  &
 798         -0.1829E-01, -0.1438E+01,0.7029,    &
 799         0.2329E+01,0.1687/
 800    data coeff2/0.1719E+03,0.2827, 0.0, -0.2549E-01, &
 801         -0.1473E+01, 0.6425, 0.2045E+01, 0.0/
 802
 803    ! if theta is out of range, then substitute default value
 804    if (theta .lt. 50.0 .or. theta .gt. 55.0) theta = 53.1
 805
 806    if (ia .eq. 1) then
 807       u = coeff1(1)
 808       do i=1,5
 809          u = u + coeff1(i+1)*tb(i)
 810       end do
 811       u = u + coeff1(7)*(theta-53.0)
 812       u = u + coeff1(8)*sst
 813    else if (ia .eq. 2) then
 814       u = coeff2(1)
 815       do i=1,5
 816          u = u + coeff2(i+1)*tb(i)
 817       end do
 818       u = u + coeff2(7)*(theta-53.0)
 819    else
 820       write (0,*) "Invalid algorithm ID"
 821    end if
 822
 823 end subroutine ualg
 824
 825
 826