C ---------------------------------------------------------------------
C     PscaleH.for
C ---------------------------------------------------------------------
! 27 Jun 2016 - Changed calling sequence from (gravity) to (gravity,coriolis)
! 30 Apr 2010 - Slight change in setting model labeling
! 30 Mar 2010 - Added CALL ppres
C 29 Mar 2010 - Added ENTRY PRES
C 29 Mar 2010 - Changed subroutine name to PscaleH
C  6 Aug 2009 - Altered some comments about the units for pressure.
C 13 Mar 2008 - Corrected length in common block /PCONST/
C      5 Dec 2006 - Added integration of pressure for acoustic dispersion
C     29 Mar 2006 - Adjusted pressure w/ or w/o integrating
C     17 Feb 2006 - Added calculation of derivatives in /pp/
C     16 Feb 2006 - Added integration to get the correct pressure
C      2 Feb 2006 - Added modp(1) and modp(2) to ENTRY ipresg()
C      2 Feb 2006 - Added CALL setppres
C     26 Jan 2006 - 
C     25 Jan 2006 - Changed name to PscaleH.for
C     24 Jan 2006 - different entry points for acoustics and gravity
C     23 Jan 2006, version for gravity waves
C ---------------------------------------------------------------------

C     SCALE HEIGHT PRESSURE MODEL
C     The pressure is in Newtons per square meter.  For purposes of a loss model 
C     subroutine that depends on pressure, this subroutine may give an 
C     approximate value based on an exponential pressure profile.    6 Aug 2009

C     However, the ratio of pressure to its gradient is always correct 
C     and will always give the correct Brunt-Vaisala frequency.      6 Aug 2009

C ---------------------------------------------------------------------
      SUBROUTINE PscaleH                                          ! 29 Mar 2010
C ---------------------------------------------------------------------

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

      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) :: pressure
      TYPE (values_and_derivs_1_2) :: temperature
      TYPE (values_and_derivs_1_2) :: molecular_weight

      TYPE :: gravity_and_derivs_1
         REAL(KIND=8) :: gs,gsr,gsth,gsph
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: prtm,prr,prth,prph
         REAL(KIND=8) :: pthtm,pthr,pthth,pthph
         REAL(KIND=8) :: pphtm,pphr,pphth,pphph
      END TYPE 
      TYPE (gravity_and_derivs_1) :: gravity
      TYPE (gravity_and_derivs_1) :: coriolis                      ! 27 Jun 2016
      REAL(KIND=8) T,g,gr,gth,gph,mw

      PARAMETER (NWARSZ=1000)
C     COMMON DECK "WW" INSERTED HERE
      COMMON/WW/ dum(10),MAXW,W(1000)
      EQUIVALENCE (EARTHR,W(1))                                   ! 25 Jan 2006

C     PRESSURE    550-574
      EQUIVALENCE (W(550),pmodel),(W(551),pform),(W(552),pid)

      REAL(KIND=8) KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)

C     COMMON DECK "CONST" INSERTED HERE
      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 "PP" INSERTED HERE
      REAL(KIND=8) p,modp
      COMMON/PP/MODP(4),p,ppt,ppr,ppth,ppph

      EQUIVALENCE (W(553),p0),(W(554),Hscale0)                    ! 25 Jan 2006
      EQUIVALENCE (W(555),z0),(W(556),eqnum)                      ! 16 Feb 2006
C ---------------------------------------------------------------------
      ENTRY pres                                                  ! 29 Mar 2010
      z = R - earthr                                              ! 25 Jan 2006
C ---------------------------------- copied from here ------- 5 Dec 2006 
c      print*,'R =',r
c      print*,'earthr =',earthr
c      print*,'z =',z
      
      CALL tempg(temperature)                                     ! 24 Jan 2006
      T = temperature%value

      CALL grav(gravity,coriolis)                                  ! 27 Jun 2016
      g = gravity%gs
      gr = gravity%gsr
      gth = gravity%gsth
      gph = gravity%gsph

      CALL molwtg(molecular_weight)                               ! 24 Jan 2006
      mw = molecular_weight%value

