C vdraft3.for 
C-----------------------------------------------------------------------------
!  6 Jan 2011 - Added entry point IPWINDg(modug)
C 17 Dec 2009 - Verified the model
C 16 Dec 2009 - Finished adding all of the necessary equations
C 11 Dec 2009 - Added subroutines psiD and psiE
C 30 Nov 2009 - Copied vdraft2 to vdraft3 and added ring vortexes at top
C-----------------------------------------------------------------------------
C  9 Nov 2009 - Added:     implicit double precision (A-H,O-Z)
C  6 Nov 2009 - Added:     double precision d1mach
C  6 Nov 2009 - Used D1mach instead of rmach for machine constants
C 19 Mar 2008 - Added the calculation of c1 and c2
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 10 Mar 2008: corrected length errors in blank common
C 6 Sept 2006: Fixed overflow error on exponential calculation (cosh)
C  8 Aug 2006: Correction
C  7 Aug 2006: Correctly did the stream function in spherical polar coordinates
C  7 Aug 2006: There is a discontinuity in the velocity at the center of the
C              updraft/downdraft
C  3 Aug 2006: Deleted some stuff.
C  3 Aug 2006: Added calculations for sin gamma = 0
C  2 Aug 2006: Using better formula for central earth angle
C  1 Aug 2006: Corrections
C 31 Jul 2006: Formulas for small gamma (GAM) and small beta
C 18 Jul 2006: Corrected PVHPG
C 17 Jul 2006: Corrected VR and PVRR
C 13 Jul 2006: Used dcosh instead of dsech
C 10 Jul 2006: Corrected $ to % (4 places)
C 10 Jul 2006: Corrected comma comment
C 27 Apr 2006: Added = to line 369
C 27 Apr 2006: vsiC1 -> psiC1, vsiC2 -> psiC2, pii -> psi, VDRAFT -> VDRAFT2
C 24 Apr 2006: Modular version broken into subroutines
C 24 Apr 2006: Version for the updraft being capped at a user-specified height
C 24 Apr 2006: Copied from vdraft.for
C-----------------------------------------------------------------------------
C 15 Jun 2005: Converted to perturbation model
C 31 May 2005: Converted back to old check for wrong model
C 28 May 2005: Commented out old programming
C 27 May 2005: Deleted old programming
C 27 May 2005: Commented out old programming
C 27 May 2005: Corrected PVHR1 and PVHR2
C 27 May 2005: Corrected SING1SQ and SING2SQ
C 26 May 2005: Corrected PCOSALPTH and PCOSALPPH
C 26 May 2005: Print out phi, theta, alpha
C 25 May 2005: Corrected cos alpha calculation plus derivatives.
C 24 May 2005: Added printing of gradients.
C 23 May 2005: Reverted to new formulas.
C 23 May 2005: Reverted to old alternate without errors
C 23 May 2005: Reverted to cos^2 alpha = 1.d-16 when negative.
C 23 May 2005: Reverted to cos^2 alpha = zero when negative.
C 23 May 2005: Reverted to cos^2 alpha = 1.d-16 when negative.
C 23 May 2005: Reverted to cos^2 alpha = 1.d-38 when negative.
C 23 May 2005: Reverted to old alternate with errors
C 19 May 2005: included both old and new formulas
C 18 May 2005: Corrected formulas for PVTHTH and PVPHTH
C 18 May 2005: Corrected formulas for PVTHPH and PVPHPH
C 17 May 2005: Added and changed some debugging statements.
C 16 May 2005: Added some debugging statements.
C 14 May 2005: Corrected formula for PVHR
C 12 May 2005: Reorganized to avoid dividing by zero.
C 11 May 2005: Added print statements at end
C 3 May 2005: Declared LAMBDA0 double precision.
C-----------------------------------------------------------------------------
      SUBROUTINE VDRAFT3  
C-----------------------------------------------------------------------------
C     WIND VELOCITY PERTURBATION MODEL 
C     An updraft plus a downdraft centered at a specified longitude and 
C     latitude, based on a stream function that is capped at the top and 
C     bottom.  
      implicit double precision (a-h,o-z)

      REAL(KIND=8) :: h, z1, delta1, z2, delta2, c1, c2
      REAL(KIND=8) :: th0, gam, singam ! 31 Jul 2006

      TYPE :: value_and_derivs_1_dim
         REAL(KIND=8) :: value
         REAL(KIND=8) :: p,pp
      END TYPE 

      TYPE (value_and_derivs_1_dim) :: psiA1, psiB1sg, psiC1 ! 7 Aug 2006
      TYPE (value_and_derivs_1_dim) :: psiA2, psiB2sg, psiC2 ! 7 Aug 2006
      TYPE (value_and_derivs_1_dim) :: psiB1
      TYPE (value_and_derivs_1_dim) :: psiB2
      TYPE (value_and_derivs_1_dim) :: psiD1sg, psiE1
      TYPE (value_and_derivs_1_dim) :: psiD2sg, psiE2
      TYPE (value_and_derivs_1_dim) :: psiD1
      TYPE (value_and_derivs_1_dim) :: psiD2

      TYPE :: value_and_derivs_2_dim
         REAL(KIND=8) :: value
         REAL(KIND=8) :: h,g
         REAL(KIND=8) :: hh,hg,gg
