C tplume.for 
C ----------------------------------------------------------------------!-------
! 26 Oct 2010: Added alternative definition for value of F
! 26 Oct 2010: Added common block WW to subroutine FUNC
! 26 Oct 2010: Put the calculation of F in a separate subroutine
! 25 Oct 2010: Generalized functional dependence on gamma
! 25 Oct 2010: Correctly initialized temperature perturbation and derivatives
!  9 Sep 2010: Increased to 2 plumes
!  9 Sep 2010: Moved the plume calculation to a separate subroutine
!              in anticipation of adding multiple plumes
C 10 Apr 2009: Compiled
C 26 Feb 2009: Copied stuff from VTANH.FOR
C 24 Feb 2009: Copied from tdraft2.for
C ----------------------------------------------------------------------!-------
      SUBROUTINE TPLUME
C ----------------------------------------------------------------------!-------
C     TEMPERATURE PERTURBATION MODEL
C     TEMPERATURE PERTURBATION OF A PLUME
C     LOCALIZED AT A SPECIFIED LONGITUDE AND LATITUDE
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)

      TYPE :: values_and_derivs_1_2
         DOUBLE PRECISION :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 
      TYPE (values_and_derivs_1_2) :: temperature

      TYPE :: model
         CHARACTER :: model_name
         REAL(KIND=8) :: data_id
         CHARACTER :: pert_name
         REAL(KIND=8) :: pert_id
      END TYPE model
      TYPE (model) :: modtg

      REAL(KIND=8) :: h, z1, delta1, z2, delta2, c1, c2

      TYPE :: value_and_derivs_1_dim
         REAL(KIND=8) :: value
         REAL(KIND=8) :: p,pp
      END TYPE 

      TYPE (value_and_derivs_1_dim) :: psiA1, psiB1sg, psiC1
      TYPE (value_and_derivs_1_dim) :: psiA2, psiB2sg, psiC2

      TYPE :: value_and_derivs_2_dim
         REAL(KIND=8) :: value
         REAL(KIND=8) :: h,g
         REAL(KIND=8) :: hh,hg,gg
C                        g means d/d cos gamma
      END TYPE 

      TYPE (value_and_derivs_2_dim) :: psi1sg, psi2sg, psisg

      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
      DOUBLE PRECISION :: KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)
C
C     COMMON DECK "TT" INSERTED HERE
      DOUBLE PRECISION MODT
      COMMON/TT/MODT(4),T,PTT,PTR,PTTH,PTPH

C     COMMON DECK "WW" INSERTED HERE
      COMMON/WW/ dum(10),MAXW,W(1000)
      EQUIVALENCE (EARTHR,W(1))
C     DELTA TEMPERATURE    225-249
      EQUIVALENCE (W(225),DTMODEL),(W(226),DTFORM),(W(227),DTID)
      EQUIVALENCE (LAMBDA0,W(228)),(PH0,W(229))
     &           ,(LAMBDA1,W(230)),(PH1,W(231))                  !  9 Sep 2010

      DOUBLE PRECISION LAMBDA0,LAMBDA1                            !  9 Sep 2010
