C-----------------------------------------------------------------------------
      SUBROUTINE VTANH  
C-----------------------------------------------------------------------------
!  2 Apr 2010 - Added second derivatives
!  1 Apr 2010 - Added structures for arguments 
!  1 Apr 2010 - Added entry point IWINDg(modug)
!  1 Apr 2010 - Added entry point WINDg(wind)
C 16 Mar 2010 - Increased dimension of ..GP.. arrays to 611
C  4 May 2009 - replaced REAL*8 with DOUBLE PRECISION
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 10 Mar  2008 - corrected length errors in blank common
c 18 June 2004 - subroutine works correctly, testing commented out
c  3 June 2004, 10 pm
C-----------------------------------------------------------------------------
C     WIND/CURRENT VELOCITY MODEL 
C     VERTICAL PROFILES FOR HORIZONTAL WIND SPEED AND DIRECTION 
C     USING TANH HYPERBOLIC TANGENT PROFILES.

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

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

C     COMMON DECK "RKAM" INSERTED HERE        
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)  
      double precision kr,kth,kph                                  ! 10 Mar 2008

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)
      EQUIVALENCE (EARTHR,W(1))
C     WIND        100-124      
C
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(611)      ! 16 Mar 2010
      EQUIVALENCE (UGP,IDSU),(ANU,UGP(11))  

C.......The changes from the Cyber version indicated in lower case 
C.......below enlarge the data arrays to accommodate 100 profile points.
      DOUBLE PRECISION ca(100),uma(99),za(99),dla(99),DLCA(100)
     1 ,ALCA(100),DCA(100)
      DOUBLE PRECISION cB(100),umB(99),zB(99),dlB(99),DLCB(100)
     1 ,ALCB(100),DCB(100)
      EQUIVALENCE (Z0A,UGP(12)),(CA,ugp(112))
      EQUIVALENCE (Z0B,UGP(312)),(CB,ugp(412))
      EQUIVALENCE (ZA,UGP(13)),(UMA,ugp(113)),(DLA,ugp(213))
      EQUIVALENCE (ZB,UGP(313)),(UMB,ugp(413)),(DLB,ugp(513))
      INTEGER UPX    ,UQTBL(10),ULTBL(10),UIRMTBL(10)       
      DATA UPX/2/                                           
      DATA UQTBL/1,11,612,7*0/    
      DATA ULTBL/1,100,8*0/       
      DATA UIRMTBL/1,2,8*0/     
      DATA RECOGU,N/10.d0,0/
      DATA AQU/0.d0/
C                             
C-----------------------------------------------------------------------------

      ENTRY SETWND            
      ANU=AQU       
      UMX=UPX                 
      CALL IMOVE(UNTBL,UQTBL,10)     
      CALL IMOVE(UITBL,ULTBL,10)     
      CALL IMOVE(UFRMTBL,UIRMTBL,10) 
      CALL SETPWN                    

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

      RETURN                         
C                                    
C-----------------------------------------------------------------------------

      ENTRY IWINDR                   
      IF(W(100) .NE. 10.) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(100) should be 10 if you are using model VTANH.'
      ENDIF

      MODU(1)=8HVTANH     
      MODU(2)= w(102)
C     IF HAD PREVIOUS CALL BUT NOTHING THIS TIME, EXIT NOW
C     RETAINING PREVIOUS TABULAR DATA COUNT
      IF(N.GT.0 .AND. ANU.EQ.0.d0) RETURN
C
      N=(ANU+1)/6 - 2
      avl=float(n)
      IF(N .LE. 0)
     1    CALL RERROR('VTANH','BAD N VALUE',avl)
C
      ANU=0.d0
C
C     CONVERT 'UMA' ARRAY INPUT(OVERLAYS 'CA' ARRAY) TO 'CA' ARRAY
      U0A=CA(1)
      CALL ITANH(N,Z0A,ZA,U0A,UMA,DLA,CA,CAVA,DLCA,ALCA,DCA)
C
C     CONVERT 'UMB' ARRAY INPUT(OVERLAYS 'CB' ARRAY) TO 'CB' ARRAY
      U0B=CB(1)
      CALL ITANH(N,Z0B,ZB,U0B,UMB,DLB,CB,CAVB,DLCB,ALCB,DCB)
C
      CALL IPWINDR    
      RETURN    
C        
C-----------------------------------------------------------------------------

      ENTRY IWINDg(modug)                   
      IF(W(100) .NE. 10.) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(100) should be 10 if you are using model VTANH.'
      ENDIF

      MODU(1)=8HVTANH     
      MODU(2)= w(102)
      modug%model_name = 'VTANH'
      modug%data_id = W(102)

C     IF HAD PREVIOUS CALL BUT NOTHING THIS TIME, EXIT NOW
C     RETAINING PREVIOUS TABULAR DATA COUNT
      IF(N.GT.0 .AND. ANU.EQ.0.d0) RETURN