C                        ,g means d/d cos gamma
      END TYPE 

      TYPE (value_and_derivs_2_dim) :: psi1sg, psi2sg, psisg ! 7 Aug 2006

C------------------------------------------------------------------------------

      TYPE :: model
         CHARACTER :: model_name
         REAL(KIND=8) :: data_id
         CHARACTER :: pert_name
         REAL(KIND=8) :: pert_id
      END TYPE model
      TYPE (model) :: modug

C------------------------------------------------------------------------------

C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA       
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10

C     COMMON DECK "RKAM" INSERTED HERE        
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,
     1 DRDT,DTHDT,DPHDT,DKRDT(17)  
      double precision kr,kth,kph                                  ! 10 Mar 2008

C     COMMON DECK "UU" INSERTED HERE     
      DOUBLE PRECISION MODU, LAMBDA0 ! 3 May 2005
      COMMON/UU/MODU(4)                   
     1  ,V0  ,PVT0  ,PVR0  ,PVTH0  ,PVPH0      
     2  ,VR0 ,PVRT0 ,PVRR0 ,PVRTH0 ,PVRPH0     
     3  ,VTH0,PVTHT0,PVTHR0,PVTHTH0,PVTHPH0    
     4  ,VPH0,PVPHT0,PVPHR0,PVPHTH0,PVPHPH0    

      COMMON/WW/ dum(10), MAXW,W(1000)
      EQUIVALENCE (EARTHR,W(1))
     1,(AZ1,W(10))
C     WIND PERTURBATION   125-149      
      EQUIVALENCE (LAMBDA0,W(128)),(PHI0,W(129))
      EQUIVALENCE (W1,W(130)),(GAMMA1,W(131)),(H1,W(132))    
      EQUIVALENCE (W2,W(133)),(GAMMA2,W(134)),(H2,W(135))    
      EQUIVALENCE (z1,W(136)),(delta1,W(137))
      EQUIVALENCE (z2,W(138)),(delta2,W(139))
      EQUIVALENCE (stren1,W(140)),(thick1,W(141))
      EQUIVALENCE (deltaz1,W(142)),(deltag1,W(143))
      EQUIVALENCE (stren2,W(144)),(thick2,W(145))
      EQUIVALENCE (deltaz2,W(146)),(deltag2,W(147))

C     COMMON DECK "B2" INSERTED HERE                        
      INTEGER DUMX,DUNTBL,DUITBL,DUFRMTBL,IDSDU(20)              
      COMMON/B2/DUMX,DUNTBL(10),DUITBL(10),DUFRMTBL(10),DUGP(10) 
      EQUIVALENCE (DUGP,IDSDU)                                
      INTEGER DUPX    ,DUQTBL(10),DULTBL(10),DUIRMTBL(10)       
      DATA DUPX/1/                                           
      DATA DUQTBL/1,11,8*0/    
      DATA DULTBL/1,9*0/       
      DATA DUIRMTBL/1,9*0/
C                             
C-----------------------------------------------------------------------------
      ENTRY SETPWN            
C-----------------------------------------------------------------------------
c      print*,'Beginning of setpwn'
      DUMX=DUPX
      CALL IMOVE(DUNTBL,DUQTBL,10)     
      CALL IMOVE(DUITBL,DULTBL,10)     
      CALL IMOVE(DUFRMTBL,DUIRMTBL,10) 
c      print*,'before calcmax'
      CALL CALCMAX
      CALL CALCMAX2
      RETURN
C                                    
C-----------------------------------------------------------------------------
      ENTRY IPWINDR
C-----------------------------------------------------------------------------
      IF(W(125).NE.4.) THEN   
C      IF(IDINT(W(125)).NE.4) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(125) should be 4 if you are using model VDRAFT3.'
      ENDIF

      MODU(3)=8HVDRAFT3  
      MODU(4)=W(127)
      EARSQ=EARTHR**2
      SING1SQ=(DSIN(GAMMA1))**2 ! 27 May 2005
      SING2SQ=(DSIN(GAMMA2))**2 ! 27 May 2005
      th0 = pid2 - Lambda0 ! 31 Jul 2006
      SINLAM0=DSIN(LAMBDA0)
      COSTH0=SINLAM0 ! 25 May 2005:
      COSLAM0=DCOS(LAMBDA0)
      SINTH0=COSLAM0 ! 25 May 2005:
      H1SQ=H1**2
      H2SQ=H2**2
      CONST1 = earthr*w1/2.d0 ! 8 Aug 2006
      CONST2 = earthr*w2/2.d0 ! 8 Aug 2006
      wv1=stren1*w1
      wv2=stren2*w2
      deltav1=thick1*delta1
      deltav2=thick2*delta2
      R1=EARTHR+Z1+DELTAZ1
      R2=EARTHR+Z2+DELTAZ2
      E1=wv1*r1/6.d0
      E2=wv2*r2/6.d0
      F1=deltav1/r1
      F2=deltav2/r2
      RETURN         
                    
C-----------------------------------------------------------------------------
      ENTRY IPWINDg(modug)                                    !  6 Jan 2011