C     DOUBLE PRECISION LATITUDE,LONGITUDE 
C
C     COMMON DECK "B6" INSERTED HERE                        
      INTEGER DTMX,DTNTBL,DTITBL,DTFRMTB,IDSDT(10)              
      COMMON/B6/DTMX,DTNTBL(10),DTITBL(10),DTFRMTB(10),DTGP(611) 
      EQUIVALENCE (DTGP,IDSDT),(ANDT,DTGP(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
     1   ca(100),dTma(99),za(99),dla(99),DLCA(100),ALCA(100),DCA(100)
     2 , cB(100),dTmB(99),zB(99),dlB(99),DLCB(100),ALCB(100),DCB(100)
      EQUIVALENCE (Z0A,dTGP(12)),(CA,dTgp(112))
      EQUIVALENCE (Z0B,dTGP(312)),(CB,dTgp(412))
      EQUIVALENCE (ZA,dTGP(13)),(dTMA,dTgp(113)),(DLA,dTgp(213))
      EQUIVALENCE (ZB,dTGP(313)),(dTMB,dTgp(413)),(DLB,dTgp(513))
      INTEGER dTPX    ,dTQTBL(10),dTLTBL(10),dTIRMTB(10)       
      DATA dTPX/2/                                           
      DATA dTQTBL/1,11,612,7*0/    
      DATA dTLTBL/1,100,8*0/       
      DATA dTIRMTB/1,2,8*0/     
      DATA RECOGdT,N/10.d0,0/
      DATA AQdT/0.d0/

C ----------------------------------------------------------------------!-------
      ENTRY SETPTMP
C ----------------------------------------------------------------------!-------
      ANdT=AQdT
      dTMX=dTPX 
      CALL IMOVE(dTNTBL,dTQTBL,10)
      CALL IMOVE(dTITBL,dTLTBL,10)
      CALL IMOVE(dTFRMTB,dTIRMTB,10)
      RETURN

C ----------------------------------------------------------------------!-------
      ENTRY IPTEMP
C ----------------------------------------------------------------------!-------
      IF(W(225).NE.5.) THEN   
C      IF(IDINT(W(225)).NE.5) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(225) should be 5 if you are using model TPLUME.'
      ENDIF
C
      MODT(3)='TPLUME'  
      MODT(4)=W(227)
      EARSQ=EARTHR**2

!      print *, ' lambda0 = ', lambda0
!      print *, ' ph0 = ', ph0
      TH0 = PID2 - LAMBDA0 
      SINLAM0=DSIN(LAMBDA0)
      COSTH0=SINLAM0
      COSLAM0=DCOS(LAMBDA0)
      SINTH0=COSLAM0

!      print *, ' lambda1 = ', lambda1
!      print *, ' ph1 = ', ph1
      TH1 = PID2 - LAMBDA1                                        !  9 Sep 2010
      SINLAM1=DSIN(LAMBDA1)
      COSTH1=SINLAM1
      COSLAM1=DCOS(LAMBDA1)
      SINTH1=COSLAM1

C     IF HAD PREVIOUS CALL BUT NOTHING THIS TIME, EXIT NOW
C     RETAINING PREVIOUS TABULAR DATA COUNT
      IF(N.GT.0 .AND. ANdT.EQ.0.d0) RETURN
C
      N=(ANdT+1)/6 - 2
      avl=float(n)
      IF(N .LE. 0)
     1    CALL RERROR('TPLUME','BAD N VALUE',avl)
C
      ANdT=0.d0
C
C     CONVERT 'dTMA' ARRAY INPUT(OVERLAYS 'CA' ARRAY) TO 'CA' ARRAY
      dT0A=CA(1)
      CALL ITANH(N,Z0A,ZA,dT0A,dTMA,DLA,CA,CAVA,DLCA,ALCA,DCA)
C
C     CONVERT 'dTMB' ARRAY INPUT(OVERLAYS 'CB' ARRAY) TO 'CB' ARRAY
      dT0B=CB(1)
      CALL ITANH(N,Z0B,ZB,dT0B,dTMB,DLB,CB,CAVB,DLCB,ALCB,DCB)
C
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY IPTEMPg(modtg) 
C ----------------------------------------------------------------------!-------
C     This subroutine is called once for each new W array.

C
      IF(W(225).NE.5.) THEN   
C      IF(IDINT(W(225)).NE.5) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(225) should be 5 if you are using model TPLUME.'
      ENDIF
C
      modtg%pert_name = 'TPLUME'
      modtg%pert_id = W(227)
      EARSQ=EARTHR**2

      TH0 = PID2 - LAMBDA0 
      SINLAM0=DSIN(LAMBDA0)
      COSTH0=SINLAM0
      COSLAM0=DCOS(LAMBDA0)
      SINTH0=COSLAM0

      TH1 = PID2 - LAMBDA1                                        !  9 Sep 2010
      SINLAM1=DSIN(LAMBDA1)
      COSTH1=SINLAM1
      COSLAM1=DCOS(LAMBDA1)
      SINTH1=COSLAM1

c NEED TO ADD SOME CALLS TO ITANH2 HERE

      RETURN

C ----------------------------------------------------------------------!-------
      ENTRY PTEMP
C ----------------------------------------------------------------------!-------
      H=R-EARTHR

      CALL plume(deltaT,delPTR,delPTTH,delPTPH,
     & H,TH,PH-PH0,TH0,SINTH0,COSLAM0,COSTH0,SINLAM0,
     & N,Z0A,ZA,DT0A,DLA,CA,CAVA,DLCA,ALCA,DCA,
     &   Z0B,ZB,DT0B,DLB,CB,CAVB,DLCB,ALCB,DCB)                   !  9 Sep 2010

         T=T+deltaT
         PTR =PTR + delPTR 
         PTTH=PTTH + delPTTH 
         PTPH=PTPH + delPTPH 

      if(ph1.eq.ph0.and.th1.eq.th0) return                        !  9 Sep 2010
!      print *, ' lambda0 = ', lambda0
!      print *, ' ph0 = ', ph0
!      print *, ' lambda1 = ', lambda1
!      print *, ' ph1 = ', ph1

      CALL plume(deltaT,delPTR,delPTTH,delPTPH,
     & H,TH,PH-PH1,TH1,SINTH1,COSLAM1,COSTH1,SINLAM1,
     & N,Z0A,ZA,DT0A,DLA,CA,CAVA,DLCA,ALCA,DCA,
     &   Z0B,ZB,DT0B,DLB,CB,CAVB,DLCB,ALCB,DCB)                   !  9 Sep 2010

         T=T+deltaT
         PTR =PTR + delPTR 
         PTTH=PTTH + delPTTH 
         PTPH=PTPH + delPTPH 

      RETURN

C ----------------------------------------------------------------------!-------
      ENTRY PTEMPg(temperature) 
C ----------------------------------------------------------------------!-------

C Still need to add above formulas plus higher derivatives for gravity version
      RETURN

      END SUBROUTINE TPLUME
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      SUBROUTINE plume(deltaT,delPTR,delPTTH,delPTPH,
     & H,TH,PH,TH0,SINTH0,COSLAM0,COSTH0,SINLAM0,
     & N,Z0A,ZA,DT0A,DLA,CA,CAVA,DLCA,ALCA,DCA,
     &   Z0B,ZB,DT0B,DLB,CB,CAVB,DLCB,ALCB,DCB)                   !  9 Sep 2010
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)

      DOUBLE PRECISION
     1   ca(100),dTma(99),za(99),dla(99),DLCA(100),ALCA(100),DCA(100)
     2 , cB(100),dTmB(99),zB(99),dlB(99),DLCB(100),ALCB(100),DCB(100)

      COSTH=DCOS(TH)
      SINTH=DSIN(TH)
      COSPH=DCOS(PH)
      SINPH=DSIN(PH)
C      COSGAM=SINLAM0*COSTH+COSLAM0*SINTH*COSPH
      call bettercea(th0,th,gam,costh0,costh,cosgam,
     1 sinth0,sinth,singam,ph) 


      SINGSQ=SINGAM**2 


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

      deltaT=0.d0                                                 ! 25 Oct 2010
      delPTR=0.d0                                                 ! 25 Oct 2010
      delPTTH=0.d0                                                ! 25 Oct 2010
      delPTPH=0.d0                                                ! 25 Oct 2010

C     CALULATE STRENGTH OF TEMPERATURE PERTURBATION
      CALL FTANH(H,delT,PdelTR,N,Z0A,ZA,DT0A,DLA,CA,CAVA,DLCA,ALCA,DCA)
         if(delT.eq.0.d0) return
C      
C     CALCULATE WIDTH OF TEMPERATURE PERTURBATION
      CALL FTANH(H,WID,PWIDR,N,Z0B,ZB,DT0B,DLB,CB,CAVB,DLCB,ALCB,DCB)
         if(wid.eq.0.d0) return
         x=gam/wid                                                ! 25 Oct 2010
         CALL FUNC(x,F,dFdx)                                      ! 26 Oct 2010
         deltaT=delT*F
         pxpr = -gam*pwidr/wid**2                                 ! 25 Oct 2010
         pFpr = dFdx*pxpr                                         ! 25 Oct 2010
         delPTR = F*pdelTr +delT*pFpr                             ! 25 Oct 2010
         IF(SINGAM.EQ.0.) RETURN
         pxpCOSgam = -1.d0/(wid*singam)                           ! 25 Oct 2010
         pFpCOSgam = dFdx*pxpCOSgam                               ! 25 Oct 2010
         delPTTH = -SGPGTH*delT*pFpCOSgam                         ! 25 Oct 2010
         delPTPH = -SGPGPH*delT*pFpCOSgam                         ! 25 Oct 2010
C
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PLUME
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      SUBROUTINE FUNC(x,F,dFdx)                                   ! 26 Oct 2010
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)

C     COMMON DECK "WW" INSERTED HERE
      COMMON/WW/ dum(10),MAXW,W(1000)
      EQUIVALENCE (EARTHR,W(1))
C     DELTA TEMPERATURE    225-249
      EQUIVALENCE (W(232),a2),(W(233),a4),(W(234),a6),(W(235),a8) ! 26 Oct 2010
C
      IF(a2.eq.0.d0.and.a4.eq.0.d0.and.a6.eq.0.d0.and.a8.eq.0.d0) then
         F = dexp(-x**2)                                          ! 25 Oct 2010
         dFdx = -2.d0*x*F                                         ! 25 Oct 2010
      ELSE
         F = dexp(-x**2*(a2+x**2*(a4+x**2*(a6+a8*x**2))))         ! 26 Oct 2010
       dFdx=F*(-2.d0*x*(a2+x**2*(2.d0*a4+x**2*(3.d0*a6+4.d0*a8*x**2))))
      ENDIF                                                       ! 26 Oct 2010

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE FUNC
C ----------------------------------------------------------------------!-------


