C     This subroutine is used to calculate the Total Water Vapor of the atmosphere in the polar 
C     regions. It is adapted from the IDL program written by Jungang Miao in April 1998.
C     This FORTRAN version HAS NOT BEEN TESTED by now.
C       --- Jungang Miao 
C       --- April 4, 2003, Bremen, Germany.
C
C     The channel numbers of SSM/T2 used here are:
C       Channel 1:   92 GHz
C       Channel 2:  150 GHz
C       Channel 3:  183.31+/-7 GHz
C       Channel 4:  183.31+/-3 GHz
C       Channel 5:  183.31+/-1 GHz
C
C     input parameters:
C       n_scan:  total scan number of SSM/T2, i.e. the size of the second dimension of 2-D 
C                arrays of ch12, ch23, ch34, ch45, and twv.
C       ch12: a 2-D array with dimensions (28, n_scan), denoting the Tb difference of channels
C             1 and 2.
C       ch23: a 2-D array with dimensions (28, n_scan), denoting the Tb difference of channels
C             2 and 3.
C       ch34: a 2-D array with dimensions (28, n_scan), denoting the Tb difference of channels
C             3 and 4.
C       ch45: a 2-D array with dimensions (28, n_scan), denoting the Tb difference of channels
C             4 and 5.
C     output parameter:
C       twv:  a 2-D array with dimensions (28, n_scan), denoting the Total Water Vapor.



      SUBROUTINE vapor_ssmt2_200304 (n_scan, ch12, ch23, ch34, ch45, twv)
  
      INTEGER*2 n_scan
      REAL*4 dt23, dt34, dt45, fact_twv
      REAL*4 ch12(1:28,1:n_scan), ch23(1:28,1:n_scan), ch34(1:28,1:n_scan), 
     *       ch45(1:28,1:n_scan), twv(1:28,1:n_scan)
      REAL*4 scan_angle(1:28)
      DATA scan_angle/47.26, 43.50, 39.83, 36.22, 32.66, 29.13, 25.64, 
     *                22.18, 18.73, 15.31, 11.89,  8.49,  5.09,  1.70,
     *                 1.70,  5.09,  8.49, 11.89, 15.31, 18.73, 22.18,   
     *                 25.64, 29.13, 32.66, 36.22, 39.83, 43.40, 47.26/
 

      do 10 i=1, 28
         scan_angle(i)=cos(scan_angle/180.0*3.1415926)
 10      continue

         do 30 i=1, 28 
            do 20 j=1, n_scan
               dt23=ch23(i,j)
               dt34=ch34(i,j)
               dt45=ch45(i,j)
               call newalgo(dt23, dt34, dt45, fact_twv)
               twv(i,j)=fact_twv*scan_angle(i)
 20            continue
 30            continue

               return
               end

C     -------- subroutine newalgo --------

      subroutine newalgo(dt23, dt34, dt45, fact_twv)

      if ((dt45.gt.0).or.(dt34.gt.0)) then
        if ((dt34.gt.0).or.(dt23.gt.0)) then
           fact_twv=9999.9
           goto 40
        else 
           call algo_234(dt23, dt34, fact_twv)
           goto 40
        endif

      else
        call algo_345(dt34, dt45, fact_twv)

        if(fact_twv.gt.1.5) then
          call algo_234(dt23, dt34, fact_twv)
          goto 40
        endif

      endif

 40   return
      end
            
C     -------- subroutine algo_234 --------


      subroutine algo_234(dt23, dt34, fact_twv)

      INTEGER*2 nr
      REAL*4 px(1:4), py(1:4), c0(1:4), c1(1:4)
      DATA px/4.066, 1.451, 5.591, 3.525/
      DATA py/2.458, 0.744, 4.754, 0.384/
      DATA c0/2.041, 1.887, 2.010, 2.414/
      DATA c1/2.275, 1.789, 2.316, 2.110/

      nr=1
      fact_twv=c0(nr)+c1(nr)*alog((dt_23-py(nr))/(dt34-px(nr)))

      return
      end


C     -------- subroutine algo_345 --------


      subroutine algo_345(dt34, dt45, fact_twv)

      INTEGER*2 nr
      REAL*4 px(1:4), py(1:4), c0(1:4), c1(1:4)
      DATA px/2.556, 1.831, 1.378, 3.380/
      DATA py/1.370, 0.901, 0.343, 3.027/
      DATA c0/0.689, 0.685, 0.671, 0.693/
      DATA c1/0.723, 0.690, 0.565, 0.753/

      nr=1
      fact_twv=c0(nr)+c1(nr)*alog((dt_34-py(nr))/(dt45-px(nr)))

      return
      end