C-----------------------------------------------------------------------------
      IF(W(125).NE.4.) THEN   
C      IF(IDINT(W(125)).NE.4) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(125) should be 4 if you are using model VDRAFT3.'
      ENDIF

      MODU(3)='VDRAFT3'  
      MODU(4)=W(127)
      modug%pert_name = 'VDRAFT3'
      modug%pert_id = W(127)
      EARSQ=EARTHR**2
      SING1SQ=(DSIN(GAMMA1))**2 ! 27 May 2005
      SING2SQ=(DSIN(GAMMA2))**2 ! 27 May 2005
      th0 = pid2 - Lambda0 ! 31 Jul 2006
      SINLAM0=DSIN(LAMBDA0)
      COSTH0=SINLAM0 ! 25 May 2005:
      COSLAM0=DCOS(LAMBDA0)
      SINTH0=COSLAM0 ! 25 May 2005:
      H1SQ=H1**2
      H2SQ=H2**2
      CONST1 = earthr*w1/2.d0 ! 8 Aug 2006
      CONST2 = earthr*w2/2.d0 ! 8 Aug 2006
      wv1=stren1*w1
      wv2=stren2*w2
      deltav1=thick1*delta1
      deltav2=thick2*delta2
      R1=EARTHR+Z1+DELTAZ1
      R2=EARTHR+Z2+DELTAZ2
      E1=wv1*r1/6.d0
      E2=wv2*r2/6.d0
      F1=deltav1/r1
      F2=deltav2/r2
      RETURN         
                    
C-----------------------------------------------------------------------------
      ENTRY PWINDR
C-----------------------------------------------------------------------------
      RSQ=R**2
      H=R-EARTHR
      HSQ=H**2
      HCUBE=H*HSQ
      COSTH=DCOS(TH)
      SINTH=DSIN(TH)
      COSPH=DCOS(PH-PHI0)
      SINPH=DSIN(PH-PHI0)
C      COSGAM=SINLAM0*COSTH+COSLAM0*SINTH*COSPH
      call bettercea(th0,th,gam,costh0,costh,cosgam,
     1 sinth0,sinth,singam,ph-phi0) !  2 Aug 2006

      CALL psiA(h1, h, psiA1)
      CALL psiB(CONST1, sing1sq, cosgam, singam, psiB1sg) ! 31 Jul 2006
      CALL psiC(z1, delta1, h, psiC1)
      CALL psiD(E1, gamma1-deltag1, F1, gam, cosgam, singam, psiD1sg)
      CALL psiE(z1+deltaz1, deltav1, h, psiE1)

      CALL psiA(h2, h, psiA2)
      CALL psiB(CONST2, sing2sq, cosgam, singam, psiB2sg) ! 31 Jul 2006
      CALL psiC(z2, delta2, h, psiC2)
      CALL psiD(E2, gamma2-deltag2, F2, gam, cosgam, singam, psiD2sg)
      CALL psiE(z2+deltaz2, deltav2, h, psiE2)

      psi1sg%value = psiA1%value*psiB1sg%value*psiC1%value 
     1             + psiD1sg%value*psiE1%value
      psi1sg%h = psiB1sg%value*(psiA1%p*psiC1%value+psiA1%value*psiC1%p)
     1         + psiD1sg%value*psiE1%p
      psi1sg%g = psiA1%value*psiB1sg%p*psiC1%value
     1         + psiD1sg%p*psiE1%value
c      print *, ' psiA1 value = ', psiA1%value
c      print *, ' psiA1 p = ', psiA1%p
c      print *, ' psiB1sg value = ', psiB1sg%value
c      print *, ' psiB1sg p = ', psiB1sg%p
c      print *, ' psiC1 value = ', psiC1%value
c      print *, ' psiC1 p = ', psiC1%p
c      print *, ' psiD1sg value = ', psiD1sg%value
c      print *, ' psiD1sg p = ', psiD1sg%p
c      print *, ' psiE1 value = ', psiE1%value
c      print *, ' psiE1 p = ', psiE1%p
c      print *, ' psi1sg g = ', psi1sg%g
      psi1sg%hh = psiB1sg%value*(psiA1%pp*psiC1%value + 
     1       2.d0*psiA1%p*psiC1%p + psiA1%value*psiC1%pp)
     2          + psiD1sg%value*psiE1%pp
c      psi1sg%gg = psiA1%value*psiB1sg%pp*psiC1%value
      psi1sg%hg = psiB1sg%p*(psiA1%p*psiC1%value+psiA1%value*psiC1%p)
     1          + psiD1sg%p*psiE1%p

      psi2sg%value = psiA2%value*psiB2sg%value*psiC2%value 
     1             + psiD2sg%value*psiE2%value
      psi2sg%h = psiB2sg%value*(psiA2%p*psiC2%value+psiA2%value*psiC2%p)
     1         + psiD2sg%value*psiE2%p
      psi2sg%g = psiA2%value*psiB2sg%p*psiC2%value
     1         + psiD2sg%p*psiE2%value
      psi2sg%hh = psiB2sg%value*(psiA2%pp*psiC2%value + 
     1       2.d0*psiA2%p*psiC2%p + psiA2%value*psiC2%pp)
     2          + psiD2sg%value*psiE2%pp
