C vdraft2.for 
C-----------------------------------------------------------------------------
C 16 Dec 2009 - Correction
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 VDRAFT2  
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_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

      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
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)  
      double precision kr,kth,kph                                  ! 10 Mar 2008
      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))
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))
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
      RETURN
C                                    
C-----------------------------------------------------------------------------
      ENTRY IPWINDR
C-----------------------------------------------------------------------------
      IF(W(125).NE.2.) THEN   
C      IF(IDINT(W(125)).NE.2) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(125) should be 2 if you are using model VDRAFT2.'
      ENDIF

      MODU(3)=8HVDRAFT2  
      MODU(4)=W(127)
      EARSQ=EARTHR**2
      SING1SQ=(DSIN(GAMMA1))**2 ! 27 May 2005
      SING2SQ=(DSIN(GAMMA2))**2 ! 27 May 2005
C      c1 = earsq*w1*sing1sq/2.d0
C      c2 = earsq*w2*sing2sq/2.d0
      c1 = 0.d0                                                    ! 19 Mar 2008
      c2 = 0.d0                                                    ! 19 Mar 2008
      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
      RETURN         
C                    
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 psiA(h2, h, psiA2)
      CALL psiB(CONST2, sing2sq, cosgam, singam, psiB2sg) ! 31 Jul 2006
      CALL psiC(z2, delta2, h, psiC2)

      psi1sg%value = psiA1%value*psiB1sg%value*psiC1%value - c1
      psi1sg%h = psiB1sg%value*(psiA1%p*psiC1%value+psiA1%value*psiC1%p)
      psi1sg%g = psiA1%value*psiB1sg%p*psiC1%value
c      print *, ' psiA1 value = ', psiA1%value
c      print *, ' psiB1sg p = ', psiB1sg%p
c      print *, ' psiC1 value = ', psiC1%value
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)
      psi1sg%gg = psiA1%value*psiB1sg%pp*psiC1%value
      psi1sg%hg = psiB1sg%p*(psiA1%p*psiC1%value+psiA1%value*psiC1%p)

      psi2sg%value = psiA2%value*psiB2sg%value*psiC2%value - c2
      psi2sg%h = psiB2sg%value*(psiA2%p*psiC2%value+psiA2%value*psiC2%p)
      psi2sg%g = psiA2%value*psiB2sg%p*psiC2%value
      psi2sg%hh = psiB2sg%value*(psiA2%pp*psiC2%value + 
     1 2.d0*psiA2%p*psiC2%p + psiA2%value*psiC2%pp)
      psi2sg%gg = psiA2%value*psiB2sg%pp*psiC2%value
      psi2sg%hg = psiB2sg%p*(psiA2%p*psiC2%value+psiA2%value*psiC2%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
      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

      PVRTH = psisg%gg*sgpgth/r
      PVRPH = psisg%gg*sgpgph/r

      IF(SINGAM.NE.0.D0) THEN

         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

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
         VTH=0.d0
         VPH=0.D0
         PVTHR=0.d0
         PVPHR=0.d0
C For South of axis, opposite sign for North
         PVTHTH=-COSGAM*(
     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)) ! 16 Dec 2009
         PVPHTH=0.D0
         PVTHPH=0.D0
C For East of axis, opposite sign for West
         PVPHPH=SINTH*PVTHTH
      ENDIF
C   
      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
      END SUBROUTINE VDRAFT2
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
      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

      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

      END SUBROUTINE psiC