DEWPOINT.f

Go to the documentation of this file.

      SUBROUTINE dewpoint( VP, TD)

       use ctlblk_mod, only: jsta, jend, im, spval, ista, iend
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!
!          NT IS THE TABLE SIZE
      integer,PARAMETER :: nt=2000
!...TRANSLATED BY FPP 3.00Z36 11/09/90  14:48:53  
!...SWITCHES: OPTON=I47,OPTOFF=VAE0
      real,intent(out) :: td(ista:iend,jsta:jend)
      real,intent(in) ::  vp(ista:iend,jsta:jend)
      real tdp(nt)
!jw
      integer nn,i,j,jnt
      real rvp1,rvp2,rt3,rvp3,rlog3,ra,rb,rapb,rtest,rnt,rdvp
      real rgs,rvp,rlvp,rn,rd,rch,rt,w1,w2
      real a,b,dntm1

      logical :: jcontinue=.true.

!          PREPARE THE TABLE (TDP=DEWPT AS FCN OF VAPOR PRESS).
!          RANGE IN CENTIBARS IS FROM RVP1 THRU RVP2
      rvp1  = 0.0001e0
      rvp2  = 10.e0
!          THE TRIPLE POINT
      rt3   = 273.16e0
!          VAPOR PRESS AT THE TRIPLE POINT
      rvp3  = 0.611e0
      rlog3 = log(rvp3)
!          (SPEC HT OF WATER -CSUBP OF VAPOR)/GAS CONST OF VAPOR.
      ra    = 5.0065e0
!          LATENT HEAT AT T3/(GAS CONST OF VAPOR * TRIPLE PT TEMP).
      rb    = 19.83923e0
      rapb  = ra + rb
!          CRITERION FOR CONVERGENCE OF NEWTON ITERATION
      rtest = 1.e-6
!MEB  RTEST=1.E-8  !  PROBABLY WON'T CONVERGE WITH 32-BIT AT THIS CRITERION
!
      rnt   = float(nt)
!          TABLE INCREMENT IN VAPOR PRESS
      rdvp  = (rvp2-rvp1)/(rnt-1.e0)
!          RGS WILL BE THE GUESSED VALUE OF (T3  /  DEWPOINT)
      rgs   = 1.e0
      rvp   = rvp1-rdvp
!
      DO 20 nn=1,nt
        rvp=rvp+rdvp
        rlvp=log(rvp)-rlog3-rapb
!     ***** ENTER NEWTON ITERATION LOOP
        jcontinue=.true.
        do while (jcontinue)
   10   rn=ra*log(rgs)-rapb*rgs-rlvp
!          THAT WAS VALUE OF FUNCTION
!          NOW GET ITS DERIVATIVE
        rd=(ra/rgs)-rapb
!          THE DESIRED CHANGE IN THE GUESS
        rch=rn/rd
        IF( abs(rch) < rtest ) jcontinue=.false.
!            NEED MORE ITERATIONS
          DO WHILE (abs(rch) >= rtest)
            rgs=rgs-rch
            EXIT
          ENDDO
        ENDDO
!          *****
!          HAVE ACCURATE ENUF VALUE OF RGS=T3/DEWPOINT.
   15   rt=rt3/rgs
        tdp(nn)=rt
!
   20 CONTINUE
!      PRINT 25,RVP1,RVP2,TDP(1),TDP(NT)
!  25  FORMAT(/'0', 'IN SUBROUTINE DEWPOINT, THE DEWPT TABLE ',
!    1             'HAS RVP1=', 1PE13.6, ', RVP2=', 1PE13.6,
!    2             ', TDP(1)=', 1PE13.6, ', AND TDP(NT)=',
!    3             1PE13.6, '.'/)
!           CONSTANTS FOR USING THE TABLE
      a     = 1./rdvp
      b     = 1. - a*rvp1
      dntm1 = float(nt) -.01
!
!X      END IF
!
!          *********** ENTER TO USE THE TABLE.  ************
!
!$omp parallel do private(i,j,w1,w2,jnt)
      DO j=jsta,jend
        DO i=ista,iend
          IF(vp(i,j)<spval)THEN
          w1  = min(max((a*vp(i,j)+b),1.0),dntm1)
          w2  = aint(w1)
          jnt = int(w2)
          td(i,j) = tdp(jnt) + (w1-w2)*(tdp(jnt+1)-tdp(jnt))
          ELSE
          td(i,j) = spval
          ENDIF
        ENDDO
      ENDDO
!
!
      RETURN
      END