c      psi2sg%gg = psiA2%value*psiB2sg%pp*psiC2%value
      psi2sg%hg = psiB2sg%p*(psiA2%p*psiC2%value+psiA2%value*psiC2%p)
     1          + psiD2sg%p*psiE2%p

      SINGSQ=SINGAM**2 ! 1 Aug 2006

      psisg%value = psi1sg%value + psi2sg%value
      psisg%h = psi1sg%h + psi2sg%h
      psisg%g = psi1sg%g + psi2sg%g
      psisg%hh = psi1sg%hh + psi2sg%hh
      psisg%hg = psi1sg%hg + psi2sg%hg
C      psisg%gg = psi1sg%gg + psi2sg%gg

      VR = - psisg%g/r ! 7 Aug 2006
c      print *, ' psisg g = ', psisg%g
c      print *, ' vr = ', vr

      PVRR = - psisg%hg/r + psisg%g/(rsq) ! 7 Aug 2006

      SGPGTH=SINLAM0*SINTH-COSLAM0*COSTH*COSPH
      SGPGPH=COSLAM0*SINTH*SINPH

      IF(SINGAM.NE.0.D0) THEN

         PGTH=SGPGTH/SINGAM 
         PGPH=SGPGPH/SINGAM 

         VH = (-psisg%h - psisg%value/r)/(singam)

         SINALP=COSLAM0*SINPH/SINGAM
         SINALPSQ=SINALP**2
         COSALP=(SINTH*COSTH0-SINTH0*COSTH*COSPH)/SINGAM ! 25 May 2005
c      cosalp = (COSTH0-COSTH*COSGAM)/(SINTH*SINGAM) ! equivalent
         COSALPSQ=COSALP*COSALP ! 25 May 2005
c         one=sinalpsq+cosalpsq ! 26 May 2005
         VTH=VH*COSALP
         VPH=VH*SINALP

         PVHR = -psisg%hh/(singam) - psisg%h/(r*singam) + 
     1   psisg%value/(rsq*singam) ! 7 Aug 2006

         PVTHR=PVHR*COSALP
         COS2ALP=1.D0-2.D0*SINALPSQ ! 17 May 2005
         PSINALPTH=-SINALP*COSGAM/SINGSQ*SGPGTH ! 18 May 2005
         PCOSALPTH=(COSTH*COSTH0+SINTH*SINTH0*COSPH)/SINGAM ! 25 May 2005
     1    -COSALP*COSGAM*SGPGTH/SINGSQ  ! 26 May 2005
         PSINALPPH=COSLAM0*COSPH/SINGAM-SINALP*COSGAM/SINGSQ*SGPGPH !18 May 2005
         PCOSALPPH=SINTH0*COSTH*SINPH/SINGAM
     1   -COSALP*COSGAM*SGPGPH/SINGSQ ! 26 May 2005

         psiB1%value = psiB1sg%value/singam
         psiB2%value = psiB2sg%value/singam
         psiB1%p = psiB1%value*cosgam/singsq + psiB1sg%p/singam
         psiB2%p = psiB2%value*cosgam/singsq + psiB2sg%p/singam
         psiD1%value = psiD1sg%value/singam
         psiD2%value = psiD2sg%value/singam
         psiD1%p = psiD1%value*cosgam/singsq + psiD1sg%p/singam
         psiD2%p = psiD2%value*cosgam/singsq + psiD2sg%p/singam

C The following 3 lines have a different sign from VDRAFT because
C the definition of PVHPG is different here.
         PVHPG = (-psisg%hg -psisg%g/r)/(singam) - 
     1   (psisg%h*cosgam + psisg%value/r)/(singam*singsq) ! 7 Aug 2006
         PVHTH=-PVHPG*SGPGTH 
         PVHPH=-PVHPG*SGPGPH

         PVTHTH=PVHTH*COSALP+VH*PCOSALPTH ! 18 May 2005
         PVTHPH=PVHPH*COSALP+VH*PCOSALPPH ! 18 May 2005
         PVPHR=PVHR*SINALP
         PVPHTH=PVHTH*SINALP+VH*PSINALPTH ! 18 May 2005
         PVPHPH=PVHPH*SINALP+VH*PSINALPPH ! 18 May 2005

      ELSE ! 3 Aug 2006
         DEN=DSQRT(DTHDT**2+(DPHDT*SINTH)**2)
         IF(DEN.NE.0.D0) THEN
            COSALP=DTHDT*DSIGN(1.D0,DPHDT)/DEN
            SINALP=SINTH*DPHDT*DSIGN(1.D0,DTHDT)/DEN
         ELSE
            COSALP=DCOS(AZ1)
            SINALP=-DSIN(AZ1)
         ENDIF
         PGTH=COSALP
         PGPH=SINTH*SINALP

         VTH=0.d0
         VPH=0.D0
         PVTHR=0.d0
         PVPHR=0.d0
         SGPVHPG = 
     1   CONST1*(psiA1%p*psiC1%value+psiA1%value*psiC1%p + 
     2   psiA1%value* psiC1%value/r) +
     3   CONST2*(psiA2%p*psiC2%value+psiA2%value*psiC2%p + 
     4   psiA2%value* psiC2%value/r)
      print *, ' SGPVHPG = ', SGPVHPG

