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