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

      SUBROUTINE VVORTX3  

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

C     WIND VELOCITY MODEL 

C     VERTICAL VORTEX WIND PERTURBATION WITH VISCOUS CORE AND

C     MULTIPLIES VELOCITY FIELD BY A GUASSIAN HEIGHT PROFILE.

      implicit double precision (a-h,o-z)

      COMMON/PCONST/CREF,RGAS,GAMMA       

      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10

      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)  

      DOUBLE PRECISION MODU

      COMMON/UU/MODU(4)                   

     1  ,V  ,PVT  ,PVR  ,PVTH  ,PVPH      

     2  ,VR ,PVRT ,PVRR ,PVRTH ,PVRPH     

     3  ,VTH,PVTHT,PVTHR,PVTHTH,PVTHPH    

     4  ,VPH,PVPHT,PVPHR,PVPHTH,PVPHPH    

      COMMON/WW/ dum(10), MAXW,W(1000)

      EQUIVALENCE (EARTHR,W(1))

C     WIND        100-124      

       EQUIVALENCE (U0,W(103)),(R0,W(104)),(TH0,W(105)),(PH0,W(106))    

       EQUIVALENCE (HWIDTH,W(107)),(HVMAX,W(108))           

C     COMMON DECK "B1" INSERTED HERE                        

      INTEGER UMX,UNTBL,UITBL,UFRMTBL,IDSU(20)              

      COMMON/B1/UMX,UNTBL(10),UITBL(10),UFRMTBL(10),UGP(10) 

      EQUIVALENCE (UGP,IDSU)                                

      INTEGER UPX    ,UQTBL(10),ULTBL(10),UIRMTBL(10)       

      DATA UPX/1/                                           

      DATA UQTBL/1,11,8*0/    

      DATA ULTBL/1,9*0/       

      DATA UIRMTBL/1,9*0/     

C                             

      ENTRY SETWND            

      UMX=UPX                 

      CALL IMOVE(UNTBL,UQTBL,10)     

      CALL IMOVE(UITBL,ULTBL,10)     

      CALL IMOVE(UFRMTBL,UIRMTBL,10) 

      CALL SETPWN                    

      RETURN                         

C                                    

      ENTRY IWINDR                   

      IF(W(100) .NE. 9.) THEN

           PRINT*,' Model check number mismatch: '

           STOP ' W(100) should be 9 if you are using model VVORTX3.'

      ENDIF



      MODU(1)=8HVVORTX3  

      MODU(2)= w(102)

       DENOM=0.d0        

       IF(HWIDTH.NE.0.d0) DENOM=1.d0/HWIDTH/HWIDTH

      CALL IPWINDR                                

      RETURN         

C                    

       ENTRY WINDR   

      CALL CLEAR(V,20)

      DR=R-EARTHR-HVMAX

       DTH = TH - (PID2-TH0)

       DPH = PH - PH0       

       RAD2 = EARTHR * dSQRT(DTH * DTH + DPH * DPH)  

       A = 1.397d0     

       B = - 1.26d0    

      EXPO=RAD2/R0     

      EXPO=B*EXPO*EXPO 

      EXB=0.d0         

       IF(EXPO .GT. -675.d0) EXB = dEXP(EXPO)   

      FX=1.d0-EXB                               

       RDUM = A * EARTHR * U0 * R0 / RAD2/rad2

      FZ=dEXP(-DR*DR*DENOM)                     

      DFDZ=-2.d0*DR*DENOM                       

      RDUM=FZ*RDUM

       DUX = FX / RAD2 + RAD2 * B * EXB / R0 * * 2

       VTH = - RDUM * FX * DPH

       VPH = RDUM * FX * DTH

      V=dSQRT(VTH*VTH + VPH*VPH)

       RDUM2=2.d0*RDUM*EARTHR*EARTHR

       PVTHTH = RDUM2 * DTH * DPH / RAD2 * DUX

       PVPHPH = - PVTHTH     

       PVTHPH = DPH*dph * RDUM2 / RAD2 * DUX -RDUM*FX

       PVPHTH = - DTH*dth* RDUM2 / RAD2* DUX + RDUM*FX

      PVTH=(VTH*PVTHTH + VPH*PVPHTH)/V              

      PVPH=(VTH*PVTHPH + VPH*PVPHPH)/V              

      PVTHR=VTH*DFDZ                                

      PVPHR=VPH*DFDZ                                

      PVR=(VTH*PVTHR+VPH*PVPHR)/V                   

      CALL PWINDR                                   

       RETURN                                       

       END                                          

                                                    



