C vvortx4.for 
C-----------------------------------------------------------------------------
! 27 Apr 2010 - Still need to add second derivatives
!  6 Apr 2010 - Moved the module to math-mod.for
!  5 Apr 2010 - Added a module for adding structures, element by element
!  2 Apr 2010 - Added structures for arguments 
!  2 Apr 2010 - Added entry point IPWINDg(modug)
!  2 Apr 2010 - Added entry point PWINDg(wind)
C  2 Jul 2009 - working correctly
C 30 Jun 2009 - copied from vdraft2.for
C-----------------------------------------------------------------------------
      SUBROUTINE VVORTX4  
C-----------------------------------------------------------------------------
C     WIND VELOCITY PERTURBATION MODEL 
C     A vortex with a viscous core and a Gaussian intensity profile in the
C     vertical.  The axis of the vortex is vertical and may be positioned
C     above any geographic latitude and longitude.  The vortex rotates 
C     anticlockwise looking down.  The core is essentially a solid-rotating
C     fluid while the outside falls off as the inverse distance from the center.
C-----------------------------------------------------------------------------

      USE v_and_g2_maths

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

      DOUBLE PRECISION :: h, z1, delta1, z2, delta2, c1, c2
      DOUBLE PRECISION :: th0, gam, singam 

      TYPE :: value_and_derivs_1_dim
         DOUBLE PRECISION :: value
         DOUBLE PRECISION :: p,pp
      END TYPE value_and_derivs_1_dim


      TYPE :: value_and_derivs_2_dim
         DOUBLE PRECISION :: value
         DOUBLE PRECISION :: h,g
         DOUBLE PRECISION :: hh,hg,gg
C                        ,g means d/d cos gamma
      END TYPE value_and_derivs_2_dim

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-----------------------------------------------------------------------------

      TYPE (vector_and_gradients_1_2) :: wind, wind1, wind2
      TYPE (gradient2) :: VHpp

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

C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA       
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3         
      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,DRDT(20)  
      double precision kr,kth,kph                     

C     COMMON DECK "UU" INSERTED HERE     
      DOUBLE PRECISION MODU, LAMBDA0 
      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 (u0,W(130)),(GAMMA0,W(131)),(Hmax,W(132))    
      EQUIVALENCE (WH,W(133))

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) 
      RETURN
C                                    
C-----------------------------------------------------------------------------
      ENTRY IPWINDR
C-----------------------------------------------------------------------------
      IF(W(125).NE.3.) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(125) should be 3 if you are using model VVORTX4.'
      ENDIF

      MODU(3)=8HVVORTX4  
      MODU(4)=W(127)
      SING0SQ=(DSIN(GAMMA0))**2 
      th0 = pid2 - Lambda0 
      SINLAM0=DSIN(LAMBDA0)
      COSTH0=SINLAM0 
      COSLAM0=DCOS(LAMBDA0)
      SINTH0=COSLAM0 
      HmSQ=Hmax**2
      WHSQ=WH**2
      DENOM = 0.D0
      IF(WHSQ.NE.0.D0) DENOM=1.D0/WHSQ
      C2 = 1.256431209
      C1 = 1.D0+1.D0/(2.D0*C2)
      C3 = C1*U0*GAMMA0
      C4 = C2/GAMMA0**2
      RETURN         
C                    
C-----------------------------------------------------------------------------
      ENTRY IPWINDg(modug)
C-----------------------------------------------------------------------------
      IF(W(125).NE.3.) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(125) should be 3 if you are using model VVORTX4.'
      ENDIF

      MODU(3)=8HVVORTX4  
      MODU(4)=W(127)
      modug%pert_name = 'VVORTX4'
      modug%pert_id = W(127)

      SING0SQ=(DSIN(GAMMA0))**2 
      th0 = pid2 - Lambda0 
      SINLAM0=DSIN(LAMBDA0)
      COSTH0=SINLAM0 
      COSLAM0=DCOS(LAMBDA0)
      SINTH0=COSLAM0 
      HmSQ=Hmax**2
      WHSQ=WH**2
      DENOM = 0.D0
      IF(WHSQ.NE.0.D0) DENOM=1.D0/WHSQ
      C2 = 1.256431209
      C1 = 1.D0+1.D0/(2.D0*C2)
      C3 = C1*U0*GAMMA0
      C4 = C2/GAMMA0**2
      RETURN         
C                    
C-----------------------------------------------------------------------------
      ENTRY PWINDR
