
C-----------------------------------------------------------------------------
      SUBROUTINE WLINEAR    
C-----------------------------------------------------------------------------
!  6 Jan 2011 - Commented out calls to IPWINDg and PWINDg
!  1 Apr 2010 - Added structures for arguments and tested them
! 31 Mar 2010 - Added entry point IWINDg(modug)
! 31 Mar 2010 - Added entry point WINDg(wind)
C-----------------------------------------------------------------------------
C     LINEAR WIND VELOCITY PROFILE   
C     PROVIDES CONSTANT RADIAL, ZONAL AND MERIDONAL WINDS  
C     EXCEPT THAT A POSSIBLE LINEAR HEIGHT GRADIENT OF THE ZONAL
C     COMPONENT IS ALLOWED.                   

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

C     COMMON DECK "RKAM" INSERTED HERE        
      DOUBLE PRECISION KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)

C     COMMON DECK "UU" INSERTED HERE  
      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)

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------------------------------------------------------------------------------
! A structure for a value plus first and second derivatives:      ! 31 Mar 2010

      TYPE :: gradient
         DOUBLE PRECISION :: tm,r,th,ph
      END TYPE gradient

      TYPE :: gradient2
         TYPE (gradient) :: tm,r,th,ph
      END TYPE gradient2

      TYPE :: value_and_gradients_1_2
         DOUBLE PRECISION :: value
         TYPE (gradient) :: p
         TYPE (gradient2) :: pp
      END TYPE value_and_gradients_1_2

! A structure for a vector plus first and second derivatives:     ! 31 Mar 2010

      TYPE :: vector_and_gradients_1_2
         TYPE  (value_and_gradients_1_2) :: size,r,th,ph
      END TYPE vector_and_gradients_1_2

      TYPE (vector_and_gradients_1_2) wind, wind0

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

C     COMMON DECK "B1" INSERTED HERE                
      INTEGER UMX,UNTBL,UITBL,UFRMTBL,IDSU(10)      
      COMMON/B1/UMX,UNTBL(10),UITBL(10),UFRMTBL(10),UGP(10)    
      EQUIVALENCE (UGP,IDSU),(earthr,w(1))
      EQUIVALENCE (W(103),UR0)            
      EQUIVALENCE (W(104),UTH0),(W(105),UPH0),(W(106),PUPHZ0)
      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                              
C-----------------------------------------------------------------------------
C                              
      ENTRY SETWND             
      UMX=UPX       
      CALL IMOVE(UNTBL,UQTBL,10)
      CALL IMOVE(UITBL,ULTBL,10)
      CALL IMOVE(UFRMTBL,UIRMTBL,10)
      CALL SETPWN                   
C                                   
! The following works, but it is long.  An equivalent, shorter method follows.

!         wind0%size%value = 0.d0
!          wind0%size%p%tm = 0.d0
!          wind0%size%p%r  = 0.d0
!          wind0%size%p%th = 0.d0
!          wind0%size%p%ph = 0.d0
!      wind0%size%pp%tm%tm = 0.d0
!      wind0%size%pp%tm%r  = 0.d0
!      wind0%size%pp%tm%th = 0.d0
!      wind0%size%pp%tm%ph = 0.d0
!       wind0%size%pp%r%tm = 0.d0
!       wind0%size%pp%r%r  = 0.d0
!       wind0%size%pp%r%th = 0.d0
!       wind0%size%pp%r%ph = 0.d0
!      wind0%size%pp%th%tm = 0.d0
!      wind0%size%pp%th%r  = 0.d0
!      wind0%size%pp%th%th = 0.d0
!      wind0%size%pp%th%ph = 0.d0
!      wind0%size%pp%ph%tm = 0.d0
!      wind0%size%pp%ph%r  = 0.d0
!      wind0%size%pp%ph%th = 0.d0
!      wind0%size%pp%ph%ph = 0.d0
!
!            wind0%r%value = 0.d0
!             wind0%r%p%tm = 0.d0
!             wind0%r%p%r  = 0.d0
!             wind0%r%p%th = 0.d0
!             wind0%r%p%ph = 0.d0
!         wind0%r%pp%tm%tm = 0.d0
!         wind0%r%pp%tm%r  = 0.d0
!         wind0%r%pp%tm%th = 0.d0
!         wind0%r%pp%tm%ph = 0.d0
!          wind0%r%pp%r%tm = 0.d0
!          wind0%r%pp%r%r  = 0.d0
!          wind0%r%pp%r%th = 0.d0
!          wind0%r%pp%r%ph = 0.d0
!         wind0%r%pp%th%tm = 0.d0
!         wind0%r%pp%th%r  = 0.d0
!         wind0%r%pp%th%th = 0.d0
!         wind0%r%pp%th%ph = 0.d0
!         wind0%r%pp%ph%tm = 0.d0
!         wind0%r%pp%ph%r  = 0.d0
!         wind0%r%pp%ph%th = 0.d0
!         wind0%r%pp%ph%ph = 0.d0
!
!           wind0%th%value = 0.d0
!            wind0%th%p%tm = 0.d0
!            wind0%th%p%r  = 0.d0
!            wind0%th%p%th = 0.d0
!            wind0%th%p%ph = 0.d0
!        wind0%th%pp%tm%tm = 0.d0
!        wind0%th%pp%tm%r  = 0.d0
!        wind0%th%pp%tm%th = 0.d0
!        wind0%th%pp%tm%ph = 0.d0
!         wind0%th%pp%r%tm = 0.d0
!         wind0%th%pp%r%r  = 0.d0
!         wind0%th%pp%r%th = 0.d0
!         wind0%th%pp%r%ph = 0.d0
!        wind0%th%pp%th%tm = 0.d0
!        wind0%th%pp%th%r  = 0.d0
!        wind0%th%pp%th%th = 0.d0
!        wind0%th%pp%th%ph = 0.d0
!        wind0%th%pp%ph%tm = 0.d0
!        wind0%th%pp%ph%r  = 0.d0
!        wind0%th%pp%ph%th = 0.d0
!        wind0%th%pp%ph%ph = 0.d0
!
!           wind0%ph%value = 0.d0
!            wind0%ph%p%tm = 0.d0
!            wind0%ph%p%r  = 0.d0
!            wind0%ph%p%th = 0.d0
!            wind0%ph%p%ph = 0.d0
!        wind0%ph%pp%tm%tm = 0.d0
!        wind0%ph%pp%tm%r  = 0.d0
!        wind0%ph%pp%tm%th = 0.d0
!        wind0%ph%pp%tm%ph = 0.d0
!         wind0%ph%pp%r%tm = 0.d0
!         wind0%ph%pp%r%r  = 0.d0
!         wind0%ph%pp%r%th = 0.d0
!         wind0%ph%pp%r%ph = 0.d0
!        wind0%ph%pp%th%tm = 0.d0
!        wind0%ph%pp%th%r  = 0.d0
!        wind0%ph%pp%th%th = 0.d0
!        wind0%ph%pp%th%ph = 0.d0
!        wind0%ph%pp%ph%tm = 0.d0
!        wind0%ph%pp%ph%r  = 0.d0
!        wind0%ph%pp%ph%th = 0.d0
!        wind0%ph%pp%ph%ph = 0.d0