c         print*,' rgas = ', rgas
      hscale = (rgas*T)/(mw*g)
c      PRINT *,'hscale pressure = ',hscale
C -------------------------------------------- to here ------ 5 Dec 2006 

      if (eqnum.gt.0.d0) then
         numeq = eqnum + 0.5                               ! 29 Mar 2006
         p = p0*exp(rkvars(numeq)-(z-z0)/Hscale0)          ! 29 Mar 2006

c         print*,'p0 = ',p0
c         print*,'numeq = ',numeq
c         print*,'rkvars(numeq) = ',rkvars(numeq)
c         print*,'z = ',z
c         print*,'z0 = ',z0
c         print*,'hscale0 = ', hscale0
c         print*,'pressure (presg) = ',p
c         print*,' p,rkvars(numeq) ' ,p,rkvars(numeq)
      else
c         p = p0*exp(-(z-z0)/Hscale0)                       ! 29 Mar 2006
c         print*,'p0 = ',p0
c         print*,'z = ',z
c         print*,'z0 = ',z0
c         print*,'hscale0 = ', hscale0
c         print*,'pressure (presg) = ',p
      endif

      constant = (p*mw)/(rgas*T)                          ! 17 Feb 2006
      ppt = 0.d0                                           ! 17 Feb 2006
      ppr = constant*gr                                   ! 17 Feb 2006
      ppth = constant*gth                                 ! 17 Feb 2006
      ppph = constant*gph                                 ! 17 Feb 2006

      CALL ppres                                                  !  2 Mar 2010

      RETURN                                               ! 25 Jan 2006

      ENTRY presg(pressure)                                ! 24 Jan 2006
      z = R - earthr                                       ! 25 Jan 2006
c      print*,'R =',r
c      print*,'earthr =',earthr
c      print*,'z =',z
      
      CALL tempg(temperature)                              ! 24 Jan 2006
      T = temperature%value

      CALL grav(gravity,coriolis)                                  ! 27 Jun 2016
      g = gravity%gs
      gr = gravity%gsr
      gth = gravity%gsth
      gph = gravity%gsph

      CALL molwtg(molecular_weight)                        ! 24 Jan 2006
      mw = molecular_weight%value

      hscale = (rgas*T)/(mw*g)
c      PRINT *,'hscale pressure = ',hscale

      if (eqnum.gt.0.d0) then
         numeq = eqnum + 0.5                               ! 29 Mar 2006
         p = p0*exp(rkvars(numeq)-(z-z0)/Hscale0)          ! 29 Mar 2006
c         print*,'p0 = ',p0
c         print*,'numeq = ',numeq
c         print*,'z = ',z
c         print*,'z0 = ',z0
c         print*,'hscale0 = ', hscale0
c         print*,'pressure (presg) = ',p
      else
         p = p0*exp(-(z-z0)/Hscale0)                       ! 29 Mar 2006
c         print*,'p0 = ',p0
c         print*,'z = ',z
c         print*,'z0 = ',z0
c         print*,'hscale0 = ', hscale0
c         print*,'pressure (presg) = ',p
      endif

      pressure%value = p                                   ! 25 Jan 2006
      constant = (p*mw)/(rgas*T)                           ! 25 Jan 2006

c      rtest = (mw*gr)/(rgas*T)
c      rtest2 = (((gr/mw)*molecular_weight%pr) -
c     & ((gr/T)*temperature%pr) + (gravity%prr))*mw/(rgas*T)

c      print *,'rtest = ',rtest
c      print *,'rtest2 = ',rtest2
c      print *,' '

      pressure%ptm = 0.d0
      pressure%pr = constant*gr
      pressure%pth = constant*gth
      pressure%pph = constant*gph