C-----------------------------------------------------------------------------
      H=R-EARTHR
      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) 


      GAMMASQ=GAM**2
      SINGSQ=SINGAM**2 



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

      VR=0.D0
      PVRR=0.D0
      PVRTH = 0.D0
      PVRPH = 0.D0

      IF(SINGAM.NE.0.D0) THEN

         A=C3/GAM
         D=DEXP(-C4*GAMMASQ)
         B=(1.D0-D)
         C=DEXP(-DENOM*(H-Hmax)**2)
         VH=A*B*C
         PAPG=-A/GAM
         PBPG=2.D0*C4*GAM*D
         PVHPG=-(A*PBPG+B*PAPG)*C/SINGAM
         PCPR =-2.D0*(H-Hmax)*DENOM*C
         PVHR=A*B*PCPR

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

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

         PVHTH=-PVHPG*SGPGTH 
         PVHPH=-PVHPG*SGPGPH 

         PVPHR= PVHR*COSALP
         PVPHTH=PVHTH*COSALP+VH*PCOSALPTH 
         PVPHPH=PVHPH*COSALP+VH*PCOSALPPH 

         PVTHR=-PVHR*SINALP
         PVTHTH=-PVHTH*SINALP-VH*PSINALPTH 
         PVTHPH=-PVHPH*SINALP-VH*PSINALPPH 

      ELSE 
         VTH=0.d0
         VPH=0.D0
         PVTHR=0.d0
         PVPHR=0.d0
         PVTHTH=0.D0
         PVPHTH=C3*C4*C
         PVTHPH=-C3*C4*C*SINTH0
         PVPHPH=0.D0
      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

      V0 = DSQRT(VSQ0)

      IF (V0 /= 0.D0) THEN
         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
      ELSE
         PVT0 = 0.d0
         PVR0 = 0.d0
         PVTH0 = 0.d0
         PVPH0 = 0.d0
      ENDIF

C      print *, ' pvr0 = ', pvr0
C      print *, ' pvth0 = ', pvth0
C      print *, ' pvph0 = ', pvph0
C
      RETURN
C-----------------------------------------------------------------------------
      ENTRY PWINDg(wind)
C-----------------------------------------------------------------------------
      H=R-EARTHR
      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) 


      GAMMASQ=GAM**2
      SINGSQ=SINGAM**2 



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

      VR=0.D0
      PVRR=0.D0
      PVRTH = 0.D0
      PVRPH = 0.D0

      IF(SINGAM.NE.0.D0) THEN

         A=C3/GAM
         D=DEXP(-C4*GAMMASQ)
         B=(1.D0-D)
         C=DEXP(-DENOM*(H-Hmax)**2)
         VH=A*B*C
         PAPG=-A/GAM
         PBPG=2.D0*C4*GAM*D
         PVHPG=-(A*PBPG+B*PAPG)*C/SINGAM
         PCPR =-2.D0*(H-Hmax)*DENOM*C
         PVHR=A*B*PCPR

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

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

         PVHTH=-PVHPG*SGPGTH 
         PVHPH=-PVHPG*SGPGPH 

         PVPHR= PVHR*COSALP
         PVPHTH=PVHTH*COSALP+VH*PCOSALPTH 
         PVPHPH=PVHPH*COSALP+VH*PCOSALPPH 

         PVTHR=-PVHR*SINALP
         PVTHTH=-PVHTH*SINALP-VH*PSINALPTH 
         PVTHPH=-PVHPH*SINALP-VH*PSINALPPH 

      ELSE 
         VTH=0.d0
         VPH=0.D0
         PVTHR=0.d0
         PVPHR=0.d0
         PVTHTH=0.D0
         PVPHTH=C3*C4*C
         PVTHPH=-C3*C4*C*SINTH0
         PVPHPH=0.D0
      ENDIF
C   
      PVRT=0.d0
      PVTHT=0.d0
      PVPHT=0.d0

            wind1%r%v =     VR 
             wind1%r%p%tm =   PVRT 
             wind1%r%p%r  =   PVRR 
             wind1%r%p%th =  PVRTH 
             wind1%r%p%ph =  PVRPH 
         wind1%r%pp%tm%tm = 0.d0
         wind1%r%pp%tm%r  = 0.d0
         wind1%r%pp%tm%th = 0.d0
         wind1%r%pp%tm%ph = 0.d0
          wind1%r%pp%r%tm = 0.d0
          wind1%r%pp%r%r  = 0.d0
          wind1%r%pp%r%th = 0.d0
          wind1%r%pp%r%ph = 0.d0
         wind1%r%pp%th%tm = 0.d0
         wind1%r%pp%th%r  = 0.d0
         wind1%r%pp%th%th = 0.d0
         wind1%r%pp%th%ph = 0.d0
         wind1%r%pp%ph%tm = 0.d0
         wind1%r%pp%ph%r  = 0.d0
         wind1%r%pp%ph%th = 0.d0
         wind1%r%pp%ph%ph = 0.d0

           wind1%th%v =    VTH 
            wind1%th%p%tm =  PVTHT 
            wind1%th%p%r  =  PVTHR 
            wind1%th%p%th = PVTHTH 
            wind1%th%p%ph = PVTHPH 
        wind1%th%pp%tm%tm = 0.d0
        wind1%th%pp%tm%r  = 0.d0
        wind1%th%pp%tm%th = 0.d0
        wind1%th%pp%tm%ph = 0.d0
         wind1%th%pp%r%tm = 0.d0
         wind1%th%pp%r%r  = 0.d0
         wind1%th%pp%r%th = 0.d0
         wind1%th%pp%r%ph = 0.d0
        wind1%th%pp%th%tm = 0.d0
        wind1%th%pp%th%r  = 0.d0
        wind1%th%pp%th%th = 0.d0
        wind1%th%pp%th%ph = 0.d0
        wind1%th%pp%ph%tm = 0.d0
        wind1%th%pp%ph%r  = 0.d0
        wind1%th%pp%ph%th = 0.d0
        wind1%th%pp%ph%ph = 0.d0

           wind1%ph%v =    VPH 
            wind1%ph%p%tm =  PVPHT 
            wind1%ph%p%r  =  PVPHR 
            wind1%ph%p%th = PVPHTH 
            wind1%ph%p%ph = PVPHPH 
        wind1%ph%pp%tm%tm = 0.d0
        wind1%ph%pp%tm%r  = 0.d0
        wind1%ph%pp%tm%th = 0.d0
        wind1%ph%pp%tm%ph = 0.d0
         wind1%ph%pp%r%tm = 0.d0
         wind1%ph%pp%r%r  = 0.d0
         wind1%ph%pp%r%th = 0.d0
         wind1%ph%pp%r%ph = 0.d0
        wind1%ph%pp%th%tm = 0.d0
        wind1%ph%pp%th%r  = 0.d0
        wind1%ph%pp%th%th = 0.d0
        wind1%ph%pp%th%ph = 0.d0
        wind1%ph%pp%ph%tm = 0.d0
        wind1%ph%pp%ph%r  = 0.d0
        wind1%ph%pp%ph%th = 0.d0
        wind1%ph%pp%ph%ph = 0.d0

