C tdraft2.for 
C-----------------------------------------------------------------------------
C  4 May 2009 - replaced REAL*8 with DOUBLE PRECISION
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 10 Nov 2006: Changing the formula to insure no perturbation above capping ht.
C  9 Aug 2006: Correction
C  8 Aug 2006: Added calculation of TH0
C  8 Aug 2006: Changed name F to Ftemp.  Function F already exists in S8.
C  8 Aug 2006: Correctly did the stream function in spherical polar coordinates
C  4 Aug 2006: Using better formula for central earth angle
C  4 Aug 2006: Formulas for small gamma (GAM) and small beta
C 20 Jul 2006: Added entry points iptempg and ptempg
C 10 Jul 2006: Corrected $ to % (4 places)
C 27 Apr 2006: vsi -> psi
C 24 Apr 2006: Modular version broken into subroutines
C 24 Apr 2006: Version for the updraft being capped at a user-specified height
C 24 Apr 2006: Copied from tdraft.for
C------------------------------------------------------------------------------
      SUBROUTINE TDRAFT2
C------------------------------------------------------------------------------
C     TEMPERATURE PERTURBATION MODEL
C     REPRESENTS THE TEMPERATURE PERTURBATION AS A FUNCTION OF THE
C     STREAM FUNCTION THAT GIVES AN UPDRAFT PLUS A DOWNDRAFT.
C-----------------------------------------------------------------------------
      implicit double precision (a-h,o-z)

      EXTERNAL Ftemp

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: 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                        ! 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
      DOUBLE PRECISION KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)
C     COMMON DECK "B6" INSERTED HERE
      INTEGER DTMX,DTNTBL,DTITBL,DTFRMTB,IDSDT(10)
      COMMON/B6/DTMX,DTNTBL(10),DTITBL(10),DTFRMTB(10),DTGP(10)
      EQUIVALENCE (DTGP,IDSDT)

C     COMMON DECK "TT" INSERTED HERE
      DOUBLE PRECISION MODT
      COMMON/TT/MODT(4),T,PTT,PTR,PTTH,PTPH
      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 (A1,W(228)), (A2,W(229)),(A,W(230)),(B,W(231))
      EQUIVALENCE (LAMBDA0,W(128)),(PHI0,W(129))
      EQUIVALENCE (W1,W(130)),(GAMMA1,W(131)),(H1,W(132))    
      EQUIVALENCE (W2,W(133)),(GAMMA2,W(134)),(H2,W(135))    
      EQUIVALENCE (z1,W(136)),(delta1,W(137))
      EQUIVALENCE (z2,W(138)),(delta2,W(139))

      DOUBLE PRECISION LAMBDA0
C     DOUBLE PRECISION LATITUDE,LONGITUDE 
C     COMMON DECK "B2" INSERTED HERE                        
      INTEGER DUMX,DUNTBL,DUITBL,DUFRMTBL,IDSDU(20)              
      INTEGER DFMX,DFNTBL(10),DFITBL(10),DFFRMTB(10)
      DATA DFMX/1/
      DATA DFNTBL/1,11,8*0/
      DATA DFITBL/1,9*0/
      DATA DFFRMTB/1,9*0/
C
C------------------------------------------------------------------------------
      ENTRY SETPTMP
C------------------------------------------------------------------------------
      DTMX=DFMX 
      CALL IMOVE(DTNTBL,DFNTBL,10)
      CALL IMOVE(DTITBL,DFITBL,10)
      CALL IMOVE(DTFRMTB,DFFRMTB,10)
      RETURN
C
C------------------------------------------------------------------------------
      ENTRY IPTEMP
C------------------------------------------------------------------------------
      IF(W(225).NE.4.) THEN   
C      IF(IDINT(W(225)).NE.4) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(225) should be 4 if you are using model TDRAFT2.'
      ENDIF
C
      MODT(3)=8HTDRAFT2  
      MODT(4)=W(227)
      EARSQ=EARTHR**2
      SING1SQ=(DSIN(GAMMA1))**2
      SING2SQ=(DSIN(GAMMA2))**2
      TH0 = PID2 - LAMBDA0 ! 8 Aug 2006
      SINLAM0=DSIN(LAMBDA0)
      COSTH0=SINLAM0
      COSLAM0=DCOS(LAMBDA0)
      SINTH0=COSLAM0
      H1SQ=H1**2
      H2SQ=H2**2
      C1=EARTHR*W1*SING1SQ/2.D0 ! 9 Aug 2006
      C2=EARTHR*W2*SING2SQ/2.D0 ! 9 Aug 2006
      CONST1 = earthr*w1/2.d0 ! 8 Aug 2006
      CONST2 = earthr*w2/2.d0 ! 8 Aug 2006
      RETURN
C
C------------------------------------------------------------------------------
      ENTRY IPTEMPg(modtg) ! 20 Jul 2006
C------------------------------------------------------------------------------
C     This subroutine is called once for each new W array.

C
      IF(W(225).NE.4.) THEN   