c      print*,'pressure (pscale) =',pressure%value
c      print*,'pressure gradient (pscale) =',pressure%pr

      ppt = pressure%ptm                                   ! 17 Feb 2006
      ppr = pressure%pr                                    ! 17 Feb 2006
      ppth = pressure%pth                                  ! 17 Feb 2006
      ppph = pressure%pph                                  ! 17 Feb 2006

      pressure%ptmtm = 0.d0
      pressure%prtm = constant*(((gr/mw)*molecular_weight%ptm) -
     & ((gr/T)*temperature%ptm) + (gravity%prtm)) +
     & (pressure%pr*pressure%ptm/pressure%value)
      pressure%pthtm = constant*(((gth/mw)*molecular_weight%ptm) -
     & ((gth/T)*temperature%ptm) + (gravity%pthtm)) +
     & (pressure%pth*pressure%ptm/pressure%value)
      pressure%pphtm = constant*(((gph/mw)*molecular_weight%ptm) -
     & ((gph/T)*temperature%ptm) + (gravity%pphtm)) +
     & (pressure%pph*pressure%ptm/pressure%value)

      pressure%ptmr = 0.d0
      pressure%prr = constant*(((gr/mw)*molecular_weight%pr) -
     & ((gr/T)*temperature%pr) + (gravity%prr*molecular_weight%pr)) +
     & (pressure%pr**2/pressure%value)
      pressure%pthr = constant*(((gth/mw)*molecular_weight%pr) -
     & ((gth/T)*temperature%pr) + (gravity%pthr)) +
     & (pressure%pth*pressure%pr/pressure%value)
      pressure%pphr = constant*(((gph/mw)*molecular_weight%pr) -
     & ((gph/T)*temperature%pr) + (gravity%pphr)) +
     & (pressure%pph*pressure%pr/pressure%value)

      pressure%ptmth = 0.d0
      pressure%prth = constant*(((gr/mw)*molecular_weight%pth) -
     & ((gr/T)*temperature%pth) + (gravity%prth)) +
     & (pressure%pr*pressure%pth/pressure%value)
      pressure%pthth = constant*(((gth/mw)*molecular_weight%pth) -
     & ((gth/T)*temperature%pth) + (gravity%pthth)) +
     & (pressure%pth**2/pressure%value)
      pressure%pphth = constant*(((gph/mw)*molecular_weight%pth) -
     & ((gph/T)*temperature%pth) + (gravity%pphth)) +
     & (pressure%pph*pressure%pth/pressure%value)

      pressure%ptmph = 0.d0
      pressure%prph = constant*(((gr/mw)*molecular_weight%pph) -
     & ((gr/T)*temperature%pph) + (gravity%prph)) +
     & (pressure%pr*pressure%pph/pressure%value)
      pressure%pthph = constant*(((gth/mw)*molecular_weight%pph) -
     & ((gth/T)*temperature%pph) + (gravity%pthph)) +
     & (pressure%pth*pressure%pph/pressure%value)
      pressure%pphph = constant*(((gph/mw)*molecular_weight%pph) -
     & ((gph/T)*temperature%pph) + (gravity%pphph)) +
     & (pressure%pph**2/pressure%value)
 
c      print *,'pressure%pr = ',pressure%pr
c      print *,'pressure%prr = ',pressure%prr
c      print *,'molecular_weight%pr = ',molecular_weight%pr
c      print *,'gravity%pr = ',gravity%pr
c      print *,'temperature%pr = ',temperature%pr
c
c      print *,'pressure%pth = ',pressure%pth
c      print *,'molecular_weight%pth = ',molecular_weight%pth
c      print *,'gravity%pth = ',gravity%pth
c      print *,'temperature%pth = ',temperature%pth
c      print *,'mw = ',mw
c      print *,'gr = ',gr
c      print *,'gth = ',gth
c      print *,'gph = ',gph
c      print *,'T = ',T

      CALL ppresg(pressure) ! 24 January 2006

      RETURN
C ---------------------------------------------------------------------
      END SUBROUTINE PscaleH
C ---------------------------------------------------------------------

C ---------------------------------------------------------------------
      SUBROUTINE setpres