c         PVPHTH=0.D0
c         PVTHPH=0.D0
C For South of axis, opposite sign for North
c         PVTHTH=-SGPVHPG
C For East of axis, opposite sign for West
c         PVPHPH=-SINTH*SGPVHPG
C in the following formula, the derivatives psiB1sg%p and psiB2sg%p are used
C only to get the value of the function divided by sin gamma.
         VHDSG=
     1      0.5D0*psiB1sg%p*(
     2 psiA1%p*psiC1%value+psiA1%value*psiC1%p+psiA1%value*psiC1%value/r
     3   ) + 0.5D0*psiB2sg%p*(
     4psiA2%p*psiC2%value+psiA2%value*psiC2%p+psiA2%value*psiC2%value/r)
      print *, ' VHDSG = ', VHDSG

         PVHTH=-SGPVHPG*PGTH 
         PVHPH=-SGPVHPG*PGPH

         PVTHTH=PVHTH*COSALP+VHDSG*(COSTH*COSTH0+SINTH*SINTH0*COSPH)
         PVTHPH=PVHPH*COSALP+VHDSG*SINTH0*COSTH*SINPH
         PVPHR=0.d0
         PVPHTH=PVHTH*SINALP+VHDSG*(-SINALP*PGTH)
         PVPHPH=PVHPH*SINALP+VHDSG*(COSLAM0*COSPH-SINALP*PGPH)

      ENDIF

      PVRTH = (psiA1%value*psiB1sg%pp*psiC1%value*sgpgth 
     1       - psiD1sg%pp*psiE1%value*pgth
     2       + psiA2%value*psiB2sg%pp*psiC2%value*sgpgth 
     3       - psiD2sg%pp*psiE2%value*pgth)/r
      PVRPH = (psiA1%value*psiB1sg%pp*psiC1%value*sgpgph 
     1       - psiD1sg%pp*psiE1%value*pgph
     2       + psiA2%value*psiB2sg%pp*psiC2%value*sgpgph 
     3       - psiD2sg%pp*psiE2%value*pgph)/r

      PVRT=0.d0
      PVTHT=0.d0
      PVPHT=0.d0

      VR0 = VR0 + VR
      PVRT0 = PVRT0 + PVRT
      PVRR0 = PVRR0 + PVRR
      PVRTH0 = PVRTH0 + PVRTH 
      PVRPH0 = PVRPH0 + PVRPH
      VTH0 = VTH0 + VTH
      PVTHT0 = PVTHT0 + PVTHT
      PVTHR0 = PVTHR0 + PVTHR
      PVTHTH0 = PVTHTH0 + PVTHTH
      PVTHPH0 = PVTHPH0 + PVTHPH
      VPH0 = VPH0 + VPH
      PVPHT0 = PVPHT0 + PVPHT
      PVPHR0 = PVPHR0 + PVPHR
      PVPHTH0 = PVPHTH0 + PVPHTH
      PVPHPH0 = PVPHPH0 + PVPHPH
C
c      print *, ' vr0 = ', vr0
c      print *, ' vth0 = ', vth0
c      print *, ' vph0 = ', vph0
      VSQ0 = VR0**2 + VTH0**2 + VPH0**2
c      IF(VSQ0.LT.0.D0) THEN
c         PRINT*, 'Velocity cannot be zero!'
c         STOP
c      ENDIF

      V0 = DSQRT(VSQ0)
      PVT0 = 0.d0
      PVR0 = 0.d0
      PVTH0 = 0.d0
      PVPH0 = 0.d0
      IF (V0.EQ.0.d0) RETURN
      PVT0 = (VR0*PVRT0 + VTH0*PVTHT0 + VPH0*PVPHT0)/V0
      PVR0 = (VR0*PVRR0 + VTH0*PVTHR0 + VPH0*PVPHR0)/V0
      PVTH0 = (VR0*PVRTH0 + VTH0*PVTHTH0 + VPH0*PVPHTH0)/V0
      PVPH0 = (VR0*PVRPH0 + VTH0*PVTHPH0 + VPH0*PVPHPH0)/V0

c      print *, ' pvr0 = ', pvr0
c      print *, ' pvth0 = ', pvth0
c      print *, ' pvph0 = ', pvph0
C
      RETURN
C-----------------------------------------------------------------------------
      END SUBROUTINE VDRAFT3
C-----------------------------------------------------------------------------
      SUBROUTINE psiA (hi,h,psi)
C-----------------------------------------------------------------------------
      implicit double precision (A-H,O-Z)                        !  9 Nov 2009
      REAL(KIND=8) :: h, hi, hsq, hisq, sum, sumsq, sumcube

      TYPE :: value_and_derivs_1_dim
         REAL(KIND=8) :: value
         REAL(KIND=8) :: p,pp
      END TYPE 

      TYPE (value_and_derivs_1_dim) :: psi

      hsq=h*h

      psi%value = 1.d0
      psi%p = 0.d0
      psi%pp = 0.d0

      IF (hi.EQ.0.d0) RETURN

      hisq = hi*hi
      sum = hisq + hsq
      sumsq = sum*sum

      psi%value = hsq/sum
      psi%p = 2.d0*hisq*h/sumsq
      sumcube = sumsq*sum
      psi%pp = 2.d0*hisq*(hisq-3.d0*hsq)/sumcube 

      RETURN