C      IF(IDINT(W(225)).NE.4) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(225) should be 4 if you are using model TDRAFT2.'
      ENDIF
C
      modtg%pert_name = 'TDRAFT2'
      modtg%pert_id = W(227)
      EARSQ=EARTHR**2
      SING1SQ=(DSIN(GAMMA1))**2
      SING2SQ=(DSIN(GAMMA2))**2
      SINLAM0=DSIN(LAMBDA0)
      COSTH0=SINLAM0
      COSLAM0=DCOS(LAMBDA0)
      SINTH0=COSLAM0
      H1SQ=H1**2
      H2SQ=H2**2
      C1=EARSQ*W1*SING1SQ/2.D0
      C2=EARSQ*W2*SING2SQ/2.D0
      CONST1 = earsq*w1/2.d0
      CONST2 = earsq*w2/2.d0
      RETURN

C------------------------------------------------------------------------------
      ENTRY PTEMP
C------------------------------------------------------------------------------
      RSQ=R**2
      H=R-EARTHR
      HSQ=H**2
      HCUBE=H*HSQ
      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) !  4 Aug 2006

      CALL psiA(h1, h, psiA1)
      CALL psiB(CONST1, sing1sq, cosgam, singam, psiB1sg) !  4 Aug 2006
      CALL psiC(z1, delta1, h, psiC1)

      CALL psiA(h2, h, psiA2)
      CALL psiB(CONST2, sing2sq, cosgam, singam, psiB2sg) !  4 Aug 2006
      CALL psiC(z2, delta2, h, psiC2)

      psi1sg%value = (psiA1%value*psiB1sg%value - c1)*psiC1%value ! 10 Nov 2006
      psi1sg%h = psiB1sg%value*(psiA1%p*psiC1%value+psiA1%value*psiC1%p)
     1 -c1*psiC1%p                                                ! 10 Nov 2006
      psi1sg%g = psiA1%value*psiB1sg%p*psiC1%value
c      psi1sg%hh = psiB1sg%value*(psiA1%pp*psiC1%value + 
c     1 2.d0*psiA1%p*psiC1%p + psiA1%value*psiC1%pp)
c     2 -c1*psiC1%pp                                              ! 10 Nov 2006
c      psi1sg%gg = psiA1%value*psiB1sg%pp*psiC1%value
c      psi1sg%hg = psiB1sg%p*(psiA1%p*psiC1%value+psiA1%value*psiC1%p)

      psi2sg%value = (psiA2%value*psiB2sg%value - c2)*psiC2%value ! 10 Nov 2006
      psi2sg%h = psiB2sg%value*(psiA2%p*psiC2%value+psiA2%value*psiC2%p)
     1 -c2*psiC2%p                                                ! 10 Nov 2006
      psi2sg%g = psiA2%value*psiB2sg%p*psiC2%value
c      psi2sg%hh = psiB2sg%value*(psiA2%pp*psiC2%value + 
c     1 2.d0*psiA2%p*psiC2%p + psiA2%value*psiC2%pp)
c     2 -c2*psiC2%pp                                              ! 10 Nov 2006
c      psi2sg%gg = psiA2%value*psiB2sg%pp*psiC2%value
c      psi2sg%hg = psiB2sg%p*(psiA2%p*psiC2%value+psiA2%value*psiC2%p)

      SINGSQ=SINGAM**2 ! 4 Aug 2006

      psisg%value = psi1sg%value + psi2sg%value
      psisg%h = psi1sg%h + psi2sg%h
      psisg%g = psi1sg%g + psi2sg%g
c      psisg%hh = psi1sg%hh + psi2sg%hh
c      psisg%hg = psi1sg%hg + psi2sg%hg
c      psisg%gg = psi1sg%gg + psi2sg%gg

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

         deltaT=A1*Ftemp(A,B*EARTHR,PSI1SG%value,FPRIME1) + 
     1          A2*Ftemp(A,B*EARTHR,PSI2SG%value,FPRIME2)
         T=T+deltaT
         PPSI1TH= -psi1sg%g *SGPGTH
         PPSI2TH= -psi2sg%g *SGPGTH
         PPSI1PH= -psi1sg%g *SGPGPH
         PPSI2PH= -psi2sg%g *SGPGPH
         PTR=PTR+A1*FPRIME1*psi1sg%h+A2*FPRIME2*psi2sg%h
         PTTH=PTTH+A1*FPRIME1*PPSI1TH+A2*FPRIME2*PPSI2TH
         PTPH=PTPH+A1*FPRIME1*PPSI1PH+A2*FPRIME2*PPSI2PH
C
      RETURN

C------------------------------------------------------------------------------
      ENTRY PTEMPg(temperature) ! 20 Jul 2006
C------------------------------------------------------------------------------

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

      END SUBROUTINE TDRAFT2
C------------------------------------------------------------------------------
      DOUBLE PRECISION FUNCTION Ftemp(A,B,X,FPRIME)
      DOUBLE PRECISION A,B,X,FPRIME
      Ftemp=A+B*X
      FPRIME=B
      RETURN
      END FUNCTION Ftemp