C-----------------------------------------------------------------------------
C don't need below now.

         wind1%A%v =    V 
          wind1%A%p%tm =  PVT 
          wind1%A%p%r  =  PVR 
          wind1%A%p%th = PVTH 
          wind1%A%p%ph = PVPH 

         wind1%A%v = 0.d0
          wind1%A%p%tm = 0.d0
          wind1%A%p%r  = 0.d0
          wind1%A%p%th = 0.d0
          wind1%A%p%ph = 0.d0

      wind1%A%pp%tm%tm = 0.d0
      wind1%A%pp%tm%r  = 0.d0
      wind1%A%pp%tm%th = 0.d0
      wind1%A%pp%tm%ph = 0.d0
       wind1%A%pp%r%tm = 0.d0
       wind1%A%pp%r%r  = 0.d0
       wind1%A%pp%r%th = 0.d0
       wind1%A%pp%r%ph = 0.d0
      wind1%A%pp%th%tm = 0.d0
      wind1%A%pp%th%r  = 0.d0
      wind1%A%pp%th%th = 0.d0
      wind1%A%pp%th%ph = 0.d0
      wind1%A%pp%ph%tm = 0.d0
      wind1%A%pp%ph%r  = 0.d0
      wind1%A%pp%ph%th = 0.d0
      wind1%A%pp%ph%ph = 0.d0

C-----------------------------------------------------------------------------
      wind2 = wind + wind1

      wind = wind2

      V0 = wind%A%v 

          VR0 =  wind%r%v 
        PVRT0 =  wind%r%p%tm 
        PVRR0 =  wind%r%p%r  
       PVRTH0 =  wind%r%p%th 
       PVRPH0 =  wind%r%p%ph 

         VTH0 =  wind%th%v 
       PVTHT0 =  wind%th%p%tm 
       PVTHR0 =  wind%th%p%r  
      PVTHTH0 =  wind%th%p%th 
      PVTHPH0 =  wind%th%p%ph 

         VPH0 =  wind%ph%v 
       PVPHT0 =  wind%ph%p%tm 
       PVPHR0 =  wind%ph%p%r  
      PVPHTH0 =  wind%ph%p%th 
      PVPHPH0 =  wind%ph%p%ph 

c      print *, ' vr0 = ', vr0
c      print *, ' vth0 = ', vth0
c      print *, ' vph0 = ', vph0

C-----------------------------------------------------------------------------
C don't need below now.

      IF (V0 /= 0.D0) THEN
         wind%A%p%tm = (VR0*PVRT0 + VTH0*PVTHT0 + VPH0*PVPHT0)/V0
         wind%A%p%r  = (VR0*PVRR0 + VTH0*PVTHR0 + VPH0*PVPHR0)/V0
         wind%A%p%th = (VR0*PVRTH0 + VTH0*PVTHTH0 + VPH0*PVPHTH0)/V0
         wind%A%p%ph = (VR0*PVRPH0 + VTH0*PVTHPH0 + VPH0*PVPHPH0)/V0
      ELSE
         wind%A%p%tm = 0.d0
         wind%A%p%r  = 0.d0
         wind%A%p%th = 0.d0
         wind%A%p%ph = 0.d0
      ENDIF
C-----------------------------------------------------------------------------

         V0 =  wind%A%v 
       PVT0 =  wind%A%p%tm 
       PVR0 =  wind%A%p%r  
      PVTH0 =  wind%A%p%th 
      PVPH0 =  wind%A%p%ph 

C      print *, ' pvr0 = ', pvr0
C      print *, ' pvth0 = ', pvth0
C      print *, ' pvph0 = ', pvph0

      RETURN
C-----------------------------------------------------------------------------
      END SUBROUTINE VVORTX4
C-----------------------------------------------------------------------------