C-----------------------------------------------------------------------------
      END SUBROUTINE psiA
C-----------------------------------------------------------------------------
      SUBROUTINE psiB (const, singisq, cosgam, singam, psisg)
C-----------------------------------------------------------------------------
      implicit double precision (A-H,O-Z)
C 6 Nov 2009      common/crmach/rmach(9) ! 31 Jul 2006
C 6 Nov 2009      REAL(KIND=8) :: rmach, right ! 31 Jul 2006
C 6 Nov 2009      equivalence (rmach(3),right) ! equals the smallest relative spacing
      double precision d1mach                                      ! 6 Nov 2009
      double precision right                                       ! 6 Nov 2009
      REAL(KIND=8) :: const, singisq, cosgam, singam, cosgsq, exbet

      TYPE :: value_and_derivs_1_dim
         REAL(KIND=8) :: value
         REAL(KIND=8) :: p,pp
      END TYPE 

      TYPE (value_and_derivs_1_dim) :: psisg, beta

c      one = cosgam**2 +singam**2
c      print *, ' one  = ', one

      right=d1mach(3)                                              ! 6 Nov 2009
c      print*,'right = ',right

      psisg%value = const*singisq
      psisg%p = 0.d0 ! = d psisg/ d cos gamma
      psisg%pp = 0.d0 ! = d2 psisg/ d cos gamma 2

      IF (singisq.EQ.0.d0) RETURN

c      cosgsq = cosgam*cosgam

c      beta%value=(1.d0 - cosgsq)/singisq
      beta%value=singam**2/singisq ! 31 Jul 2006
      exbet =DEXP(-beta%value)

      beta%p = - 2.d0*cosgam/singisq ! = d beta / d cos gamma
      beta%pp = - 2.d0/singisq ! = d 2 beta / d cos gamma 2

      psisg%p = const*singisq*exbet*beta%p
      psisg%pp = const*singisq*exbet*(beta%pp-beta%p**2)

c      print *, ' right  = ', right
      if(beta%value.eq.0.D0) THEN
      psisg%value = 0.D0
      error = 0.d0
      error2 = 0.d0
c      print *, ' beta  = ', beta%value
c      print *, ' psisg   = ' , psisg%value
c      print *, ' error = ' , error, error2, 0
      return
      ENDIF

      IF(beta%value/2.D0 .LT. right) THEN
         psisg%value = const*singisq*beta%value
         error = const*singisq*beta%value**2/2.d0
         error2 = error/psisg%value
c      print *, ' beta  = ', beta%value
c      print *, ' psisg   = ' , psisg%value
c      print *, ' error = ' , error, error2, 1
c      print *, ' error = ' , beta%value, error2, 1
         return

      else IF (beta%value**2 /6.D0 .LT. right) THEN
         psisg%value = const*singisq*beta%value*
     2    (1.D0-BETA%VALUE/2.d0)
         error = const*singisq*beta%value**3/6.d0
         error2 = error/psisg%value
c      print *, ' beta  = ', beta%value
c      print *, ' psisg   = ' , psisg%value
c      print *, ' error = ' , error, error2, 2
c      print *, ' error = ' , beta%value, error2, 2
         return

      else IF (beta%value**3 /24.D0 .LT. right) THEN
         psisg%value = const*singisq*beta%value*
     2    (1.D0-BETA%VALUE/2.d0*
     3    (1.D0-BETA%VALUE/3.d0))
         error = const*singisq*beta%value**4/24.d0
         error2 = error/psisg%value
c      print *, ' beta  = ', beta%value
c      print *, ' psisg  = ' , psisg%value
c      print *, ' error = ' , error, error2, 3
         return

      else IF (beta%value**4 /120.D0 .LT. right) THEN
         psisg%value = const*singisq*beta%value*
     2    (1.D0-BETA%VALUE/2.d0*
     3    (1.D0-BETA%VALUE/3.d0*
     4    (1.D0-BETA%VALUE/4.d0)))
         error = const*singisq*beta%value**5/120.d0
         error2 = error/psisg%value
c      print *, ' beta  = ', beta%value
c      print *, ' psisg   = ' , psisg%value
c      print *, ' error = ' , error, error2, 4
         return

      else IF (beta%value**5 /720.D0 .LT. right) THEN
         psisg%value = const*singisq*beta%value*
     2    (1.D0-BETA%VALUE/2.d0*
     3    (1.D0-BETA%VALUE/3.d0*
     4    (1.D0-BETA%VALUE/4.d0*
     5    (1.D0-BETA%VALUE/5.d0))))
         error = const*singisq*beta%value**6/720.d0
         error2 = error/psisg%value