! The following is a shorter equivalent way to zero out the structure.

      wind0%size=value_and_gradients_1_2
     &           (0.d0,gradient(0.d0,0.d0,0.d0,0.d0)
     & , gradient2(gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)))

      wind0%r=value_and_gradients_1_2
     &           (0.d0,gradient(0.d0,0.d0,0.d0,0.d0)
     & , gradient2(gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)))

      wind0%th=value_and_gradients_1_2
     &           (0.d0,gradient(0.d0,0.d0,0.d0,0.d0)
     & , gradient2(gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)))

      wind0%ph=value_and_gradients_1_2
     &           (0.d0,gradient(0.d0,0.d0,0.d0,0.d0)
     & , gradient2(gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)
     &            ,gradient(0.d0,0.d0,0.d0,0.d0)))

! The following was used for testing.

!      wind0%size=value_and_gradients_1_2
!     &           (1.d0,gradient(2.d0,3.d0,4.d0,5.d0)
!     & , gradient2(gradient(11.d0,12.d0,13.d0,14.d0)
!     &            ,gradient(21.d0,22.d0,23.d0,24.d0)
!     &            ,gradient(31.d0,32.d0,33.d0,34.d0)
!     &            ,gradient(41.d0,42.d0,43.d0,44.d0)))
!
!      wind0%r=value_and_gradients_1_2
!     &           (1,gradient(2,3,4,5)
!     & , gradient2(gradient(11,12,13,14)
!     &            ,gradient(21,22,23,24)
!     &            ,gradient(31,32,33,34)
!     &            ,gradient(41,42,43,44)))

!      wind = wind0

!      PRINT*, 'wind = ', wind                                 ! testing types

      RETURN                        

C-----------------------------------------------------------------------------
      ENTRY IWINDR                  
C-----------------------------------------------------------------------------
C     This subroutine is called once for each new W array.

      IF(W(100) .NE. 1.) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(100) should be 1 if you are using model WLINEAR.'
      endif

      MODU(1)=7HWLINEAR                               
      MODU(2)=w(102)
      CALL IPWINDR                                    
      RETURN                                          
C                                                     
C-----------------------------------------------------------------------------
      ENTRY IWINDg(modug)                  
C-----------------------------------------------------------------------------
C     This subroutine is called once for each new W array.

      IF(W(100) .NE. 1.) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(100) should be 1 if you are using model WLINEAR.'
      endif

      MODU(1)=7HWLINEAR                               
      MODU(2)=w(102)
      modug%model_name = 'WLINEAR'
      modug%data_id = W(102)
c      CALL IPWINDg(modug)                                       !  6 Jan 2011
      RETURN                                          
C                                                     
C-----------------------------------------------------------------------------
       ENTRY WINDR                                    
C-----------------------------------------------------------------------------
       H = R - EARTHR                
      CALL CLEAR(V,20)               
       VR = UR0                      
       VTH = UTH0                    
       VPH = (UPH0 + PUPHZ0 * H)     
      V=DSQRT(VR*VR+VTH*VTH+VPH*VPH) 
      PVR=PUPHZ0                     
      IF(V.NE.0.D0) PVR=VPH/V*PUPHZ0 
      PVPHR = PUPHZ0                 
      CALL PWINDR                    
       RETURN                        

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

       ENTRY WINDg(wind)                                    

      wind = wind0

       H = R - EARTHR                
      CALL CLEAR(V,20)               
       VR = UR0                      
       VTH = UTH0                    
       VPH = (UPH0 + PUPHZ0 * H)     
      V=DSQRT(VR*VR+VTH*VTH+VPH*VPH) 
      PVR=PUPHZ0                     
      IF(V.NE.0.D0) PVR=VPH/V*PUPHZ0 
      PVPHR = PUPHZ0                 

         wind%size%value = V
          wind%size%p%r  = PVR
            wind%r%value = VR
           wind%th%value = VTH
           wind%ph%value = VPH
            wind%ph%p%r  = PVPHR

! Notice that all second derivatives are zero for this model.

c      CALL PWINDg(wind)                                         !  6 Jan 2011
       RETURN                        
C-----------------------------------------------------------------------------
       END SUBROUTINE WLINEAR                              
C-----------------------------------------------------------------------------
