C------------------------------------------------------------------------------
C tblob2.for
C------------------------------------------------------------------------------
C  4 May 2009 - replaced REAL*8 with DOUBLE PRECISION
C 13 Mar 2008 - Corrected length in common block /PCONST/
C     12 Jul 2006 - version for gravity waves, but still need higher derivatives
C------------------------------------------------------------------------------
      SUBROUTINE TBLOB2
C------------------------------------------------------------------------------
C     TEMPERATURE PERTURBATION MODEL
C     MULTIPLICATIVE PERTURBATION WITH EXPONENTIALLY DECAYING
C     EFFECT IN ALL THREE DIRECTIONS. GIVE LATITUDE
C     INSTEAD OF CO-LATITUDE.

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

      TYPE :: values_and_derivs_1_2            ! 12 Jul 2006
         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                          ! 12 Jul 2006
         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
C     DELTA TEMPERATURE    225-249
      EQUIVALENCE (W(225),DTMODEL),(W(226),DTFORM),(W(227),DTID)
C
      EQUIVALENCE (C0,W(228)), (Z0,W(229)),(TBTH0,W(230))
      EQUIVALENCE (PH0,W(231)),(WZ,W(232)),(WTH,W(233)),(WPH,W(234))
      INTEGER DFMX,DFNTBL(10),DFITBL(10),DFFRMTB(10)
      DATA RECOGDT/2.d0/
      DATA DFMX/1/
      DATA DFNTBL/1,11,8*0/
      DATA DFITBL/1,9*0/
      DATA DFFRMTB/1,9*0/
C
      ENTRY SETPTMP
      DTMX=DFMX 
      CALL IMOVE(DTNTBL,DFNTBL,10)
      CALL IMOVE(DTITBL,DFITBL,10)
      CALL IMOVE(DTFRMTB,DFFRMTB,10)
      RETURN
C
      ENTRY IPTEMP
      IF(idint(RECOGDT) .NE. idint(DTMODEL))
     1    CALL RERROR('DTEMP  ','WRNG MODEL',RECOGDT)
      MODT(3)=6HTBLOB2
      MODT(4)=DTID
C
      FWZ=0.d0
      FWTH=0.d0
      FWPH=0.d0
      TH0= PID2-TBTH0
      IF(WZ.NE.0.d0) FWZ=2.d0/WZ/WZ
      IF(WTH.NE.0.d0) FWTH=2.d0/WTH/WTH
      IF(WPH.NE.0.d0) FWPH=2.d0/WPH/WPH
      RETURN
C
      ENTRY IPTEMPg(modtg)                ! 12 Jul 2006
C     This subroutine is called once for each new W array.

      IF(W(225) .NE. 2.0) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(225) should be 2 if you are using TBLOB2'
      ENDIF
C
      modtg%pert_name = 'NPTEMP'
      modtg%pert_id = W(227)
      RETURN

      ENTRY PTEMP
      IF(C0.EQ.0.d0) RETURN
      DZ=R-EARTHR-Z0
      DTH=TH-TH0
      DPH=PH-PH0
      DEXPO=0.d0
      EXPO=-0.5d0*(DZ*DZ*FWZ+DTH*DTH*FWTH+DPH*DPH*FWPH)
      IF(EXPO .GT. -200.d0) DEXPO=C0*dEXP(EXPO)
      DEL=1.d0+DEXPO
      PTR=PTR*DEL-T*DEXPO*FWZ*DZ
      PTTH=PTTH*DEL-T*DEXPO*FWTH*DTH
      PTPH=PTPH*DEL-T*DEXPO*FWPH*DPH
      T=T*DEL
      RETURN

      ENTRY PTEMPg(temperature)               ! 12 Jul 2006
      IF(C0.EQ.0.d0) RETURN
      DZ=R-EARTHR-Z0
      DTH=TH-TH0
      DPH=PH-PH0
      DEXPO=0.d0
      EXPO=-0.5d0*(DZ*DZ*FWZ+DTH*DTH*FWTH+DPH*DPH*FWPH)
      IF(EXPO .GT. -200.d0) DEXPO=C0*dEXP(EXPO)
      DEL=1.d0+DEXPO
      temperature%pr=temperature%pr*DEL
     1 -temperature%value*DEXPO*FWZ*DZ
      temperature%pth=temperature%pth*DEL
     1 -temperature%value*DEXPO*FWTH*DTH
      temperature%pph=temperature%pph*DEL
     1 -temperature%value*DEXPO*FWPH*DPH
      temperature%value=temperature%value*DEL
C Still need to add higher derivatives for gravity version.  ! 12 Jul 2006
      RETURN

      END
