C Twave.for 
C ----------------------------------------------------------------------!-------
! 27 Dec 2010 - Subroutine is working
! 15 Dec 2010 - First try
!  3 Dec 2010 - Copied from Ttanh6.for and Tplume.for
C ----------------------------------------------------------------------!-------
      SUBROUTINE Twave
C ----------------------------------------------------------------------!-------
C     TEMPERATURE PERTURBATION MODEL
C     Twave is a temperature perturbation model that has height profiles 
C     composed of three segments with hyperbolic tangent transitions 
C     that are the same as those in ttanh5 and ttanh6.  The parameters for
C     the profiles vary harmonically with latitude to emulate a gravity wave.
C     The lower segment is meant to approximately cancel the gradient of the 
C     background profile so that the total temperature in the lower segment 
C     is constant in both height and latitude.  The upper segment is meant to 
C     give approximately zero temperature perturation to the background 
C     temperature.  The middle segment gives a sharp transition from the 
C     lower segment to the upper segment.  The height of the transition 
C     varies sinusoidally with latitude.
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)

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

      TYPE :: model
         CHARACTER :: model_name
         DOUBLE PRECISION :: data_id
         CHARACTER :: pert_name
         DOUBLE PRECISION :: pert_id
      END TYPE model
      TYPE (model) :: modtg

      DOUBLE PRECISION :: h, z1, delta1, z2, delta2, c1, c2

      TYPE :: value_and_derivs_1_dim
         DOUBLE PRECISION :: value
         DOUBLE PRECISION :: 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
         DOUBLE PRECISION :: value
         DOUBLE PRECISION :: h,g
         DOUBLE PRECISION :: 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 (h0,W(228)),(deltaH,W(229))
     &           ,(deltaT,W(230)),(delta1,W(231))                
     &           ,(delta2,W(232)),(slope,W(233))                
     &           ,(A,W(234)),(kgrav0,W(235)),(phase,W(236))                

      DOUBLE PRECISION :: kgrav0,kgravTH

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

      DOUBLE PRECISION ::
     1   c(100),dTm(99),z(99),dl(99),DLC(100),ALC(100),DC(100),TC(100)
     2  ,pdelTz(0:101),pcz(0:101)

      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.6.) THEN   
C      IF(IDINT(W(225)).NE.6) THEN   
           PRINT*,' Model check number mismatch: '
           STOP ' W(225) should be 6 if you are using model Twave.'
      ENDIF
C
      MODT(3)='Twave'  
      MODT(4)=W(227)

      z0=0.d0
      z(3)=1000.0d0  
      dTm(1)=-deltaT
      dTm(2)=0.d0
      dTm(3)=0.d0
      dT0=dTm(1)+(z0-h0)*slope
      dl(1)=delta1
      dl(2)=delta2
      dl(3)=0.d0
      kgravTH=earthr*kgrav0
      n=2

      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY IPTEMPg(modtg) 
C ----------------------------------------------------------------------!-------
C     This subroutine is called once for each new W array.

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

      RETURN

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

      z(1)=h0-deltaH/2.d0+A*dcos(kgravTH*(pi/2.d0-TH)+phase)
      z(2)=h0+deltaH/2.d0+A*dcos(kgravTH*(pi/2.d0-TH)+phase)

      CALL ITANH3(N,Z0,Z,dT0,dTM,DL,C,CAV,DLC,ALC,DC
     & ,TC,PCZ)
 

      CALL ftanh3(h,delT,delpTr,pdelTz,n,z0,z,dT0,dl,c,CAV,DLC,ALC,DC
     & ,TC,PCZ)

      pzipTH = A*kgravTH*dsin(kgravTH*(pi/2.d0-TH)+phase)
 
c      PRINT *,'delT = ',delT

      delpTt = 0.d0 
      delpTth = pzipTH*(pdelTz(1)+pdelTz(2)) 
      delpTph = 0.d0 


         T=T+delT
         PTR =PTR + delPTR 
         PTTH=PTTH + delpTTH 
!         PTPH=PTPH + delPTPH 

      RETURN

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

      h = r - earthr

c NEED TO ADD SOME CALLS TO ITANH3 HERE
! (Actually, for gravity waves, we do need the second derivative of T with 
! respect to height and latitude, so we would need to generalize ftanh2 
! and ftanh3.  Neither would be enough.)

      CALL ftanh3
     & (h,delT,delpTr,pdelTz,n,z0,z,delT0,dl,c,CAV,DLC,ALC,DC
     & ,PCZ)

c      PRINT *,' at tempg, delT0 = ', delT0
      temperature%value = temperature%value + delT

      temperature%pr = temperature%pr + delpTr
      temperature%prr = temperature%prr + delpTrr

c      PRINT *,'delT = ',delT
c      PRINT *,'delpTr = ',delpTr
c      PRINT *,'delpTrr = ',temperature%prr

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

C ----------------------------------------------------------------------!-------
      END SUBROUTINE Twave
C ----------------------------------------------------------------------!-------