c      print *, ' beta  = ', beta%value
c      print *, ' psisg   = ' , psisg%value
c      print *, ' error = ' , error, error2, 5
         return

      else IF (beta%value**6 /5040.D0 .LT. right) THEN
         psisg%value = const*singisq*beta%value*
     2    (1.D0-BETA%VALUE/2.d0*
     3    (1.D0-BETA%VALUE/3.d0*
     4    (1.D0-BETA%VALUE/4.d0*
     5    (1.D0-BETA%VALUE/5.d0*
     6    (1.D0-BETA%VALUE/6.d0)))))
         error = const*singisq*beta%value**7/5040.d0
         error2 = error/psisg%value
c      print *, ' beta  = ', beta%value
c      print *, ' psisg   = ' , psisg%value
c      print *, ' error = ' , error, error2, 6
         return

      else IF (beta%value**7 /40320.D0 .LT. right) THEN
         psisg%value = const*singisq*beta%value*
     2    (1.D0-BETA%VALUE/2.d0*
     3    (1.D0-BETA%VALUE/3.d0*
     4    (1.D0-BETA%VALUE/4.d0*
     5    (1.D0-BETA%VALUE/5.d0*
     6    (1.D0-BETA%VALUE/6.d0*
     7    (1.D0-BETA%VALUE/7.d0))))))
         error = const*singisq*beta%value**8/40320.d0
         error2 = error/psisg%value
c      print *, ' beta  = ', beta%value
c      print *, ' psisg   = ' , psisg%value
c      print *, ' error = ' , error, error2, 7
         return

      else IF (beta%value**8 /362880.D0 .LT. right) THEN
         psisg%value = const*singisq*beta%value*
     2    (1.D0-BETA%VALUE/2.d0*
     3    (1.D0-BETA%VALUE/3.d0*
     4    (1.D0-BETA%VALUE/4.d0*
     5    (1.D0-BETA%VALUE/5.d0*
     6    (1.D0-BETA%VALUE/6.d0*
     7    (1.D0-BETA%VALUE/7.d0*
     8    (1.D0-BETA%VALUE/8.d0)))))))
         error = const*singisq*beta%value**9/362880.d0
         error2 = error/psisg%value
c      print *, ' beta  = ', beta%value
c      print *, ' psisg   = ' , psisg%value
c      print *, ' error = ' , error, error2, 8
         return

      else IF (beta%value**9 /3628800.D0 .LT. right) THEN
         psisg%value = const*singisq*beta%value*
     2    (1.D0-BETA%VALUE/2.d0*
     3    (1.D0-BETA%VALUE/3.d0*
     4    (1.D0-BETA%VALUE/4.d0*
     5    (1.D0-BETA%VALUE/5.d0*
     6    (1.D0-BETA%VALUE/6.d0*
     7    (1.D0-BETA%VALUE/7.d0*
     8    (1.D0-BETA%VALUE/8.d0*
     9    (1.D0-BETA%VALUE/9.d0))))))))
         error = const*singisq*beta%value**10/3628800.d0
         error2 = error/psisg%value
c      print *, ' beta  = ', beta%value
c      print *, ' psisg   = ' , psisg%value
c      print *, ' error = ' , error, error2, 9
         return

      ELSE
         psisg%value = const*singisq*(1.d0-exbet)
         error = right
         error2 = error/psisg%value
c      print *, ' beta  = ', beta%value
c      print *, ' psisg   = ' , psisg%value
c      print *, ' error = ' , error, error2, 20
         RETURN
      endif

C-----------------------------------------------------------------------------
      END SUBROUTINE psiB
C-----------------------------------------------------------------------------
      SUBROUTINE psiC (h0, delta, h, psi)
C-----------------------------------------------------------------------------
      implicit double precision (A-H,O-Z)                        !  9 Nov 2009
C 6 Nov 2009      common/crmach/rmach(9) ! 6 Sept 2006
C 6 Nov 2009      REAL(KIND=8) :: rmach, large ! 6 Sept 2006
C 6 Nov 2009      equivalence (rmach(2),large) ! equals the largest magnitude
      double precision d1mach                                      ! 6 Nov 2009
      double precision large                                       ! 6 Nov 2009
      REAL(KIND=8) :: h0, delta, h, arg, tanharg, secharg, sechargsq
      REAL(KIND=8) :: maxlarge

      TYPE :: value_and_derivs_1_dim
         REAL(KIND=8) :: value
         REAL(KIND=8) :: p,pp
      END TYPE 

      TYPE (value_and_derivs_1_dim) :: psi

      psi%value = 1.d0
      psi%p = 0.d0 ! = d psi/ d h
      psi%pp = 0.d0 ! = d2 psi/ d h 2
c      print *, ' psi value = ', psi%value

      IF (delta.LE.0.d0) RETURN

      arg = (h-h0)/delta
      tanharg = dtanh(arg)

      if (dabs(arg).ge.maxlarge) then              ! 6 Sept 2006
         secharg = 0.d0
      else
         secharg = 1.d0/dcosh(arg)  ! 13 Jul 2006
      endif   

c      if (arg.gt.1250) then
c         print*,'h = ',h
c         print*,'h0 = ',h0
c         print*,'delta = ',delta
c         print*,'arg 1= ',arg
c      endif
c      secharg = 1.d0/dcosh (arg)  ! 13 Jul 2006
c      if (arg.gt.1250) print*,'arg 2= ',arg