C
      N=(ANU+1)/6 - 2
      avl=float(n)
      IF(N .LE. 0)
     1    CALL RERROR('VTANH','BAD N VALUE',avl)
C
      ANU=0.d0
C
C     CONVERT 'UMA' ARRAY INPUT(OVERLAYS 'CA' ARRAY) TO 'CA' ARRAY
      U0A=CA(1)
      CALL ITANH2(N,Z0A,ZA,U0A,UMA,DLA,CA,CAVA,DLCA,ALCA,DCA)
C
C     CONVERT 'UMB' ARRAY INPUT(OVERLAYS 'CB' ARRAY) TO 'CB' ARRAY
      U0B=CB(1)
      CALL ITANH2(N,Z0B,ZB,U0B,UMB,DLB,CB,CAVB,DLCB,ALCB,DCB)
C
      CALL IPWINDg(modug)    
      RETURN    
C        
C-----------------------------------------------------------------------------

       ENTRY WINDR   
      CALL CLEAR(V,20)
C
      H = R - EARTHR
C     CALULATE WIND SPEED
      CALL FTANH(H,V,PVR,N,Z0A,ZA,U0A,DLA,CA,CAVA,DLCA,ALCA,DCA)
C      
C     CALCULATE WIND DIRECTION
C alpha is an angle measured clockwise from north.  It gives the direction 
c from which the wind is blowing.  
c VTH is positive if the wind is blowing toward the south.
c VPH is positive if the wind is blowing toward the east.
      CALL FTANH(H,ALPHA,PALPHAR,N,Z0B,ZB,U0B,DLB,CB,CAVB,DLCB,ALCB,DCB)
      COSALP=DCOS(ALPHA)
      SINALP=DSIN(ALPHA)
      VTH = V*COSALP
      VPH = -V*SINALP
      PVTHR = PVR*COSALP+VPH*PALPHAR
      PVPHR = -PVR*SINALP-VTH*PALPHAR
C
      CALL PWINDR  
C****************************************************************************
C This allows ALPHA and its height derivative to be plotted in the PROFILE !*
C program as Southward and Eastward components of velocity, but it cannot  !*
C be used for ray tracing.                                                 !*
C      VTH = ALPHA*1.d-10 ! 18 June 2004                                   !*
C      VPH = PALPHAR*1.d-10 ! 18 June 2004                                 !*
C****************************************************************************
       RETURN   
C-----------------------------------------------------------------------------

       ENTRY WINDg(wind)   
      CALL CLEAR(V,20)
C
      wind = wind0

      H = R - EARTHR

C     CALULATE WIND SPEED
      CALL FTANH2(H,V,PVR,pvprr,N,Z0A,ZA,U0A,DLA,CA,CAVA,DLCA,ALCA,DCA)
C      
C     CALCULATE WIND DIRECTION
C alpha is an angle measured clockwise from north.  It gives the direction 
c from which the wind is blowing.  
c VTH is positive if the wind is blowing toward the south.
c VPH is positive if the wind is blowing toward the east.
      CALL FTANH2(H,ALPHA,PALPHAR,palpharr,
     &     N,Z0B,ZB,U0B,DLB,CB,CAVB,DLCB,ALCB,DCB)
      COSALP=DCOS(ALPHA)
      SINALP=DSIN(ALPHA)
      VTH = V*COSALP
      VPH = -V*SINALP
      PVTHR = PVR*COSALP+VPH*PALPHAR
      PVPHR = -PVR*SINALP-VTH*PALPHAR
C
         wind%size%value = V
          wind%size%p%r  = PVR
           wind%th%value = VTH
           wind%ph%value = VPH
            wind%th%p%r  = PVTHR
            wind%ph%p%r  = PVPHR

       wind%size%pp%r%r  = PVPRR                                  !  2 Apr 2010
         wind%th%pp%r%r  = PVPRR*COSALP-PVR*SINALP*PALPHAR 
     & +PVPHR*PALPHAR+VPH*PALPHARR                                !  2 Apr 2010
         wind%ph%pp%r%r  = -PVPRR*SINALP-PVR*COSALP*PALPHAR 
     & -PVTHR*PALPHAR-VTH*PALPHARR                                !  2 Apr 2010

      CALL PWINDg(wind)  
C****************************************************************************
C This allows ALPHA and its height derivative to be plotted in the PROFILE !*
C program as Southward and Eastward components of velocity, but it cannot  !*
C be used for ray tracing.                                                 !*
C      VTH = ALPHA*1.d-10 ! 18 June 2004                                   !*
C      VPH = PALPHAR*1.d-10 ! 18 June 2004                                 !*
C****************************************************************************
       RETURN   
C-----------------------------------------------------------------------------
       END SUBROUTINE VTANH  
C-----------------------------------------------------------------------------
