!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CPTC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE cptc(nx,ny,nz,ibgn,iend,jbgn,jend, & 3
zp,roufns,wspd,ptsfc,pt1, c_ptneu,c_pt,c_u)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Compute c_pt (a nondimensional temperature scale)
!
! The quantity c_pt is used by the subroutine SFCFLX to obtain
! surface fluxes for both the unstable and stable cases.
!
!-----------------------------------------------------------------------
!
! AUTHORS: V. Wong, X. Song and N. Lin
! 9/10/1993
!
! For the unstable case (Bulk Richardson number bulkri < 0 ), the
! formulation is based on the paper "On the Analytical Solutions of
! Flux-Profile relationships for the Atmospheric Surface Layer" by
! D.W. Byun in J. of Applied Meteorlogy, July 1990, pp. 652-657.
!
! For stable case, the formulation is based on the paper "A Short
! History of the Operational PBL - Parameterization at ECMWF" by
! J.F.Louis, M. Tiedtke and J.F. Geleyn in "Workshop on Planetary
! Boundary Layer Parameterization", a publication by European Centre
! for Medium Range Weather Forecasts, 25-27 Nov. 1981.
!
! MODIFICATION HISTORY:
!
! 9/04/1994 (K. Brewster, V. Wong and X. Song)
! Facelift.
!
! 2/27/95 (V. Wong and X. Song)
!
! 02/07/96 (V.Wong and X.Song)
! Set an upper limiter, z1limit, for depth of the surface layer z1.
!
! 05/01/97 (V. Wong and X. Tan)
! Changed the computation of temperature scale
!
! 05/29/97 (V. Wong and X. Tan)
! Modified the formulation considering the height of the surface
! layer z1 may equal zero.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! ibgn i-index where evaluation begins.
! iend i-index where evaluation ends.
! jbgn j-index where evaluation begins.
! jend j-index where evaluation ends.
!
! zp The physical height coordinate defined at w-point of
! staggered grid.
! wspd Surface wind speed (m/s), defined as
! sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
! roufns Surface roughness
!
! ptsfc Potential temperature at the ground level (K)
! pt1 Potential temperature at the 1st scalar point above
! ground level, (K)
!
! c_ptneu Temperature scale (K) at neutral state, defined by
! surface heat flux / friction velocity
!
! OUTPUT:
!
! c_pt Temperature scale (K), defined by
! surface heat flux / friction velocity
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! The number grid points in 3 directions
INTEGER :: ibgn ! i-index where evaluation begins
INTEGER :: iend ! i-index where evaluation ends
INTEGER :: jbgn ! j-index where evaluation begins
INTEGER :: jend ! j-index where evaluation ends
REAL :: zp (nx,ny,nz) ! The physical height coordinate
! defined at w-point of staggered grid.
REAL :: roufns(nx,ny) ! Surface roughness length
REAL :: wspd (nx,ny) ! Surface wind speed (m/s)
REAL :: ptsfc(nx,ny) ! Potential temperature at the
! ground level (K)
REAL :: pt1 (nx,ny) ! Potential temperature at the
! 1st scalar
! point above ground level, (K)
REAL :: c_ptneu(nx,ny) ! Normalized temperature scale (K)
! at neutral state
REAL :: c_pt (nx,ny) ! Temperature scale (K)
REAL :: c_u (nx,ny) ! Friction velocity (m/s)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j
REAL :: z1 ! The height of 1st scalar point above
! ground level (m)
REAL :: dz0 ! z1-roufns, where roufns is defined as
! surface roughness length
REAL :: bulkri ! Bulk Richardson number
REAL :: stabp ! Monin-Obukhov STABility Parameter
! (zeta)
REAL :: y1,y0,psih ! Intermediate parameters needed
REAL :: z1drou,qb3pb2
REAL :: c7,c8
REAL :: sb,qb,pb,thetab,tb ! During computations
REAL :: a,b,c,d
REAL :: tempan,sqrtqb
REAL :: zt ! Thermal roughness length
REAL :: dzt,ztdrou
REAL :: z1droup,ztdroup
REAL, PARAMETER :: epsilon = 1.0E-6
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
INCLUDE 'globcst.inc'
INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Following Byun (1990).
!
!-----------------------------------------------------------------------
!
DO j=jbgn,jend
DO i=ibgn,iend
z1 = 0.5*(zp(i,j,3)-zp(i,j,2))
z1 = MIN(z1,z1limit)
! zt = ztz0*roufns(i,j)
!---------------------------------------------------------------------
! Test of new thermal roughness equation (Chen/Dudhia 2001) (JAB)
!---------------------------------------------------------------------
zt = (0.4*c_u(i,j)/0.000024) + 100.0
IF (abs(zt) > epsilon ) THEN
zt = 1.0/zt
ELSE IF (abs(zt) <= epsilon ) THEN
zt = ztz0*roufns(i,j)
END IF
!----------------------------------------------------------------------
dzt = z1-zt
ztdrou = z1/zt
ztdroup = (z1+zt)/zt
dz0 = z1-roufns(i,j)
z1drou = z1/roufns(i,j)
z1droup = (z1+roufns(i,j))/roufns(i,j)
bulkri = (g/ptsfc(i,j))*(pt1(i,j)-ptsfc(i,j))*dz0/ &
(wspd(i,j)*wspd(i,j))
!-------------------------------------------------------------------
! Soilmodel_option = 2 (Implicit)
!--------------------------------------------------------------------
! IF (soilmodel_option == 2) THEN
!
! IF (bulkri <= 0.0) THEN
!---------------------------------------------------------------------
! Unstable case
!---------------------------------------------------------------------
!
! stabp = 1.0 - ( (15.0 * bulkri) / &
! ( 1.0 + 7.5 * (10.0*kv*kv/(LOG(z1drou)*LOG(ztdrou))) * &
! sqrt(-bulkri * z1drou) ) )
!
! ELSE
!---------------------------------------------------------------------
! Stable case
!---------------------------------------------------------------------
!
! stabp = EXP(-bulkri)
!
! END IF
!
! c_pt(i,j) = kv*kv*stabp / (LOG(z1drou)*LOG(ztdrou))
!
! END IF !soilmodel_option = 2
!-------------------------------------------------------------------
!-------------------------------------------------------------------
! Soilmodel_option = 1 (Force-Restore)
!------------------------------------------------------------------
! IF (soilmodel_option == 1) THEN
IF (bulkri <= 0.0) THEN
!
!-----------------------------------------------------------------------
!
! Unstable case: See equations (28)-(34) in Byun (1990).
!
!-----------------------------------------------------------------------
!
bulkri = MAX (bulkri,-10.0)
sb =bulkri/prantl0l
qb=oned9*(c1l+c2l*sb*sb)
pb=oned54*(c3l+c4l*sb*sb)
qb3pb2=qb**3-pb*pb
c7 = (z1*dzt*LOG(z1droup)*LOG(z1droup))/(dz0*dz0*LOG(ztdroup))
IF( qb3pb2 >= 0.0 ) THEN
sqrtqb = SQRT(qb)
tempan = MAX( -1.0, MIN( 1.0, pb/(sqrtqb**3) ) )
thetab=ACOS(tempan)
stabp =c7*(-2.0*SQRT(qb)*COS(thetab/3.0)+c5l)
ELSE
tb =(SQRT(-qb3pb2)+ABS(pb))**oned3
stabp =c7*(-(tb+qb/tb)+c5l)
END IF
!
!-----------------------------------------------------------------------
!
! According to equation (15) in Byun (1990).
!
!-----------------------------------------------------------------------
!
c8=gammahl*stabp
y1=SQRT(1.0 - c8)
y0=SQRT(1.0 - c8/ztdrou)
psih=2.0*LOG((y1+1.0)/(y0+1.0))
!
!-----------------------------------------------------------------------
!
! Compute c_pt via equation (11) in Byun (1981).
!
!-----------------------------------------------------------------------
!
c_pt(i,j)=kv / (prantl0l*(LOG(ztdroup)-psih))
ELSE
!
!-----------------------------------------------------------------------
!
! Stable case: See Louis et al (1981).
!
!-----------------------------------------------------------------------
!
a=kv/LOG(ztdroup)
b=5.0
d=5.0
c=SQRT(1.0+d*bulkri)
c_pt(i,j) = SQRT(1.0+2.0*b*bulkri/c)
c_pt(i,j) = a*c_pt(i,j)/(prantl0l*(1.0+3.0*b*bulkri*c))
END IF
! END IF
END DO
END DO
RETURN
END SUBROUTINE cptc
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CPTCWTR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE cptcwtr(nx,ny,nz,ibgn,iend,jbgn,jend,zp,soiltyp,wspd, & 3
ptsfc,pt1,c_ptwtrneu,c_ptwtr)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Compute c_ptwtr (the product of ustar and ptstar) at sea case.
! Note: a temperature scale defined as surface heat flux divided
! by the friction velocity.
!
! The quantity c_ptwtr is used by the subroutine SFCFLX to obtain
! surface fluxes for both the unstable and stable cases.
!
!-----------------------------------------------------------------------
!
! AUTHORS: V. Wong and X. Song
! 8/04/1994
!
! MODIFICATION HISTORY:
!
! 2/27/95 (V.W. and X.S.)
!
! 1/12/96 (V.W. and X.S.)
! Changed the calculation related to zo over the sea.
! Added kvwtr to denote the Von Karman constant over the sea;
! Set a lower limiter for zo, zolimit, and an upper limiter for z1, z1limit.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! iend i-index where evaluation ends.
! jend j-index where evaluation ends.
!
! zp The physical height coordinate defined at w-point of
! staggered grid.
! wspd Surface wind speed (m/s), defined as
! sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
! ptsfc Potential temperature at the ground level (K)
! pt1 Potential temperature at the 1st scalar point above
! ground level, (K)
!
!
! OUTPUT:
!
! c_ptwtr The product of ustar and ptstar. ptstar is temperature
! scale (K), defined by surface heat flux / friction velocity
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! The number grid points in 3 directions
INTEGER :: ibgn ! i-index where evaluation begins.
INTEGER :: iend ! i-index where evaluation ends.
INTEGER :: jbgn ! j-index where evaluation begins.
INTEGER :: jend ! j-index where evaluation ends.
REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at
! w-point of staggered grid.
INTEGER :: soiltyp(nx,ny) ! Soil type at each point.
REAL :: wspd (nx,ny) ! Surface wind speed (m/s)
REAL :: ptsfc(nx,ny) ! Potential temperature at the ground level
! (K)
REAL :: pt1 (nx,ny) ! Potential temperature at the 1st scalar
! point above ground level, (K)
REAL :: c_ptwtrneu(nx,ny) ! Product of ustar and ptstar
REAL :: c_ptwtr(nx,ny) ! Product of ustar and ptstar
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j
REAL :: z1 ! The height of 1st scalar point above
! ground level (m)
REAL :: zo ! Defined as surface momentum roughness
! length
REAL :: zt ! Defined as surface heat transfer
! roughness length
REAL :: dzo ! z1-zo
REAL :: dzt ! z1-zt
REAL :: z1dzo ! z1/zo
REAL :: z1dzt ! z1/zt
REAL :: xcdn ! Cdn (sea)
REAL :: xcdh ! Hot-wired value of Cdh (sea)
REAL :: bulkri ! Bulk Richardson number
REAL :: stabp ! Monin-Obukhov STABility Parameter
! (zeta)
REAL :: y1,y0,psih ! Intermediate parameters needed
REAL :: qb3pb2
REAL :: c7,c8
REAL :: sb,qb,pb,thetab,tb ! during computations
REAL :: a,b,c,d
REAL :: tempan,sqrtqb
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
INCLUDE 'globcst.inc'
INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
xcdh = 1.1E-3
DO j=jbgn,jend
DO i=ibgn,iend
IF ( soiltyp(i,j) == 13) THEN
xcdn = (0.4+0.07*wspd(i,j)) * 1.e-3
z1 = 0.5*(zp(i,j,3)-zp(i,j,2))
z1 = MIN(z1,z1limit)
!
! calculate zo and zt
!
zo =0.032*xcdn*wspd(i,j)*wspd(i,j)/9.8
zo = MAX(zo,zolimit)
zt = z1 * EXP( -kv*SQRT(xcdn)/(prantl0w*xcdh) )
dzo = z1-zo
z1dzo = z1/zo
dzt = z1-zt
z1dzt = z1/zt
bulkri = (g/ptsfc(i,j))*(pt1(i,j)-ptsfc(i,j))*dzo*dzo/ &
(dzt*wspd(i,j)*wspd(i,j))
IF (bulkri <= 0.0) THEN
!
!-----------------------------------------------------------------------
!
! Unstable case: A modified formulation, which is similar to
! equations (28)-(34) in Byun (1990).
!
!-----------------------------------------------------------------------
!
bulkri = MAX (bulkri,-10.0)
sb =bulkri/prantl0w
qb=oned9*(c1w+c2w*sb*sb)
pb=oned54*(c3w+c4w*sb*sb)
qb3pb2=qb**3-pb*pb
c7 = (z1*dzt/(dzo*dzo))*(ALOG(z1dzo)*ALOG(z1dzo)/ALOG(z1dzt))
IF( qb3pb2 >= 0.0 ) THEN
sqrtqb = SQRT(qb)
tempan = MAX( -1.0, MIN( 1.0, pb/(sqrtqb**3) ) )
thetab=ACOS(tempan)
stabp =c7*(-2.0*SQRT(qb)*COS(thetab/3.0)+c5w)
ELSE
tb =(SQRT(-qb3pb2)+ABS(pb))**oned3
stabp =c7*(-(tb+qb/tb)+c5w)
END IF
!
!-----------------------------------------------------------------------
!
! According to a modified equation, which is similar to equation (14)
! in Byun (1990).
!
!-----------------------------------------------------------------------
!
c8=gammamw*stabp
y1=SQRT(1.0 - c8)
y0=SQRT(1.0 - c8/z1dzt)
psih=2.0*ALOG((y1+1.0)/(y0+1.0))
!
!-----------------------------------------------------------------------
!
! Compute ptstar via equation (11) in Byun (1981).
!
!-----------------------------------------------------------------------
!
c_ptwtr(i,j)=kv / (prantl0w*(ALOG(z1dzt)-psih))
ELSE
!
!-----------------------------------------------------------------------
!
! Stable case: With the modified formulation in Louis et al (1981).
!
!-----------------------------------------------------------------------
!
a=kv*kv/(prantl0w*LOG(z1dzo)*LOG(z1dzt))
b=5.0
d=5.0
c=SQRT(1.0+d*bulkri)
c_ptwtr(i,j) = SQRT(1.0+2.0*b*bulkri/c)
c_ptwtr(i,j) = a*c_ptwtr(i,j)/(prantl0l*(1.0+3.0*b*bulkri*c))
END IF
END IF
END DO
END DO
RETURN
END SUBROUTINE cptcwtr
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE HINTRP1 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE hintrp1(nx,ny,nz, kbgn,kend,a3din,z3d,zlevel, a2dout) 26
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Interpolate a 3-D array to horizontal level z=zlevel.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! Based on original SECTHRZ.
! 12/10/98.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
! kbgn
! kend
!
! a3din 3-d input array
! z3d z-coordinate of data in a3din
! zlevel Level to which data is interpolated.
!
! OUTPUT:
! a2dout 2-d output array interpolated to zlevel
!
!-----------------------------------------------------------------------
!
! Parameters of output
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
INTEGER :: kbgn, kend
REAL :: a3din(nx,ny,nz) ! 3-d input array
REAL :: z3d (nx,ny,nz) ! z-coordinate of data in a3din
REAL :: zlevel ! Level to which data is interpolated.
REAL :: a2dout(nx,ny) ! 2-d output array interpolated to zlevel
INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
! Find index for interpolation
!
!-----------------------------------------------------------------------
DO i=1,nx-1
DO j=1,ny-1
IF(zlevel <= z3d(i,j,kbgn)) GO TO 11
IF(zlevel >= z3d(i,j,kend)) GO TO 12
DO k=kbgn,kend-1
IF(zlevel >= z3d(i,j,k).AND.zlevel < z3d(i,j,k+1)) GO TO 15
END DO
11 k=kbgn
GO TO 15
12 k=kend-1
GO TO 15
15 a2dout(i,j)=a3din(i,j,k)+(a3din(i,j,k+1)-a3din(i,j,k))* &
(zlevel-z3d(i,j,k))/(z3d(i,j,k+1)-z3d(i,j,k))
!-----------------------------------------------------------------------
!
! If the data point is below the ground level, set the
! data value to the missing value.
!
!-----------------------------------------------------------------------
IF( zlevel < z3d(i,j,kbgn) ) a2dout(i,j) = -9999.0
IF( zlevel > z3d(i,j,kend) ) a2dout(i,j) = -9999.0
END DO
END DO
RETURN
END SUBROUTINE hintrp1
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE TEMPER ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE temper ( nx,ny,nz,theta, ppert, pbar, t ) 9
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Using a version of Poisson's formula, calculate temperature.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Joe Bradley
! 12/05/91
!
! MODIFICATIONS:
! Modified by Ming Xue so that arrays are only defined at
! one time level.
! 6/09/92 Added full documentation and phycst include file for
! rddcp=Rd/Cp (K. Brewster)
!
!-----------------------------------------------------------------------
!
! INPUT:
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! theta Potential temperature (degrees Kelvin)
! ppert Perturbation pressure (Pascals)
! pbar Base state pressure (Pascals)
!
! OUTPUT:
!
! t Temperature (degrees Kelvin)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx,ny,nz
!
REAL :: theta(nx,ny,nz) ! potential temperature (degrees Kelvin)
REAL :: ppert(nx,ny,nz) ! perturbation pressure (Pascals)
REAL :: pbar (nx,ny,nz) ! base state pressure (Pascals)
!
REAL :: t (nx,ny,nz) ! temperature (degrees Kelvin)
!
!-----------------------------------------------------------------------
!
! Include file
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Calculate the temperature using Poisson's formula.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
t(i,j,k) = theta(i,j,k) * &
(((ppert(i,j,k) + pbar(i,j,k)) / p0) ** rddcp)
END DO
END DO
END DO
RETURN
END SUBROUTINE temper
SUBROUTINE arps_be1(nx,ny,nz, & 1,11
pres,hgt,tk,wmr, &
lcl,lfc,el,twdf,li,cape,mcape,cin,mcin,tcap, &
p1d,ht1d,t1d,tv1d,td1d,wmr1d, &
partem,buoy,wload,mbuoy,pbesnd,mbesnd,tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the lifting condensation level (lcl), level of free
! convection (lfc), equilibrium level (el), max wet-bulb potential
! temperature difference (twdf), lifted index (LI), Convective
! Available Potential Energy (CAPE), Moist Convective Potential
! Energy (MCAPE, includes water loading), convective inhibition
! (CIN) and lid strength (tcap) over the ARPS grid.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Keith Brewster
! April, 1995 Based on OLAPS, hence LAPS, version of same.
!
! MODIFICATION HISTORY:
! 3/11/1996 (Keith Brewster)
! Cleaned-up, removed OLAPS artifacts.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
!
!-----------------------------------------------------------------------
!
! 3-D input variables
!
!-----------------------------------------------------------------------
!
REAL :: pres(nx,ny,nz)
REAL :: hgt(nx,ny,nz)
REAL :: tk(nx,ny,nz)
REAL :: wmr(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
! Output variables (2d)
!
!-----------------------------------------------------------------------
!
REAL :: lcl(nx,ny)
REAL :: lfc(nx,ny)
REAL :: el(nx,ny)
REAL :: twdf(nx,ny)
REAL :: li(nx,ny)
REAL :: cape(nx,ny)
REAL :: mcape(nx,ny)
REAL :: cin(nx,ny)
REAL :: mcin(nx,ny)
REAL :: tcap(nx,ny)
!
!-----------------------------------------------------------------------
!
! Scratch space for calculations
!
!-----------------------------------------------------------------------
!
REAL :: p1d(nz),ht1d(nz)
REAL :: t1d(nz),tv1d(nz),td1d(nz),wmr1d(nz)
REAL :: partem(nz),buoy(nz),wload(nz)
REAL :: mbuoy(nz),pbesnd(nz),mbesnd(nz)
REAL :: tem1(nx,ny)
!
!-----------------------------------------------------------------------
!
! Functions
!
!-----------------------------------------------------------------------
!
REAL :: wmr2td,tctotv
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,n,imid,jmid,nlevel
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
imid=(nx/2) + 1
jmid=(ny/2) + 1
!
!-----------------------------------------------------------------------
!
! Loop over all columns
!
!-----------------------------------------------------------------------
!
DO j=1,ny-1
DO i=1,nx-1
!
DO k=2,nz-1
n = k-1
p1d(n) = 0.01 * pres(i,j,k)
ht1d(n) = hgt(i,j,k)*1000.
wmr1d(n)= 1000. * wmr(i,j,k)
t1d(n) = tk(i,j,k) - 273.15
td1d(n) = wmr2td(p1d(n),wmr1d(n))
tv1d(n) = tctotv(t1d(n),wmr1d(n))
END DO
!
nlevel = nz-2
!
CALL sindex1
(nz,nlevel,p1d,ht1d,t1d,tv1d,td1d,wmr1d, &
partem,buoy,wload,mbuoy,pbesnd,mbesnd, &
lcl(i,j),lfc(i,j),el(i,j),twdf(i,j), &
li(i,j),cape(i,j),mcape(i,j), &
cin(i,j),mcin(i,j),tcap(i,j))
IF(i == imid .AND. j == jmid) THEN
WRITE(6,*) &
' n p t tpar buoy wmr be'
DO n = 1,nlevel
WRITE(6,801) n,p1d(n),t1d(n),partem(n),buoy(n), &
wmr1d(n),pbesnd(n)
801 FORMAT(1X,i2,f8.1,f10.2,f10.2,f10.4,f10.4,f10.1)
END DO
END IF
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Smooth the output arrays
!
!-----------------------------------------------------------------------
!
CALL smooth9p( lcl,nx,ny,1,nx-1,1,ny-1,0,tem1)
CALL smooth9p( lfc,nx,ny,1,nx-1,1,ny-1,0,tem1)
CALL smooth9p( el,nx,ny,1,nx-1,1,ny-1,0,tem1)
CALL smooth9p( twdf,nx,ny,1,nx-1,1,ny-1,0,tem1)
CALL smooth9p( li,nx,ny,1,nx-1,1,ny-1,0,tem1)
CALL smooth9p( cape,nx,ny,1,nx-1,1,ny-1,0,tem1)
CALL smooth9p(mcape,nx,ny,1,nx-1,1,ny-1,0,tem1)
CALL smooth9p( cin,nx,ny,1,nx-1,1,ny-1,0,tem1)
CALL smooth9p( mcin,nx,ny,1,nx-1,1,ny-1,0,tem1)
CALL smooth9p( tcap,nx,ny,1,nx-1,1,ny-1,0,tem1)
!
RETURN
END SUBROUTINE arps_be1
SUBROUTINE potbe1(maxlev,nlevel,pmean,tmean,wmean, k0, & 2,3
blthte,plcl,tlcl,lcl, &
p,ht,t,tv,td,w, &
partem,buoy,wload,mbuoy,pbesnd,mbesnd, &
pbeneg,velneg,pos_max, &
mbeneg,mvelneg,mpos_max,lfc,el,tcap)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate Convective Available Potential Energy (CAPE)
! in each column of ARPS grid.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Keith Brewster
! February, 1994 Based on OLAPS, hence LAPS, version of same.
! from FSL, by Steve Albers 1991
!
! MODIFICATION HISTORY:
! 3/11/1996 Cleaned-up, removed OLAPS artifacts.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: maxlev,nlevel
REAL :: pmean,tmean,wmean,blthte,plcl,tlcl,lcl
REAL :: p(maxlev),ht(maxlev)
REAL :: t(maxlev),tv(maxlev),td(maxlev),w(maxlev)
REAL :: partem(maxlev),buoy(maxlev),wload(maxlev),mbuoy(maxlev)
REAL :: pbesnd(maxlev),mbesnd(maxlev)
REAL :: pbeneg,velneg,pos_max
REAL :: mbeneg,mvelneg,mpos_max,lfc,el,tcap
!
! Parameters
!
REAL :: g,gamma
PARAMETER (g=9.80665, &
gamma = .009760) ! Dry Adiabatic Lapse Rate Deg/m
!
! Functions
!
REAL :: tsa_fast,tctotv
!
! Misc internal variables
!
INTEGER :: n,nel,k0
REAL :: deltah,delta_ht_dry,delta_ht_wet
REAL :: sntlcl,buoy_lcl,wsat,partv
REAL :: nbe_min,pbe_wet,pbe_dry,pos_area
REAL :: wlow,htzero,adjeng
!
!-----------------------------------------------------------------------
!
! Function f_mrsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_mrsat
!fpp$ expand (f_mrsat)
!!dir$ inline always f_mrsat
!*$* inline routine (f_mrsat)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
! Reset output variables.
!
! These should be the same as what is assigned
! no positive area found. They are limited in
! range to allow for contouring when positive
! areas exist in some columns of a domain and
! not in others.
!
pbeneg=-400.
velneg=20.
pos_max=0.
mbeneg=-400.
mvelneg=20.
mpos_max=0.
lfc=10000.
el=0.
tcap=0.
!
! Inxtialize parcel path arrays
!
partem(k0) = t(k0)
buoy(k0) = 0.
wload(k0) = 0.
mbuoy(k0) = 0.
pbesnd(k0) = 0.
mbesnd(k0) = 0.
! WRITE(6,810)pmean,tmean,wmean,plcl,tlcl,lcl
! 810 format(' pmean,tmean,wmean,plcl,tlcl,lcl',2F10.2,F10.5,2F10.2
! + ,F5.1)
DO n=k0+1,nlevel
deltah = ht(n) - ht(n-1)
IF(plcl < p(n-1))THEN ! lower level is below LCL
IF(plcl < p(n))THEN ! upper level is below LCL
! WRITE(6,*)' DRY CASE'
partem(n)=partem(n-1)-gamma*deltah
partv=tctotv(partem(n),w(k0))
buoy(n)=(partv-tv(n))/tv(n)
pbesnd(n)=pbesnd(n-1)+g*0.5*(buoy(n)+buoy(n-1))*deltah
wload(n)=0.
mbuoy(n)=buoy(n)
mbesnd(n)=pbesnd(n)
IF((p(k0)-p(n)) < 300.) tcap=AMAX1(tcap,(tv(n)-partv))
ELSE ! Upper level is above LCL
!
! BRACKETING CASE AROUND lcl - DRY ADIABATIC PART
!
! WRITE(6,*)' DRY ADIABATIC PART'
delta_ht_dry = lcl - ht(n-1)
! WRITE(6,307)tlcl
!307 format(' PARCEL TEMP AT lcl= ',F10.3)
CALL intrpr(maxlev,nlevel,p,tv,plcl,sntlcl)
partv=tctotv(tlcl,w(k0))
buoy_lcl=(partv-sntlcl)/sntlcl
pbe_dry=g*0.5*(buoy_lcl+buoy(n-1))*delta_ht_dry
IF((p(k0)-plcl) < 300.) tcap=AMAX1(tcap,(sntlcl-partv))
! WRITE(6,777)N,P(N),tlcl,sntlcl,buoy_lcl
!# ,buoy(N-1),delta_ht_dry,HT(N),pbesnd(N-1)+pbe_dry
!
! MOIST ADIABATIC PART
!
! WRITE(6,*)' MOIST ADIABATIC PART'
delta_ht_wet=deltah-delta_ht_dry
partem(n) = tsa_fast(blthte,p(n))
wsat=1000.*f_mrsat( p(n)*100., partem(n)+273.15 )
partv=tctotv(partem(n),wsat)
buoy(n)=(partv-tv(n))/tv(n)
pbe_wet = g*0.5*(buoy(n)+buoy_lcl)*delta_ht_wet
pbesnd(n)=pbesnd(n-1) + pbe_dry + pbe_wet
!
wload(n)=0.001*(w(k0)-wsat)
mbuoy(n)=buoy(n) - wload(n)
pbe_wet = g*0.5*(mbuoy(n)+buoy_lcl)*delta_ht_wet
mbesnd(n)=mbesnd(n-1) + pbe_dry + pbe_wet
IF((p(k0)-plcl) < 300.) tcap=AMAX1(tcap,(tv(n)-partv))
END IF ! Upper level below LCL (Dry or bracket)
ELSE ! Lower Level is above LCL
! WRITE(6,*)' GETTING PARCEL TEMPERATURE FOR MOIST CASE'
partem(n) = tsa_fast(blthte,p(n))
wsat=1000.*f_mrsat( p(n)*100., partem(n)+273.15 )
partv=tctotv(partem(n),wsat)
buoy(n)=(partv-tv(n))/tv(n)
pbesnd(n)=pbesnd(n-1)+g*0.5*(buoy(n)+buoy(n-1))*deltah
!
wload(n)=0.001*(w(k0)-wsat)
mbuoy(n)=buoy(n) - wload(n)
mbesnd(n)=mbesnd(n-1)+g*0.5*(mbuoy(n)+mbuoy(n-1))*deltah
IF((p(k0)-p(n)) < 300.) tcap=AMAX1(tcap,(tv(n)-partv))
END IF
! WRITE(6,777)N,P(N),partem(N),T(N),(buoy(n)*1000.),pbesnd(n)
!777 format(' PBE: P,partem,t,b,pbe=',I3,F6.1,4F8.2)
END DO
!
! DETERMINE ENERGY EXTREMA
! Find heights with nuetral buoyancy
!
pos_area=0.
nbe_min=0.
DO n=k0+1,nlevel
! WRITE(6,940)N
!940 format(
! :' LOOKING FOR NEUTRAL BUOYANCY - ENERGY EXTREMUM, LEVEL',I3)
IF((buoy(n)*buoy(n-1)) < 0.)THEN
wlow=buoy(n)/(buoy(n)-buoy(n-1))
htzero=ht(n)*(1.-wlow) + wlow*ht(n-1)
deltah=htzero-ht(n-1)
adjeng=pbesnd(n-1)+g*0.5*buoy(n-1)*deltah
!
IF (p(n) >= 500.) THEN
nbe_min=AMIN1(adjeng,nbe_min)
END IF
!
pos_area=adjeng-nbe_min
pos_max=AMAX1(pos_area,pos_max)
END IF
END DO
! WRITE(6,464)ICP,ICT,N1,NLEVEL
!464 format(' ICP,ICT,N1,NLEVEL',4I5)
!
! Case when equlibrium level is above top of domain
!
pos_area=pbesnd(nlevel)-nbe_min
pos_max=AMAX1(pos_area,pos_max)
!
! At least one region of positive area in sounding
! Make sure there is at least 1 J/kg to avoid some
! round-off errors esp near LCL.
!
IF(pos_max > 1.0)THEN
pbeneg=AMAX1(nbe_min,-400.)
velneg=SQRT(2.0*ABS(pbeneg))
velneg=AMIN1(velneg,20.)
ELSE ! Case when no positive area exists anywhere in sounding
pos_max=0.0
pbeneg =-400.
velneg = 20.
END IF
! WRITE(6,485)pos_max,PBENEG,VELNEG
!485 format(' pos_max',F10.1,' PBENEG',F10.1,' VELNEG',F10.1)
!
! DETERMINE ENERGY EXTREMA FOR MOIST BUOYANCY
! Find heights with nuetral buoyancy
!
pos_area=0.
nbe_min=0.
DO n=k0+1,nlevel
! WRITE(6,940)N
IF((mbuoy(n)*mbuoy(n-1)) < 0.)THEN
wlow=mbuoy(n)/(mbuoy(n)-mbuoy(n-1))
htzero=ht(n)*(1.-wlow) + wlow*ht(n-1)
deltah=htzero-ht(n-1)
adjeng=mbesnd(n-1)+g*0.5*mbuoy(n-1)*deltah
!
IF (p(n) >= 500.) THEN
nbe_min=AMIN1(adjeng,nbe_min)
END IF
!
pos_area=adjeng-nbe_min
mpos_max=AMAX1(pos_area,mpos_max)
END IF
END DO
! WRITE(6,464)ICP,ICT,N1,NLEVEL
!
! Case when equlibrium level is above top of domain
!
pos_area=mbesnd(nlevel)-nbe_min
mpos_max=AMAX1(pos_area,mpos_max)
!
! At least one region of positive area in sounding
! Make sure there is at least 1 J/kg to
! spurious pos energy due to round off.
!
IF(mpos_max > 1.0)THEN
mbeneg=AMAX1(nbe_min,-400.)
mvelneg=SQRT(2.0*ABS(pbeneg))
mvelneg=AMIN1(mvelneg,20.)
ELSE ! Case when no positive area exists anywhere in sounding
mpos_max=0.0
mbeneg =-400.
mvelneg = 20.
END IF
! WRITE(6,486)mpos_max,PBENEG,VELNEG
!486 format(' Mpos_max',F10.1,' MBENEG',F10.1,' mVELNEG',F10.1)
!
! Case when equlibrium level is above top of domain
!
mpos_max = MAX(mpos_max,(mbesnd(nlevel) - nbe_min))
!
! Find EL and LFC
! Unxts are set to km ASL
!
IF(pos_max > 1.0) THEN
IF(buoy(nlevel) > 0.) THEN
nel=nlevel
el=0.001*ht(nlevel)
ELSE
DO n=nlevel-1,k0+1,-1
IF(buoy(n) > 0.) EXIT
END DO
! 1201 CONTINUE
nel=n
wlow=buoy(n+1)/(buoy(n+1)-buoy(n))
el=0.001 * (ht(n+1)*(1.-wlow) + ht(n)*wlow)
END IF
!
DO n=nel,k0,-1
IF(buoy(n) < 0.) EXIT
END DO
! 1301 CONTINUE
IF(n > 0) THEN
wlow=buoy(n+1)/(buoy(n+1)-buoy(n))
lfc=ht(n+1)*(1.-wlow) + ht(n)*wlow
ELSE
lfc=ht(1)
END IF
ELSE
el=0.
lfc=10000.
END IF
lfc=AMIN1(lfc,10000.)
RETURN
END SUBROUTINE potbe1
SUBROUTINE sindex1(maxlev,nlevel,p,ht,t,tv,td,w, & 1,8
partem,buoy,wload,mbuoy,pbesnd,mbesnd, &
lcl_pbe,lfc,el,twdf,li,cape,mcape,cin,mcin,tcap)
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate Convective Available Potential Energy (CAPE)
! in each column of ARPS grid.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Keith Brewster
! April, 1995 Based on OLAPS, hence LAPS, version of same.
!
! MODIFICATION HISTORY:
! 3/11/1996 Cleaned-up, removed OLAPS artifacts.
! 5/13/1996 Added cap strength.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: maxlev,nlevel,k,k0
!
REAL :: p(maxlev),ht(maxlev),t(maxlev),tv(maxlev), &
td(maxlev),w(maxlev)
REAL :: partem(maxlev),buoy(maxlev),wload(maxlev),mbuoy(maxlev)
REAL :: pbesnd(maxlev),mbesnd(maxlev)
!
! Returned from sindx
!
REAL :: lfc,el,twdf,li,cape,mcape,cin,mcin,tcap
!
! Potbe variables
!
REAL :: plcl_pbe,tlcl_pbe,lcl_pbe,thepcl,thepcl0
REAL :: velneg,mvelneg
REAL :: mplcl_pbe,mtlcl_pbe,mlcl_pbe
REAL :: tmp1(maxlev),tmp2(maxlev),tmp3(maxlev),tmp4(maxlev),tmp5(maxlev),tmp6(maxlev)
REAL :: tmp7,tmp8,tmp9,tmp10,tmp11,tmp12,tmp13
!
! Functions
!
REAL :: oe
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
thepcl0=oe(t(1),td(1),p(1))
do k=2,maxlev
if(p(1)-p(k) > 300.0) then
k0=k-1
exit
end if
thepcl=oe(t(k),td(k),p(k))
if(thepcl > thepcl0) thepcl0=thepcl
enddo
!
thepcl=oe(t(1),td(1),p(1))
! print *, ' theta-e of parcel: ',thepcl
!
CALL ptlcl(p(1),t(1),td(1),plcl_pbe,tlcl_pbe)
!
! Find height of LCL
!
CALL intrpr(maxlev,nlevel,p,ht,plcl_pbe,lcl_pbe)
!
! Calculate the CAPE and such
!
CALL potbe1(maxlev,nlevel,p(1),t(1),w(1), 1, &
thepcl,plcl_pbe,tlcl_pbe,lcl_pbe, &
p,ht,t,tv,td,w, &
partem,buoy,wload,mbuoy,pbesnd,mbesnd, &
cin,velneg,cape,mcin,mvelneg,mcape,lfc,el,tcap)
! Calculate Lifted Index
!
CALL calcli(maxlev,nlevel,thepcl,p,t,li)
!
! Calculate max and min wet bulb potential temperature
!
CALL thwxn(maxlev,nlevel,p,ht,t,td,ht(1),twdf)
! Calculate MUCAPE and MUCINS (using k0 level!!!!)
CALL ptlcl(p(k0),t(k0),td(k0),mplcl_pbe,mtlcl_pbe)
CALL intrpr(maxlev,nlevel,p,ht,mplcl_pbe,mlcl_pbe)
CALL potbe1(maxlev,nlevel,p(k0),t(k0),w(k0), k0, &
thepcl0,mplcl_pbe,mtlcl_pbe,mlcl_pbe, &
p,ht,t,tv,td,w, &
tmp1,tmp2,tmp3,tmp4,tmp5,tmp6, &
tmp7,tmp8,tmp9,mcin,tmp10,mcape,tmp11,tmp12,tmp13)
if(mcape < cape) mcape=cape
if(mcin < cin ) mcin =cin
!
RETURN
END SUBROUTINE sindex1
SUBROUTINE enswrtbin (nx,ny, filename, lens, & 1,30
ctime,year,month,day,hour,minute, &
iprgem,nprgem,ihtgem,nhtgem, &
hgtp,uwndp,vwndp,wwndp,tmpp,sphp, &
uwndh,vwndh,wwndh, &
psf,mslp,vort500,temp2m,dewp2m,qv2m, &
u10m,v10m,hgtsfc,refl1km,refl4km,cmpref, &
accppt,convppt,cape,mcape,cin,mcin,lcl, &
srh01,srh03,uh25,sh01,sh06,thck,li,brn,pw)
INCLUDE 'mp.inc'
INTEGER :: nx,ny
CHARACTER (LEN=*) :: filename
CHARACTER (LEN=256) :: filnam
CHARACTER (LEN=40) :: q3dname,q3dunit
INTEGER :: nprgem,nhtgem
INTEGER :: iprgem(nprgem),ihtgem(nhtgem)
REAL :: hgtp(nx,ny,nprgem),uwndp(nx,ny,nprgem),vwndp(nx,ny,nprgem)
REAL :: wwndp(nx,ny,nprgem),tmpp(nx,ny,nprgem),sphp(nx,ny,nprgem)
REAL :: uwndh(nx,ny,nhtgem),vwndh(nx,ny,nhtgem),wwndh(nx,ny,nhtgem)
REAL :: psf(nx,ny),mslp(nx,ny),vort500(nx,ny)
REAL :: temp2m(nx,ny),dewp2m(nx,ny),qv2m(nx,ny)
REAL :: u10m(nx,ny),v10m(nx,ny),hgtsfc(nx,ny)
REAL :: refl1km(nx,ny),refl4km(nx,ny),cmpref(nx,ny)
REAL :: accppt(nx,ny),convppt(nx,ny)
REAL :: cape(nx,ny),mcape(nx,ny),cin(nx,ny),mcin(nx,ny),lcl(nx,ny)
REAL :: srh01(nx,ny),srh03(nx,ny),uh25(nx,ny),sh01(nx,ny),sh06(nx,ny)
REAL :: thck(nx,ny),li(nx,ny),brn(nx,ny),pw(nx,ny)
INTEGER :: i,j,k,klev,iret,igdfil
INTEGER :: year,month,day,hour,minute
INTEGER :: iyr,ifhr,ifmin,ifile,ihd
REAL :: ctime
CHARACTER (LEN=6) :: hgtptag(5),uwndptag(5)
CHARACTER (LEN=6) :: vwndptag(5),wwndptag(5)
CHARACTER (LEN=6) :: tmpptag(5),qvptag(5)
CHARACTER (LEN=6) :: uwndhtag(6),vwndhtag(6),wwndhtag(6)
!
DATA hgtptag / 'hgt850','hgt700','hgt600','hgt500','hgt250'/
DATA uwndptag / 'u850__','u700__','u600__','u500__','u250__'/
DATA vwndptag / 'v850__','v700__','v600__','v500__','v250__'/
DATA wwndptag / 'w850__','w700__','w600__','w500__','w250__'/
DATA tmpptag / 'tmp850','tmp700','tmp600','tmp500','tmp250'/
DATA qvptag / 'sph850','sph700','sph600','sph500','sph250'/
!
DATA uwndhtag /'u1km__','u2km__','u3km__','u4km__','u5km__','u6km__'/
DATA vwndhtag /'v1km__','v2km__','v3km__','v4km__','v5km__','v6km__'/
DATA wwndhtag /'w1km__','w2km__','w3km__','w4km__','w5km__','w6km__'/
DO i=0,nprocs-1,max_fopen
IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
!-----------------------------------------------------------------------
!
! Output one level fields
!
!-----------------------------------------------------------------------
!
IF(myproc==0) PRINT *, ' Writing surface pressure'
filnam = filename(1:lens-9)//'sfpres'// &
filename(lens-5:lens)
q3dname='SFC PRES'
q3dunit='hPa'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,psf)
IF(myproc==0) PRINT *, ' Writing MSLP pressure'
filnam = filename(1:lens-9)//'mspres'// &
filename(lens-5:lens)
q3dname='MSLP'
q3dunit='hPa'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,mslp)
IF(myproc==0) PRINT *, ' Writing 1-h total accumulated rainfall'
filnam = filename(1:lens-9)//'accppt'// &
filename(lens-5:lens)
q3dname='RAIN'
q3dunit='mm'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,accppt)
IF(myproc==0) PRINT *, ' Writing precipitable water'
filnam = filename(1:lens-9)//'pwat__'// &
filename(lens-5:lens)
q3dname='PWAT'
q3dunit='mm'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,pw)
IF(myproc==0) PRINT *, ' Writing 2-m temperature (F)'
filnam = filename(1:lens-9)//'temp2m'// &
filename(lens-5:lens)
q3dname='temp2m'
q3dunit='F'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,temp2m)
IF(myproc==0) PRINT *, ' Writing 2-m dew point temperature (F)'
filnam = filename(1:lens-9)//'dewp2m'// &
filename(lens-5:lens)
q3dname='dewp2m'
q3dunit='F'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,dewp2m)
IF(myproc==0) PRINT *, ' Writing 2-m specific humidity'
filnam = filename(1:lens-9)//'qv2m__'// &
filename(lens-5:lens)
q3dname='qv2m'
q3dunit='g/kg'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,qv2m)
IF(myproc==0) PRINT *, ' Writing 10-m u velocity'
filnam = filename(1:lens-9)//'u10m__'// &
filename(lens-5:lens)
q3dname='u10m'
q3dunit='m/s'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,u10m)
IF(myproc==0) PRINT *, ' Writing 10-m v velocity'
filnam = filename(1:lens-9)//'v10m__'// &
filename(lens-5:lens)
q3dname='v10m'
q3dunit='m/s'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,v10m)
IF(myproc==0) PRINT *, ' Writing surface geopotential height'
filnam = filename(1:lens-9)//'hgtsfc'// &
filename(lens-5:lens)
q3dname='terrain'
q3dunit='m'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,hgtsfc)
! IF(myproc==0) PRINT *, ' Writing 500 hPa vorticity'
! filnam = filename(1:lens-9)//'vrt500'// &
! filename(lens-5:lens)
! q3dname='vort500'
! q3dunit='m/s'
! CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,vort500)
IF(myproc==0) PRINT *, ' Writing 1-km AGL reflectivity'
filnam = filename(1:lens-9)//'ref1km'// &
filename(lens-5:lens)
q3dname='refl1km'
q3dunit='dBZ'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,refl1km)
IF(myproc==0) PRINT *, ' Writing 4-km AGL reflectivity'
filnam = filename(1:lens-9)//'ref4km'// &
filename(lens-5:lens)
q3dname='refl4km'
q3dunit='dBZ'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,refl4km)
IF(myproc==0) PRINT *, ' Writing composite reflectivity'
filnam = filename(1:lens-9)//'cmpref'// &
filename(lens-5:lens)
q3dname='cmpref'
q3dunit='dBZ'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,cmpref)
IF(myproc==0) PRINT *, ' Writing surface-based CAPE'
filnam = filename(1:lens-9)//'sbcape'// &
filename(lens-5:lens)
q3dname='SBCAPE'
q3dunit='J/kg'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,cape)
IF(myproc==0) PRINT *, ' Writing moist unstable CAPE'
filnam = filename(1:lens-9)//'mucape'// &
filename(lens-5:lens)
q3dname='MUCAPE'
q3dunit='J/kg'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,mcape)
IF(myproc==0) PRINT *, ' Writing surface-based CIN'
filnam = filename(1:lens-9)//'sbcins'// &
filename(lens-5:lens)
q3dname='SBCINS'
q3dunit='J/kg'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,cin)
IF(myproc==0) PRINT *, ' Writing moist unstable CIN'
filnam = filename(1:lens-9)//'mucins'// &
filename(lens-5:lens)
q3dname='MUCINS'
q3dunit='J/kg'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,mcin)
IF(myproc==0) PRINT *, ' Writing surface-based LCL'
filnam = filename(1:lens-9)//'sblcl_'// &
filename(lens-5:lens)
q3dname='SBLCL'
q3dunit='m'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,lcl)
IF(myproc==0) PRINT *, ' Writing 0-1 km AGL storm-relative helicity'
filnam = filename(1:lens-9)//'srh01_'// &
filename(lens-5:lens)
q3dname='SRH01'
q3dunit='m^2/s^2'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,srh01)
IF(myproc==0) PRINT *, ' Writing 0-3 km AGL storm-relative helicity'
filnam = filename(1:lens-9)//'srh03_'// &
filename(lens-5:lens)
q3dname='SRH03'
q3dunit='m^2/s^2'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,srh03)
IF(myproc==0) PRINT *, ' Writing updarft helicity'
filnam = filename(1:lens-9)//'vhel__'// &
filename(lens-5:lens)
q3dname='VHEL'
q3dunit='m^2/s^2'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,uh25)
IF(myproc==0) PRINT *, ' Writing 0-1 km AGL wind shear'
filnam = filename(1:lens-9)//'shr01_'// &
filename(lens-5:lens)
q3dname='SHR01'
q3dunit='1/s'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,sh01)
IF(myproc==0) PRINT *, ' Writing 0-6 km AGL wind shear'
filnam = filename(1:lens-9)//'shr06_'// &
filename(lens-5:lens)
q3dname='SHR06'
q3dunit='1/s'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,sh06)
!
!-----------------------------------------------------------------------
!
! Output constant pressure level data
!
!-----------------------------------------------------------------------
IF( nprgem > 0 ) THEN
DO klev=1,nprgem
IF(myproc==0) &
PRINT *, ' Writing binary data at pr= ',iprgem(klev),' hPa'
filnam = filename(1:lens-9)//hgtptag(klev)// &
filename(lens-5:lens)
q3dname='HGHT'
q3dunit='m'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,hgtp(1,1,klev))
filnam = filename(1:lens-9)//uwndptag(klev)// &
filename(lens-5:lens)
q3dname='UREL'
q3dunit='m/s'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,uwndp(1,1,klev))
filnam = filename(1:lens-9)//vwndptag(klev)// &
filename(lens-5:lens)
q3dname='VREL'
q3dunit='m/s'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,vwndp(1,1,klev))
filnam = filename(1:lens-9)//wwndptag(klev)// &
filename(lens-5:lens)
q3dname='WREL'
q3dunit='m/s'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,wwndp(1,1,klev))
filnam = filename(1:lens-9)//tmpptag(klev)// &
filename(lens-5:lens)
q3dname='TMPC'
q3dunit='C'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,tmpp(1,1,klev))
filnam = filename(1:lens-9)//qvptag(klev)// &
filename(lens-5:lens)
q3dname='QV'
q3dunit='g/kg'
CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,sphp(1,1,klev))
END DO ! nprgem-loop
END IF
!
!-----------------------------------------------------------------------
!
! Output AGL height level data
!
!-----------------------------------------------------------------------
!
! IF( nhtgem > 0 ) THEN
! DO klev=1,nhtgem
! IF(myproc == 0) &
! PRINT *, ' Writing binary data at z= ',ihtgem(klev),' km AGL'
! filnam = filename(1:lens-9)//uwndhtag(klev)// &
! filename(lens-5:lens)
! q3dname='UREL'
! q3dunit='m/s'
! CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,uwndh(1,1,klev))
! filnam = filename(1:lens-9)//vwndhtag(klev)// &
! filename(lens-5:lens)
! q3dname='VREL'
! q3dunit='m/s'
! CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,vwndh(1,1,klev))
! filnam = filename(1:lens-9)//wwndhtag(klev)// &
! filename(lens-5:lens)
! q3dname='WREL'
! q3dunit='m/s'
! CALL bindump2d(nx,ny,trim(filnam),q3dname,q3dunit,wwndh(1,1,klev))
! END DO ! nhtgem-loop
! END IF
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
RETURN
950 CONTINUE
IF(myproc == 0) WRITE(6,'(a)') ' Error setting up binary file'
STOP
RETURN
END SUBROUTINE enswrtbin