forked from 3rdparty/wrf-python
You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
264 lines
7.5 KiB
264 lines
7.5 KiB
c====================================================================== |
|
c |
|
c !IROUTINE: WETBULBCALC -- Calculate wet bulb temperature (C) |
|
c |
|
c !DESCRIPTION: |
|
c |
|
c Calculates wet bulb temperature in C, given pressure in |
|
c temperature in K and mixing ratio in kg/kg. |
|
c |
|
c !INPUT: |
|
c nx - index for x dimension |
|
c ny - index for y dimension |
|
c nz - index for z dimension |
|
c prs - pressure (mb) |
|
c tmk - temperature (K) |
|
c qvp - water vapor mixing ratio (kg/kg) |
|
c |
|
c !OUTPUT: |
|
c twb - Wet bulb temperature (C) |
|
c |
|
c !ASSUMPTIONS: |
|
c |
|
c !REVISION HISTORY: |
|
c 2009-March - Mark T. Stoelinga - from RIP4.5 |
|
c 2010-August - J. Schramm |
|
c 2014-March - A. Jaye - modified to run with NCL and ARW wrf output |
|
c |
|
c !INTERFACE: |
|
c ------------------------------------------------------------------ |
|
C NCLFORTSTART |
|
subroutine wetbulbcalc(prs,tmk,qvp,twb,nx,ny,nz,psafile) |
|
implicit none |
|
integer nx, ny, nz |
|
double precision prs(nz,ny,nx) |
|
double precision tmk(nz,ny,nx) |
|
double precision qvp(nz,ny,nx) |
|
double precision twb(nz,ny,nx) |
|
character*(*) psafile |
|
C NCLEND |
|
integer i,j,k |
|
integer jtch,jt,ipch,ip |
|
double precision q, t, p, e, tlcl, eth |
|
double precision fracip,fracip2,fracjt,fracjt2 |
|
double precision PSADITHTE(150),PSADIPRS(150),PSADITMK(150,150) |
|
double precision tonpsadiabat |
|
double precision eps,tlclc1,tlclc2,tlclc3,tlclc4,gamma |
|
double precision gammamd,thtecon1,thtecon2,thtecon3,celkel |
|
double precision rgas,rgasmd,cp,cpmd |
|
|
|
c |
|
c Before looping, set lookup table for getting temperature on |
|
c a pseudoadiabat. |
|
c |
|
CALL DLOOKUP_TABLE(PSADITHTE,PSADIPRS,PSADITMK,psafile) |
|
|
|
c Define constants |
|
|
|
rgas=287.04 !J/K/kg |
|
rgasmd=.608 ! rgas_moist=rgas*(1.+rgasmd*qvp) |
|
cp=1004. ! J/K/kg Note: not using Bolton's value of 1005.7 |
|
cpmd=.887 ! cp_moist=cp*(1.+cpmd*qvp) |
|
eps=0.622 |
|
tlclc1=2840. |
|
tlclc2=3.5 |
|
tlclc3=4.805 |
|
tlclc4=55. |
|
gamma=rgas/cp |
|
gammamd=rgasmd-cpmd ! gamma_moist=gamma*(1.+gammamd*qvp) |
|
thtecon1=3376. ! K |
|
thtecon2=2.54 |
|
thtecon3=.81 |
|
celkel=273.15 |
|
|
|
DO k=1,nx |
|
DO j=1,ny |
|
DO i=1,nz |
|
q=dmax1(qvp(i,j,k),1.d-15) |
|
t=tmk(i,j,k) |
|
p=prs(i,j,k)/100. |
|
e=q*p/(eps+q) |
|
tlcl=tlclc1/(dlog(t**tlclc2/e)-tlclc3)+tlclc4 |
|
eth=t*(1000./p)**(gamma*(1.+gammamd*q))* |
|
& exp((thtecon1/tlcl-thtecon2)*q*(1.+thtecon3*q)) |
|
|
|
|
|
c |
|
c Now we need to find the temperature (in K) on a moist adiabat |
|
c (specified by eth in K) given pressure in hPa. It uses a |
|
c lookup table, with data that was generated by the Bolton (1980) |
|
c formula for theta_e. |
|
c |
|
c First check if pressure is less than min pressure in lookup table. |
|
c If it is, assume parcel is so dry that the given theta-e value can |
|
c be interpretted as theta, and get temperature from the simple dry |
|
c theta formula. |
|
c |
|
|
|
if (p.le.psadiprs(150)) then |
|
tonpsadiabat=eth*(p/1000.)**gamma |
|
else |
|
c |
|
c Otherwise, look for the given thte/prs point in the lookup table. |
|
c |
|
do jtch=1,150-1 |
|
if (eth.ge.psadithte(jtch).and.eth.lt. |
|
& psadithte(jtch+1)) then |
|
jt=jtch |
|
goto 213 |
|
endif |
|
enddo |
|
jt=-1 |
|
213 continue |
|
do ipch=1,150-1 |
|
if (p.le.psadiprs(ipch).and.p.gt.psadiprs(ipch+1)) then |
|
ip=ipch |
|
goto 215 |
|
endif |
|
enddo |
|
ip=-1 |
|
215 continue |
|
if (jt.eq.-1.or.ip.eq.-1) then |
|
print*, |
|
& 'Outside of lookup table bounds. prs,thte=',p,eth |
|
stop |
|
endif |
|
fracjt=(eth-psadithte(jt))/(psadithte(jt+1)-psadithte(jt)) |
|
fracjt2=1.-fracjt |
|
fracip=(psadiprs(ip)-p)/(psadiprs(ip)-psadiprs(ip+1)) |
|
fracip2=1.-fracip |
|
if (psaditmk(ip,jt).gt.1e9.or.psaditmk(ip+1,jt).gt.1e9.or. |
|
& psaditmk(ip,jt+1).gt.1e9.or. |
|
& psaditmk(ip+1,jt+1).gt.1e9) then |
|
print*, |
|
& 'Tried to access missing tmperature in lookup table.' |
|
print*, |
|
& 'Prs and Thte probably unreasonable. prs,thte=' |
|
& ,p,eth |
|
stop |
|
endif |
|
tonpsadiabat=fracip2*fracjt2*psaditmk(ip ,jt )+ |
|
& fracip *fracjt2*psaditmk(ip+1,jt )+ |
|
& fracip2*fracjt *psaditmk(ip ,jt+1)+ |
|
& fracip *fracjt *psaditmk(ip+1,jt+1) |
|
endif |
|
|
|
twb(i,j,k)=tonpsadiabat |
|
|
|
ENDDO |
|
ENDDO |
|
ENDDO |
|
|
|
c |
|
return |
|
end |
|
c====================================================================== |
|
c |
|
c !IROUTINE: omgcalc -- Calculate omega (dp/dt) |
|
c |
|
c !DESCRIPTION: |
|
c |
|
c Calculate approximate omega, based on vertical velocity w (dz/dt). |
|
c It is approximate because it cannot take into account the vertical |
|
c motion of pressure surfaces. |
|
c |
|
c !INPUT: |
|
c mx - index for x dimension |
|
c my - index for y dimension |
|
c mx - index for vertical dimension |
|
c qvp - water vapor mixing ratio (kg/kg) |
|
c tmk - temperature (K) |
|
c www - vertical velocity (m/s) |
|
c prs - pressure (Pa) |
|
c |
|
c !OUTPUT: |
|
c omg - omega (Pa/sec) |
|
c |
|
c !ASSUMPTIONS: |
|
c |
|
c !REVISION HISTORY: |
|
c 2009-March - Mark T. Stoelinga - from RIP4.5 |
|
c 2010-August - J. Schramm |
|
c 2014-March - A. Jaye - modified to run with NCL and ARW wrf output |
|
c |
|
c ------------------------------------------------------------------ |
|
c NCLFORTSTART |
|
subroutine omgcalc(qvp,tmk,www,prs,omg,mx,my,mz) |
|
implicit none |
|
integer mx, my, mz |
|
double precision qvp(mz,my,mx) |
|
double precision tmk(mz,my,mx) |
|
double precision www(mz,my,mx) |
|
double precision prs(mz,my,mx) |
|
double precision omg(mz,my,mx) |
|
c NCLEND |
|
c Local variables |
|
integer i, j, k |
|
double precision grav,rgas,eps |
|
c |
|
c Constants |
|
c |
|
grav=9.81 ! m/s**2 |
|
rgas=287.04 !J/K/kg |
|
eps=0.622 |
|
|
|
do k=1,mx |
|
do j=1,my |
|
do i=1,mz |
|
omg(i,j,k)=-grav*prs(i,j,k)/ |
|
& (rgas*((tmk(i,j,k)*(eps+qvp(i,j,k)))/ |
|
& (eps*(1.+qvp(i,j,k)))))*www(i,j,k) |
|
enddo |
|
enddo |
|
enddo |
|
c |
|
return |
|
end |
|
c====================================================================== |
|
c |
|
c !IROUTINE: VIRTUAL_TEMP -- Calculate virtual temperature (K) |
|
c |
|
c !DESCRIPTION: |
|
c |
|
c Calculates virtual temperature in K, given temperature |
|
c in K and mixing ratio in kg/kg. |
|
c |
|
c !INPUT: |
|
c NX - index for x dimension |
|
c NY - index for y dimension |
|
c NZ - index for z dimension |
|
c RATMIX - water vapor mixing ratio (kg/kg) |
|
c TEMP - temperature (K) |
|
c |
|
c !OUTPUT: |
|
c TV - Virtual temperature (K) |
|
c |
|
c !ASSUMPTIONS: |
|
c |
|
c !REVISION HISTORY: |
|
c 2009-March - Mark T. Stoelinga - from RIP4.5 |
|
c 2010-August - J. Schramm |
|
c 2014-March - A. Jaye - modified to run with NCL and ARW wrf output |
|
c |
|
c ------------------------------------------------------------------ |
|
C NCLFORTSTART |
|
SUBROUTINE VIRTUAL_TEMP(TEMP,RATMIX,TV,NX,NY,NZ) |
|
IMPLICIT NONE |
|
INTEGER NX,NY,NZ |
|
DOUBLE PRECISION TEMP(NZ,NY,NX) |
|
DOUBLE PRECISION RATMIX(NZ,NY,NX) |
|
DOUBLE PRECISION TV(NZ,NY,NX) |
|
C NCLEND |
|
INTEGER I,J,K |
|
DOUBLE PRECISION EPS |
|
EPS = 0.622D0 |
|
DO K=1,NX |
|
DO J=1,NY |
|
DO I=1,NZ |
|
TV(I,J,K) = TEMP(I,J,K)* (EPS+RATMIX(I,J,K))/ |
|
& (EPS* (1.D0+RATMIX(I,J,K))) |
|
ENDDO |
|
ENDDO |
|
ENDDO |
|
RETURN |
|
END
|
|
|