! !################################################################## !################################################################## !###### ###### !###### 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