C vdraft.for 
C-----------------------------------------------------------------------------
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 10 Mar 2008 - corrected length errors in blank common
C 24 Apr 2006: changed 1 to 1.d0 and 0 to 0.d0
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 VDRAFT  
C-----------------------------------------------------------------------------
C     WIND VELOCITY MODEL 
C     An updraft plus a downdraft centered at a specified longitude and 
C     latitude, based on a stream function. 
      implicit double precision (a-h,o-z)
      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))    
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                             
      ENTRY SETPWN            
      DUMX=DUPX
      CALL IMOVE(DUNTBL,DUQTBL,10)     
      CALL IMOVE(DUITBL,DULTBL,10)     
      CALL IMOVE(DUFRMTBL,DUIRMTBL,10) 
      RETURN                         
C                                    
      ENTRY IPWINDR
      IF(W(125).NE.1.) THEN   
C      IF(IDINT(W(125)).NE.1) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(125) should be 1 if you are using model VDRAFT.'
      ENDIF

      MODU(3)=8HVDRAFT  
      MODU(4)=W(127)
      EARSQ=EARTHR**2
      SING1SQ=(DSIN(GAMMA1))**2 ! 27 May 2005
      SING2SQ=(DSIN(GAMMA2))**2 ! 27 May 2005
      SINLAM0=DSIN(LAMBDA0)
      COSTH0=SINLAM0 ! 25 May 2005:
      COSLAM0=DCOS(LAMBDA0)
      SINTH0=COSLAM0 ! 25 May 2005:
      H1SQ=H1**2
      H2SQ=H2**2

C     CALL IPWINDR                                
      RETURN         
C                    
      ENTRY PWINDR
      V=0.d0 ! 24 Apr 2006 and below
      PVT=0.d0
      PVR=0.d0
      PVTH=0.d0
      PVPH=0.d0
      VR=0.d0
      PVRT=0.d0
      PVRR=0.d0
      PVRTH=0.d0
      PVRPH=0.d0
      VTH=0.d0
      PVTHT=0.d0
      PVTHR=0.d0
      PVTHTH=0.d0
      PVTHPH=0.d0
      VPH=0.d0
      PVPHT=0.d0
      PVPHR=0.d0
      PVPHTH=0.d0
      PVPHPH=0.d0
      VH=0.D0
      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)
      COSGAM=SINLAM0*COSTH+COSLAM0*SINTH*COSPH
      SINGSQ=1.d0-COSGAM**2 ! 24 Apr 2006
      IF(SINGSQ.LT.0.D0) THEN
         write (3,*) 'singsq', singsq
         PRINT*, 'singsq', singsq ! 16 May 2005
         stop ! 16 May 2005
      ENDIF ! 16 May 2005
      SINGAM=DSQRT(SINGSQ)
      BETA1=SINGSQ/SING1SQ
      BETA2=SINGSQ/SING2SQ
      EXBET1=DEXP(-BETA1)
      EXBET2=DEXP(-BETA2)
      SUM1=H1SQ+HSQ
      SUM2=H2SQ+HSQ
      SUM1SQ=SUM1**2
      SUM2SQ=SUM2**2
      TEMP=EARSQ*H*COSGAM/RSQ
      VR1=TEMP*(W1*EXBET1/SUM1)
      VR2=TEMP*(W2*EXBET2/SUM2)
      VR=(VR1+VR2)*H
      PVRR1=2.D0*VR1/R*(EARTHR*H1SQ-HCUBE)/SUM1
      PVRR2=2.D0*VR2/R*(EARTHR*H2SQ-HCUBE)/SUM2
      PVRR=PVRR1+PVRR2
      SECGAM=1.D0/COSGAM
      SGPGTH=SINLAM0*SINTH-COSLAM0*COSTH*COSPH
      SGPGPH=COSLAM0*SINTH*SINPH
      PVRG1=-VR1*H*(SECGAM+2.D0*COSGAM/SING1SQ)
      PVRG2=-VR2*H*(SECGAM+2.D0*COSGAM/SING2SQ)
      PVRTH=SGPGTH*(PVRG1+PVRG2)
      PVRPH=SGPGPH*(PVRG1+PVRG2)
      IF(SINGAM.NE.0.D0) THEN
         TEMPB=-EARSQ/(R*SINGAM)
         TEMP1=1.D0-EXBET1
         TEMP2=1.D0-EXBET2
         VH1=TEMPB*W1*H1SQ*SING1SQ/SUM1SQ*TEMP1
         VH2=TEMPB*W2*H2SQ*SING2SQ/SUM2SQ*TEMP2
         VH=(VH1+VH2)*H
         SINALP=COSLAM0*SINPH/SINGAM
         SINALPSQ=SINALP**2
         COSALP=(SINTH*COSTH0-SINTH0*COSTH*COSPH)/SINGAM ! 25 May 2005
         COSALPSQ=COSALP*COSALP ! 25 May 2005
         one=sinalpsq+cosalpsq ! 26 May 2005
         VTH=VH*COSALP
         VPH=VH*SINALP
         PVHR1=VH1/R*(EARTHR*H1SQ-3.D0*EARTHR*HSQ-4.D0*HCUBE)/SUM1 ! 27 May 2005
         PVHR2=VH2/R*(EARTHR*H2SQ-3.D0*EARTHR*HSQ-4.D0*HCUBE)/SUM2 ! 27 May 2005
         PVHR=PVHR1+PVHR2
         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
         PVHPG1=VH1*H*(-1.D0+2.D0*BETA1*EXBET1/TEMP1)*COSGAM/SINGSQ !18 May 2005
         PVHPG2=VH2*H*(-1.D0+2.D0*BETA2*EXBET2/TEMP2)*COSGAM/SINGSQ !18 May 2005
         PVHPG=PVHPG1+PVHPG2
         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
      ENDIF
C      VSQ=VR**2+VH**2
C      IF(VSQ.Lt.0.D0) stop
C      V=DSQRT(VSQ)
C      IF(V.NE.0.D0) THEN
C         PVTH=(VR*PVRTH + VTH*PVTHTH + VPH*PVPHTH)/V   
C         PVPH=(VR*PVRPH + VTH*PVTHPH + VPH*PVPHPH)/V   
C         PVR=(VR*PVRR + VTH*PVTHR + VPH*PVPHR)/V
C      ENDIF
C     CALL PWINDR
C
      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
      VSQ0 = VR0**2 + VTH0**2 + VPH0**2
      IF(VSQ0.LT.0.D0) THEN
         PRINT*, 'Velocity cannot be zero!'
         STOP
      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
      RETURN
      END SUBROUTINE VDRAFT
