C --------------------------------------------------------------------- C PscaleH.for C --------------------------------------------------------------------- ! 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 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) 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) 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 ---------------------------------------------------------------------