C ---------------------------------------------------------------------
C     This subroutine is called once.

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

C     COMMON DECK "CB19" INSERTED HERE
      INTEGER PRMX,PRNTBL,PRITBL,PRFRMTB,IDSPR(10)
      COMMON/B19/PRMX,PRNTBL(10),PRITBL(10),PRFRMTB(10),PRGP(11)
      EQUIVALENCE (PRGP,IDSPR),(ANP,PRGP(11))
      INTEGER PFMX,PFNTBL(10),PFITBL(10),PFFRMTB(10)
      DATA PFMX/1/
      DATA PFNTBL/1,11,8*0/
      DATA PFITBL/1,9*0/
      DATA PFFRMTB/1,9*0/
      DATA RECOGP/1.d0/

      PRMX=PFMX
      CALL imove(PRNTBL,PFNTBL,10)  
      CALL imove(PRITBL,PFITBL,10)  
      CALL imove(PRFRMTB,PFFRMTB,10)

      CALL setppres ! 2 Feb 2006
      RETURN                        
C ---------------------------------------------------------------------
      END SUBROUTINE setpres
C ---------------------------------------------------------------------

C ---------------------------------------------------------------------
      SUBROUTINE ipres ! 24 January 2006
C ---------------------------------------------------------------------
C     This subroutine is called once for each new W array.
      IMPLICIT double precision (a-h,o-z)

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

      PARAMETER (NWARSZ=1000)
C     COMMON DECK "WW" INSERTED HERE
      COMMON/WW/dum(10),MAXW,W(NWARSZ)

C     PRESSURE    550-574
      EQUIVALENCE (W(550),PMODEL),(W(551),PFORM),(W(552),PID)

C     COMMON DECK "PP" INSERTED HERE
      REAL(KIND=8) p,modp
      COMMON/PP/MODP(4),p,ppt,ppr,ppth,ppph

C     COMMON DECK "B19" INSERTED HERE
      INTEGER PRMX,PRNTBL,PRITBL,PRFRMTB,IDSPR(10)
      COMMON/B19/PRMX,PRNTBL(10),PRITBL(10),PRFRMTB(10),PRGP(11)
      EQUIVALENCE (PRGP,IDSPR),(ANP,PRGP(11))
      INTEGER PFMX,PFNTBL(10),PFITBL(10),PFFRMTB(10)
      DATA PFMX/1/
      DATA PFNTBL/1,11,8*0/
      DATA PFITBL/1,9*0/
      DATA PFFRMTB/1,9*0/
      DATA RECOGP/1.d0/

      IF (pmodel /= 2.d0) THEN
           PRINT *,' Model check number mismatch: '
           STOP' W(550) should be 2 if you are using model PscaleH.'!25 Jan 2006
      ENDIF

      modp(1) = 'PscaleH' ! 25 Jan 2006
      modp(2) = pid ! 24 January 2006

      CALL ippres ! 24 January 2006
      RETURN

      ENTRY ipresg(modpg) ! 24 January 2006
C     This subroutine is called once for each new W array.

      IF (pmodel /= 2.d0) THEN
           PRINT *,' Model check number mismatch: '
           STOP' W(550) should be 2 if you are using model PscaleH.'!25 Jan 2006
      ENDIF

      modpg%model_name = 'PscaleH' ! 25 Jan 2006
      modpg%data_id = pid
      modpg%pert_name = ''
      modpg%pert_id = 0.d0

!      modp(1) = modpg%model_name                    ! 2 Feb 2006 ! 30 Apr 2010
!      modp(2) = modpg%data_id                       ! 2 Feb 2006 ! 30 Apr 2010
      modp(1) = 'PscaleH'                                         ! 30 Apr 2010
      modp(2) = pid                                               ! 30 Apr 2010

      CALL ippresg(modpg) ! 24 January 2006

      RETURN
C ---------------------------------------------------------------------
      END SUBROUTINE ipres
C ---------------------------------------------------------------------