c      psi%value = (1.d0 - tanharg)/2.d0

c      print*,'maxlarge = ',maxlarge
      if (2.d0*arg.ge.maxlarge) then         ! 6 Sept 2006
         psi%value = 0.d0
      else
         psi%value = 1.d0/(dexp(2.d0*arg) + 1.d0)  ! 1 Aug 2006
      endif   
c      print *, ' psi value = ', psi%value

      sechargsq = secharg**2
      psi%p = - sechargsq/(2.d0*delta)
      psi%pp = sechargsq*tanharg/(delta**2)

      RETURN

      ENTRY CALCMAX
c      large = 2
c      do 237 i = 2,1500
c      large = 2*large
c      print*,'large = ',large
      large=d1mach(2)                                              ! 6 Nov 2009
      maxlarge = dlog(large)
c      print*,'large = ',large
c      print*,'maxlarge = ',maxlarge

c      stop
      RETURN

C-----------------------------------------------------------------------------
      END SUBROUTINE psiC
C-----------------------------------------------------------------------------
      SUBROUTINE psiD (D, CAPGAM, DELTA, gamma, cosgam, singam, psisg)
C-----------------------------------------------------------------------------
      implicit double precision (A-H,O-Z)

      TYPE :: value_and_derivs_1_dim
         REAL(KIND=8) :: value
         REAL(KIND=8) :: p,pp
      END TYPE 

      TYPE (value_and_derivs_1_dim), intent (out) :: psisg
      double precision, intent (in)::D,CAPGAM,DELTA,gamma,cosgam,singam

      double precision d1mach 
      double precision right  ! equals the smallest relative spacing

c      one = cosgam**2 +singam**2
c      print *, ' one  = ', one

      right=d1mach(3) ! equals the smallest relative spacing
c      print*,'right = ',right 

      psisg%value = 0.d0 
      psisg%p = 0.d0 ! = d psisg/ d cos gamma
      psisg%pp = 0.d0 ! = d/d gamma d psisg/ d cos gamma 

c      print*,'D = ',D 
c      print*,'DELTA = ',DELTA
      IF(D.EQ.0.D0.OR.DELTA.EQ.0.D0) RETURN

c      print*,' singam = ', singam
      singsq = singam*singam
      cosgsq = cosgam*cosgam
      DELTAsq = DELTA*DELTA
C      print*,' CAPGAM = ', CAPGAM
C      print*,' GAMMA = ', GAMMA
      arg = (gamma-CAPGAM)/DELTA
      temp = D/DELTA *dexp(-arg*arg) 
      psisg%value = temp *singam**3
c      print*,'ARG = ',ARG 
c      print*,'temp = ',temp
c      print*,' psisg%value = ', psisg%value
      psisg%p = -temp * (-2.d0*arg/DELTA*singam +3.d0*cosgam)*singam
      psisg%pp = -temp * (4.d0*arg*arg/DELTAsq *singsq 
     1 -10.d0*arg/DELTA*singam*cosgam 
     2 -2.d0/DELTAsq *singsq 
     3 +3.d0*cosgsq -3.d0*singsq) ! = d/d gamma d psisg/ d cos gamma 

      RETURN

C-----------------------------------------------------------------------------
      END SUBROUTINE psiD
C-----------------------------------------------------------------------------
      SUBROUTINE psiE (h0, delta, h, psi)
C-----------------------------------------------------------------------------
      implicit double precision (A-H,O-Z)                        !  9 Nov 2009

      TYPE :: value_and_derivs_1_dim
         double precision value
         double precision p,pp
      END TYPE 

      TYPE (value_and_derivs_1_dim), intent (out) :: psi
      double precision, intent(in):: h0, delta, h

      double precision d1mach 
      double precision large ! equals the largest magnitude
      double precision maxlarge

      psi%value = 0.d0
      psi%p = 0.d0 ! = d psi/ d h
      psi%pp = 0.d0 ! = d2 psi/ d h 2
c      print *, ' psi value = ', psi%value
c      print *, ' delta = ', delta

      IF (delta.LE.0.d0) RETURN

      arg = (h-h0)/delta

      if (dabs(arg).ge.maxlarge) RETURN

      if (arg.gt.1250) then
C         print*,'h = ',h
C         print*,'h0 = ',h0
C         print*,'delta = ',delta
C         print*,'arg 1= ',arg
      endif

c      print*,'maxlarge = ',maxlarge
         psi%value = dexp(-arg*arg)
c      print *, ' psi value = ', psi%value

      psi%p = - 2.d0 *arg *psi%value/delta
      psi%pp = - 2.d0 *(1-2.d0*arg*arg) *psi%value/delta**2

      RETURN

      ENTRY CALCMAX2
c      large = 2
c      do 237 i = 2,1500
c      large = 2*large
c      print*,'large = ',large
      large=d1mach(2)                                              ! 6 Nov 2009
      maxlarge = dlog(large)
c      print*,'large = ',large
c      print*,'maxlarge = ',maxlarge

c      stop
      RETURN

C-----------------------------------------------------------------------------
      END SUBROUTINE psiE
C-----------------------------------------------------------------------------