C tdraft.for 
C-----------------------------------------------------------------------------
C  4 May 2009 - replaced REAL*8 with DOUBLE PRECISION
C 13 Mar 2008 - Corrected length in common block /PCONST/
C  9 Aug 2006: Changed name F to Ftemp.  Function F already exists in S8.
C 28 Jul 2006: Added entry points iptempg and ptempg
C 24 Apr 2006: changed 1 to 1.d0 
C------------------------------------------------------------------------------
      SUBROUTINE TDRAFT
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.
      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

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

      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 IPTEMPg(modtg) ! 28 Jul 2006
C------------------------------------------------------------------------------
C     This subroutine is called once for each new W array.

      IF(W(225).NE.3.) THEN   
C      IF(IDINT(W(225)).NE.3) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(225) should be 3 if you are using model TDRAFT.'
      ENDIF
C
      modtg%pert_name = 'TDRAFT'
      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
      RETURN
C
C------------------------------------------------------------------------------
      ENTRY IPTEMP
C------------------------------------------------------------------------------
C     This subroutine is called once for each new W array.

      IF(W(225).NE.3.) THEN   
C      IF(IDINT(W(225)).NE.3) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(225) should be 3 if you are using model TDRAFT.'
      ENDIF
C
      MODT(3)=8HTDRAFT  
      MODT(4)=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
      RETURN
C
C------------------------------------------------------------------------------
      ENTRY PTEMP
C------------------------------------------------------------------------------
      V=0.d0
      PVT=0.d0
      PVR=0.d0
      PVTH=0.d0
      PVPH=0.d0
      VR=0.d0
      PVRT=0.d0
      PVRR=0.d0
      PVRTH=0.d0
      PVRPH=0.d0
      VTH=0.d0
      PVTHT=0.d0
      PVTHR=0.d0
      PVTHTH=0.d0
      PVTHPH=0.d0
      VPH=0.d0
      PVPHT=0.d0
      PVPHR=0.d0
      PVPHTH=0.d0
      PVPHPH=0.d0
      VH=0.d0
      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)
      COSGAM=SINLAM0*COSTH+COSLAM0*SINTH*COSPH
      SINGSQ=1.d0-COSGAM**2 ! 24 Apr 2006
      IF(SINGSQ.LT.0.D0) THEN
         WRITE (3,*) 'SINGSQ = ', SINGSQ
         PRINT*, 'SINGSQ', SINGSQ
         STOP
      ENDIF
      SINGAM=DSQRT(SINGSQ)
      BETA1=SINGSQ/SING1SQ
      BETA2=SINGSQ/SING2SQ
      EXBET1=DEXP(-BETA1)
      EXBET2=DEXP(-BETA2)
      SUM1=H1SQ+HSQ
      SUM2=H2SQ+HSQ
      SUM1SQ=SUM1**2
      SUM2SQ=SUM2**2
      TEMP=EARSQ*H*COSGAM/RSQ
      VR1=TEMP*(W1*EXBET1/SUM1)
      VR2=TEMP*(W2*EXBET2/SUM2)
      VR=(VR1+VR2)*H
      PVRR1=2.D0*VR1/R*(EARTHR*H1SQ-HCUBE)/SUM1
      PVRR2=2.D0*VR2/R*(EARTHR*H2SQ-HCUBE)/SUM2
      PVRR=PVRR1+PVRR2
      SECGAM=1.D0/COSGAM
      SGPGTH=SINLAM0*SINTH-COSLAM0*COSTH*COSPH
      SGPGPH=COSLAM0*SINTH*SINPH
      PVRG1=-VR1*H*(SECGAM+2.D0*COSGAM/SING1SQ)
      PVRG2=-VR2*H*(SECGAM+2.D0*COSGAM/SING2SQ)
      PVRTH=SGPGTH*(PVRG1+PVRG2)
      PVRPH=SGPGPH*(PVRG1+PVRG2)
      TEMP1=1.D0-EXBET1
      TEMP2=1.D0-EXBET2
      IF(SINGAM.NE.0.D0) THEN
         TEMPB=-EARSQ/(R*SINGAM)
         VH1=TEMPB*W1*H1SQ*SING1SQ/SUM1SQ*TEMP1
         VH2=TEMPB*W2*H2SQ*SING2SQ/SUM2SQ*TEMP2
         VH=(VH1+VH2)*H
         SINALP=COSLAM0*SINPH/SINGAM
         SINALPSQ=SINALP**2
         COSALP=(SINTH*COSTH0-SINTH0*COSTH*COSPH)/SINGAM
         COSALPSQ=COSALP*COSALP
         VTH=VH*COSALP
         VPH=VH*SINALP
         PVHR1=VH1/R*(EARTHR*H1SQ-3.D0*EARTHR*HSQ-4.D0*HCUBE)/SUM1
         PVHR2=VH2/R*(EARTHR*H2SQ-3.D0*EARTHR*HSQ-4.D0*HCUBE)/SUM2
         PVHR=PVHR1+PVHR2
         PVTHR=PVHR*COSALP
         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
         PVHPG1=VH1*H*(-1.D0+2.D0*BETA1*EXBET1/TEMP1)*COSGAM/SINGSQ
         PVHPG2=VH2*H*(-1.D0+2.D0*BETA2*EXBET2/TEMP2)*COSGAM/SINGSQ
         PVHPG=PVHPG1+PVHPG2
         PVHTH=PVHPG*SGPGTH
         PVHPH=PVHPG*SGPGPH
         PVTHTH=PVHTH*COSALP+VH*PCOSALPTH
         PVTHPH=PVHPH*COSALP+VH*PCOSALPPH
         PVPHR=PVHR*SINALP
         PVPHTH=PVHTH*SINALP+VH*PSINALPTH
         PVPHPH=PVHPH*SINALP+VH*PSINALPPH
      ENDIF
         PSI1=EARSQ*W1*HSQ*SING1SQ/SUM1*TEMP1/2.D0-C1
         PSI2=EARSQ*W2*HSQ*SING2SQ/SUM2*TEMP2/2.D0-C2
         PSI=A1*Ftemp(A,B,PSI1,FPRIME1)+A2*Ftemp(A,B,PSI2,FPRIME2)
         T=T+PSI
         CSCGPPSI1G=R**2*H*VR1
         CSCGPPSI2G=R**2*H*VR2
         PPSI1R=H*W1*H1SQ*SING1SQ*TEMP1*EARSQ/SUM1SQ
         PPSI2R=H*W2*H2SQ*SING2SQ*TEMP2*EARSQ/SUM2SQ
         PPSI1TH=CSCGPPSI1G*SGPGTH
         PPSI2TH=CSCGPPSI2G*SGPGTH
         PPSI1PH=CSCGPPSI1G*SGPGPH
         PPSI2PH=CSCGPPSI2G*SGPGPH
         PTR=PTR+A1*FPRIME1*PPSI1R+A2*FPRIME2*PPSI2R
         PTTH=PTTH+A1*FPRIME1*PPSI1TH+A2*FPRIME2*PPSI2TH
         PTPH=PTPH+A1*FPRIME1*PPSI1PH+A2*FPRIME2*PPSI2PH
C
      RETURN
C------------------------------------------------------------------------------
      ENTRY PTEMPg(temperature) ! 28 Jul 2006
C------------------------------------------------------------------------------

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

      END SUBROUTINE TDRAFT
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
