!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE QPFGRID ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE qpfgrid(nx,ny,nz,prcrate,pprt,ptprt,qv,pbar,ptbar, & 1
rhostr,zp,j3,j3inv,raing,temxy1,temxy2,pi,dqv)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the source/sink terms in temperature and moisture
! equations as a results of large-scale (as opposed to subgird scale)
! condensation. Apply the source terms to ptprt and qv.
!
! The surface gridscale rainfall is accumulated in raing.
!
!-----------------------------------------------------------------------
!
! AUTHOR: V. Wong, L. Zhao and X. Song
! 03/12/95
!
! MODIFICATION HISTORY:
!
! 05/10/95 (L. Zhao)
! Correction made to the calculation of accumulated rainfall
! for restart runs.
!
! 8/12/95 (M. Xue)
! Rearranged the argument list to follow ARPS convention.
!
! 01/31/1996 (V. Wong and X. Song)
! Related the accumulated grid-scale precipitation to the parameter,
! qpfgfrq, that controls the the frequency of calling QPFGRID.
!
!-----------------------------------------------------------------------
!
! 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
!
! pprt Perturbation pressure (Pascal)
! ptprt Perturbation potential temperature at a given time
! level (K)
! qv Water vapor specific humidity at a given time level
! (kg/kg)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
!
! OUTPUT:
!
! raing Accumulated gridscale rainfall (mm).
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! The number grid points in 3 directions
REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal)
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: qv (nx,ny,nz) ! Water vapor specific humidity (kg/kg)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: rhostr (nx,ny,nz) ! Base state density rhobar times j3.
REAL :: zp (nx,ny,nz) ! The height of the terrain.
REAL :: j3 (nx,ny,nz) ! Coordinate transformation
! Jacobian d(zp)/d(z)
REAL :: j3inv (nx,ny,nz) ! Coordinate transformation
! Jacobian d(zp)/d(z)
REAL :: raing (nx,ny) ! Accumulated grid precipitation (mm)
REAL :: prcrate(nx,ny) ! precipitation rate (kg/(m**2*s))
REAL :: temxy1 (nx,ny) ! 2-D work array
REAL :: temxy2 (nx,ny) ! 2-D work array
REAL :: pi (nx,ny) ! 2-D work array
REAL :: dqv (nx,ny) ! 2-D work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
REAL :: rhmin,rhmax,al2ovcr,alovcp
DATA rhmin/0.05/
DATA rhmax/1.0/
DATA al2ovcr/1.346946E+07/ !Lv**2/(CP*Rv)
DATA alovcp/2.487458E+03/ !Lv/Cp
REAL :: es,temp2df,pres,tdf,theta
REAL :: p0inv,pterm,tema,temb
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
p0inv=1./p0
!
!-----------------------------------------------------------------------
!
! Initializing the working array for each (i,j) point at beginning.
!
!-----------------------------------------------------------------------
!
DO i=1,nx-1
DO j=1,ny-1
temxy1(i,j)=0.0
END DO
END DO
!
!-----------------------------------------------------------------------
!
DO k=nz-2,2,-1
!
!-----------------------------------------------------------------------
!
DO i=1,nx-1
DO j=1,ny-1
!
!-----------------------------------------------------------------------
!
! Computer needed variables from model inputs. p0 and rddcp are
! included in 'physct.inc'
!
!-----------------------------------------------------------------------
!
pres=pprt(i,j,k) + pbar(i,j,k) ! p
pi(i,j)=(pres*p0inv)**rddcp ! pi
tdf=(ptbar(i,j,k)+ptprt(i,j,k))*pi(i,j) ! temperature
!
!-----------------------------------------------------------------------
!
! Calculate saturation vapor pressure(ES) at (P,T).
!
!-----------------------------------------------------------------------
!
pterm = 0.5 + SIGN(0.5, tdf-260.0)
tema = pterm*17.2694 + (1.-pterm)*21.87456
temb = pterm*35.86 + (1.-pterm)*7.66
es = 610.78*EXP(tema*(tdf-273.16)/(tdf-temb))
!
!-----------------------------------------------------------------------
!
! Computer the saturation specific humidity.
!
!-----------------------------------------------------------------------
!
dqv(i,j)=0.622*es/(pres-0.378*es)
!
!-----------------------------------------------------------------------
!
! Set a lower limit to 'qv field' to avoid negative values.
!
!-----------------------------------------------------------------------
!
qv(i,j,k) = MAX( qv(i,j,k), rhmin*dqv(i,j) )
!
!-----------------------------------------------------------------------
!
! Set the upper limit of 'qv field' such that the RHMAX is equal to
! or less than 1.0.
!
! Modified by Zuwen He (04/2002)
! adopt new naming convention: rhsat
!
! original code: dqv(i,j)=rhmax*dqv(i,j)
! new code: dqv(i,j)=rhsat*dqv(i,j)
!
!-----------------------------------------------------------------------
!
dqv(i,j)=rhsat*dqv(i,j)
!
!-----------------------------------------------------------------------
!
! Compute Dq/Dt term.
!
!-----------------------------------------------------------------------
!
temp2df=1.0+al2ovcr*dqv(i,j)/(tdf*tdf)
!
!-----------------------------------------------------------------------
!
! Check for condensation.
!
!-----------------------------------------------------------------------
!
dqv(i,j)=(qv(i,j,k)-dqv(i,j))/temp2df
IF (dqv(i,j) >= 0.0) temxy1(i,j)=temxy1(i,j)+dqv(i,j) &
*rhostr(i,j,k)*(zp(i,j,k+1)-zp(i,j,k))*j3inv(i,j,k)
!
!-----------------------------------------------------------------------
!
! Check for evaporation.
!
!-----------------------------------------------------------------------
!
temp2df=temxy1(i,j)+dqv(i,j) &
*rhostr(i,j,k)*(zp(i,j,k+1)-zp(i,j,k))*j3inv(i,j,k)
IF(temp2df < 0) THEN
temp2df=0.0
dqv(i,j)=-temxy1(i,j) &
/(rhostr(i,j,k)*(zp(i,j,k+1)-zp(i,j,k))*j3inv(i,j,k))
END IF
!
IF(dqv(i,j) < 0)temxy1(i,j)=temp2df
temxy2(i,j)=tdf
temxy1(i,j)=MAX(0.0,temxy1(i,j)) !keep it positive
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Modify T and q fields at the Kth level.
!
!-----------------------------------------------------------------------
!
IF (curtim > dtbig) THEN
DO i=1,nx-1
DO j=1,ny-1
qv(i,j,k)=MAX(0.0,qv(i,j,k)-dqv(i,j))
tdf=temxy2(i,j) + alovcp*dqv(i,j)
theta=tdf/pi(i,j)
ptprt(i,j,k)=theta-ptbar(i,j,k)
END DO
END DO
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Add in condensation to stable preiciptation array (accumulated
! grid precipitation). The value is cumulated until the cumulus
! parameterization is turn on in the mode, which is control by
! 'confrg'. The unit for rainfall is 'mm'.
!
!-----------------------------------------------------------------------
!
! IF (mod(curtim+0.001,confrq).lt.(0.5*dtbig)) THEN
DO i=1,nx-1
DO j=1,ny-1
tema=temxy1(i,j)/rhow
! raing(i,j)=raing(i,j) +
! : 1000.0* temxy1(i,j) * ( pprt(i,j,2)+pbar(i,j,2) -
! : pprt(i,j,3)-pbar(i,j,3) ) /(rhow*g)
raing(i,j)=raing(i,j) + 1000.0*tema
prcrate(i,j)=tema*rhow/qpfgfrq
END DO
END DO
! ENDIF
RETURN
END SUBROUTINE qpfgrid
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE QPFCUMS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE qpfcums(nx,ny,nz,prcrate,qvsflx, & 1,8
u,v,w,pprt,ptprt,qv,pbar,ptbar,qvbar,rhostr,zp,j3, &
ptcumsrc,qcumsrc,rainc,nca,kfraincv, &
cldefi,xland,bmjraincv, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9,tem10)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the source/sink terms in temperature and moisture
! equations as a results of subgrid scale (as opposed to grid scale)
! cumulus convection.
!
! Kuo, Kain-Fritsch, and WRF BMJ schemes are used here.
!
! The surface convective rainfall is accumulated in rainc.
!
!-----------------------------------------------------------------------
!
! AUTHOR: V. Wong, L. Zhao and X. Song
! 3/06/95
!
! MODIFICATION HISTORY:
!
! 8/12/95 (M. Xue)
! Rearranged the argument list to follow ARPS conventions.
!
! 08/01/97 (Zonghui Huo)
! Added Kain-fritsch cumulus parameterization scheme.
!
! 13 March 2002 (Eric Kemp)
! Added WRF BMJ cumulus parameterization scheme.
!
! April 2002 (Fanyou Kong)
! Added WRF new Kain-Fritsch (April 2002 version: KF_ETA) scheme
!
!-----------------------------------------------------------------------
!
! 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
!
! u x component of velocity at a given time level (m/s)
! v y component of velocity at a given time level (m/s)
! w * z component of velocity at a given time level (m/s)
! for KUO scheme
! * time-average vertical velocity (m/s)
! for Kain Fritsch scheme
!
! pprt Perturbation pressure (Pascal)
! ptprt Perturbation potential temperature (K)
! qv Water vapor specific humidity (kg/kg)
! pbar Base state pressure (Pascal)
! ptbar Base state potential temperature (K)
! qvbar Base state water vapor specific humidity (kg/kg)
! rhostr Base state density rhobar times j3 (kg/m**3)
!
! zp The physical height coordinate defined at w-point
! j3 Coordinate transformation Jacobian, d(zp)/d(z)
!
! OUTPUT:
!
! ptcumsrc Potential temperature source term.
! qcumsrc Water specific humidity source term
! rainc Accumulated precipitation by cumulus convection (mm)
! prcrate precipitatioon rate (mm/s)
! kfraincv K-F convective rainfall (cm)
! nca K-F counter for CAPE release
! cldefi BMJ cloud efficiency
! xland BMJ land/sea mask
! bmjraincv BMJ convective rainfall (cm)
!
! WORK ARRAYS:
!
! tem1 Temporary work array
! tem2 Temporary work array
! tem3 Temporary work array
! tem4 Temporary work array
! tem5 Temporary work array
! tem6 Temporary work array
! tem7 Temporary work array
! tem8 Temporary work array
! tem9 Temporary work array
! tem10 Temporary work array
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! The number grid points in 3 directions
REAL :: u (nx,ny,nz) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz) ! Total v-velocity (m/s)
REAL :: w (nx,ny,nz) ! Total w-velocity (m/s)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal)
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: qv (nx,ny,nz) ! Water vapor specific humidity (kg/kg)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific
REAL :: rhostr (nx,ny,nz) ! Base state density rhobar times j3.
REAL :: zp (nx,ny,nz) ! The height of the terrain.
REAL :: j3 (nx,ny,nz) ! Coordinate transformation
! Jacobian d(zp)/d(z)
REAL :: ptcumsrc(nx,ny,nz) ! Source term in pt equation.
REAL :: qcumsrc(nx,ny,nz,5) ! Source term in water equations due
! to cumulus parameterization:
! qcumsrc(1,1,1,1) for qv equation
! qcumsrc(1,1,1,2) for qc equation
! qcumsrc(1,1,1,3) for qr equation
! qcumsrc(1,1,1,4) for qi equation
! qcumsrc(1,1,1,5) for qs equation
REAL :: rainc(nx,ny) ! Accumulated precipitation by convection
REAL :: prcrate(nx,ny) ! precipitation rate (kg/(m**2*s))
REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s))
REAL :: kfraincv (nx,ny) ! K-F convective rainfall (cm)
INTEGER :: nca (nx,ny) ! K-F counter for CAPE release
!EMK BMJ
REAL,INTENT(INOUT) :: cldefi(nx,ny) ! BMJ cloud efficiency
REAL,INTENT(IN) :: xland(nx,ny) ! BMJ land mask
! (1.0 = land, 2.0 = sea)
REAL,INTENT(INOUT) :: bmjraincv(nx,ny) ! BMJ convective rainfall (cm)
!EMK END
REAL :: tem1 (nx,ny,nz) ! Temporary work array
REAL :: tem2 (nx,ny,nz) ! Temporary work array
REAL :: tem3 (nx,ny,nz) ! Temporary work array
REAL :: tem4 (nx,ny,nz) ! Temporary work array
REAL :: tem5 (nx,ny,nz) ! Temporary work array
REAL :: tem6 (nx,ny,nz) ! Temporary work array
REAL :: tem7 (nx,ny,nz) ! Temporary work array
REAL :: tem8 (nx,ny,nz) ! Temporary work array
REAL :: tem9 (nx,ny,nz) ! Temporary work array
REAL :: tem10(nx,ny,nz) ! Temporary work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
REAL :: p0inv,dthmax
REAL :: kqmax,jqmax,iqmax
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'cumucst.inc'
INCLUDE 'mp.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
p0inv=1./p0
IF(cnvctopt == 3) THEN
CALL set_acct
(kfcum_acct)
IF (myproc == 0) &
WRITE(6,'(1x,a)') 'Use Kain-Fritcsh Scheme,Enter kfinterfc'
CALL kfinterfc(nx,ny,nz,u,v,w,pprt,ptprt,qv,pbar,ptbar,zp, &
ptcumsrc,qcumsrc,rainc,prcrate, &
nca,kfraincv, &
tem1(1,1,1),tem1(1,1,2),tem1(1,1,3), &
tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9,tem10)
!KFY KF_ETA
ELSE IF (cnvctopt == 5) THEN
CALL set_acct(kfcum_acct)
IF (myproc == 0) &
WRITE(6,'(1x,a)') 'Use WRF Kain-Fritcsh Scheme'
CALL interface_wrf_kfetadrv(nx,ny,nz,u,v,w, &
pprt,ptprt,qv,pbar,ptbar,zp, &
ptcumsrc,qcumsrc,prcrate, &
nca,kfraincv)
!KFY END
!EMK BMJ
ELSE IF (cnvctopt == 4) THEN
CALL set_acct(bmjcum_acct)
IF (myproc == 0) &
WRITE(6,'(1x,a)') 'Use WRF Betts-Miller-Janjic Scheme'
CALL interface_wrf_bmjdrv(nx,ny,nz,pprt,ptprt,qv,pbar,ptbar,zp, &
ptcumsrc,qcumsrc,bmjraincv,prcrate,cldefi,&
xland)
!EMK END
ELSE ! use Kuo-scheme
CALL set_acct(kuocum_acct)
IF (myproc == 0) &
WRITE(6,'(1x,a)') 'Use KUO cumulus scheme'
!
!-----------------------------------------------------------------------
!
! For each grid point, add the convective parameterization.
!
!-----------------------------------------------------------------------
!
dthmax=0.0
DO j=1,ny-1
DO i=1,nx-1
msflux=-qvsflx(i,j)
DO k=1,nz-1
ucon(k)=(u(i,j,k)+u(i+1,j,k))/2.0
vcon(k)=(v(i,j,k)+v(i,j+1,k))/2.0
wcon(k)=w(i,j,k)
picon(k)=cp*((pprt(i,j,k) + pbar(i,j,k))*p0inv)**rddcp
thtcon(k)=ptbar(i,j,k)+ptprt(i,j,k)
tmpcon(k)=thtcon(k)*picon(k)/cp
!
! DNCON = rhobar
!
dncon(k)=rhostr(i,j,k)
qvcon(k)=qv(i,j,k)
ptcumsrc(i,j,k)=0. !initializing forcing terms
qcumsrc(i,j,k,1)=0.
!
zzcon(k) =zp(i,j,k)
zcon(k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
END DO
!
!-----------------------------------------------------------------------
!
! Calculate enviromental parameters.
!
!-----------------------------------------------------------------------
!
CALL environ(nz-1)
!
!-----------------------------------------------------------------------
!
! IF there is convective heating, Kuo's cumulus parameterization
! scheme is used.
!
!-----------------------------------------------------------------------
!
IF(igo /= 0) THEN
CALL kuocp
END IF
!
!-----------------------------------------------------------------------
!
! If there is convective heating, adjust and output the forcnig
! terms into the model.
!
!-----------------------------------------------------------------------
!
IF(igo /= 0) THEN
CALL cp2mod (nz-1)
!
DO k=2,nz-2
ptcumsrc(i,j,k)=ftcon(k)
qcumsrc(i,j,k,1)=frcon(k)
END DO
! IF (mod(curtim+0.001,confrq).lt.(0.5*dtbig))THEN
rainc(i,j)=rainc(i,j)+1000.0*cprecip*confrq/rhow
prcrate(i,j)=cprecip
! ENDIF
!
!-----------------------------------------------------------------------
!
! Find the location of maximum convective heating rate.
!
!-----------------------------------------------------------------------
!
DO k=2,nz-2
IF(ptcumsrc(i,j,k) > dthmax) THEN
dthmax=ptcumsrc(i,j,k)
iqmax=i
jqmax=j
kqmax=k
END IF
END DO
END IF
END DO
END DO
END IF ! End of cumulus schemes
RETURN
END SUBROUTINE qpfcums
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE KUOCP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE kuocp 1,2
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Performs the Kuo's convective cumulus parameterization.
!
!-----------------------------------------------------------------------
!
! AUTHOR: V. Wong, L. Zhao and X. Song
! 3/06/95
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: k
INTEGER :: kdiv,klfs,kover,kdet
REAL :: supply
! real SUPPLYW
REAL :: zdetr,bkuo,dzdet,wtgnd,wtdiv,preff,envshr,vhint,overmax, &
factr,vdint,vmint,avtdiff,avgmin,wtlcl, &
anegl,dzdiv,dddt,wtlfs,anegh, &
apos
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'cumucst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Downdraft flag - 0 - no downdrafts
! 1 - simple downdraft model
!
!-----------------------------------------------------------------------
!
! Initializing the data array for theta and qv forcing terms.
!
!-----------------------------------------------------------------------
!
DO k=1,nkp
ftcon(k)=0.
frcon(k)=0.
END DO
!
!-----------------------------------------------------------------------
!
! Compute vertical moisture convergence into the cloud layer.
! Neglect the horizontal advection of specific humidity since
! it can be shown to be relatively small. Vertical flux out cloud
! top is also assumed small.
!
!-----------------------------------------------------------------------
!
! SUPPLYW=RHOE(KLCL)*RVE(KLCL)*(WPE(KLCL)+WPE(KLCL-1))*.5
!
! SUPPLY=MSFLUX-0.125*(RHOE(2)+RHOE(1))
! : *(WPE(2)+WPE(1))*(RVE(2)-RVE(1))
supply=msflux-0.125*(rhoe(klcl+1)+rhoe(klcl)) &
*(wpe(klcl+1)+wpe(klcl)) &
*(rve(klcl+1)-rve(klcl))
! DO 105 K=3,KMT-1
DO k=klcl+2,kmt-1
supply=supply-0.25*(rhoe(k)+rhoe(k-1)) &
*(wpe(k)+wpe(k-1))*(rve(k)-rve(k-1))
END DO
supply=supply-0.125*(rhoe(kmt)+rhoe(kmt-1)) &
*(wpe(kmt)+wpe(kmt-1))*(rve(kmt)-rve(kmt-1))
! SUPPLY=SUPPLYW
!
!-----------------------------------------------------------------------
!
! If water supply is not large enough, no convective is allowed.
!
!-----------------------------------------------------------------------
!
IF (supply <= 1.0E-3) THEN
igo=0
RETURN
END IF
!
!-----------------------------------------------------------------------
!
! This is the cloud model. Updraft is constant THETA e and
! saturated with respect to water. There is no ice.
! Cloud top is one level above ETL.
!
! THETA e of the updraft
!
!-----------------------------------------------------------------------
!
theu(klcl)=(3.376/tlcl-0.00254)*1000*rve(kcon) &
*(1+0.81*rve(kcon))
theu(klcl)=the(kcon)*EXP(theu(klcl))
PRINT*,'THEU(KLCL)=',theu(klcl),' KLCL=',klcl
!
!-----------------------------------------------------------------------
!
! Equilibrium Temperature Level of the source level air.
!
!-----------------------------------------------------------------------
!
igo=0 ! index variable for convection
DO k=klcl,kmt
CALL the2t(theu(klcl),pe(k),thu(k),tu(k),rsu(k))
IF(thu(k) > the(k).AND.igo == 0) THEN
igo=1
klfc=k ! level of free convection
END IF
IF(thu(k) <= the(k).AND.igo == 1) GO TO 500 !cloud top
END DO
IF (igo == 0) RETURN
PRINT*,' Convection beyond model top - THup, THenv ',thu(kmt) &
,the(kmt)
k=kmt-1
500 CONTINUE
!
! KETL=MIN(K,KMT)
! KCT=MIN(KETL+1,KMT)
!
ketl=MIN(k-1,kmt)
kct=MIN(ketl,kmt)
CALL the2t(theu(klcl),pe(kct),thu(kct),tu(kct),rsu(kct))
DO k=1,klfc-1
thu(k)=the(k)
END DO
!
!-----------------------------------------------------------------------
!
! If the cloud is not at least CDZMIN deep or cloud top is
! under 500 mb, no convection is allowed.
!
!-----------------------------------------------------------------------
!
IF (ze(ketl)-ze(klfc) < cdzmin.OR.pe(kct) > 50000.) THEN
igo=0
RETURN
END IF
!
!-----------------------------------------------------------------------
!
! Require the positive area be 50% greater than the negative
! area below the LFC and 5% greater in total.
!
!-----------------------------------------------------------------------
!
anegl=0. ! negative area
DO k=klcl,klfc-1
anegl=anegl+(thu(k)-the(k))*(zc(k)-zc(k-1))
END DO
!
apos=0. ! positive area
DO k=klfc,ketl-1
apos=apos+(thu(k)-the(k))*(zc(k)-zc(k-1))
END DO
!
anegh=0.
DO k=ketl,kct
anegh=anegh+(thu(k)-the(k))*(zc(k)-zc(k-1))
END DO
!
IF(apos < ABS(anegl)*1.5.OR.apos < ABS(anegl+anegh)*1.05) THEN
igo=0
RETURN
END IF
!
!-----------------------------------------------------------------------
!
! The downdraft model - starts at THETA e minimum (LFS).
! Downdraft is 2 degrees colder than
! environment at cloud base increasing to 5 degrees
! colder at the ground.
!
! Find LFS as THETA e minimum
!
!-----------------------------------------------------------------------
!
IF (idownd == 1) THEN
!
DO k=kct,2,-1
IF (thee(k) < thee(k+1).AND.thee(k) < thee(k-1)) GO TO 510
END DO
k=2
510 CONTINUE
klfs=k
IF(klfs <= klcl)klfs=klcl+1
thd(klfs)=the(klfs)
!
!-----------------------------------------------------------------------
!
! Limit dd deficit at the ground to the maximum of positive
! temperature difference of updraft if less than 2.5 degrees.
!
!-----------------------------------------------------------------------
!
dddt=0.
DO k=klcl,kct
dddt=MAX(dddt,thu(k)-the(k))
END DO
IF(dddt > 2.5) dddt=5.
!
thd(2)=the(2)-dddt
thd(klcl)=the(klcl)-dddt*.2
DO k=klcl,klfs
thd(k)=thd(klcl)+(thd(klfs)-thd(klcl))/(ze(klfs)-ze(klcl)) &
*(ze(k)-ze(klcl))
END DO
DO k=3,klcl-1
thd(k)=thd(2)+(thd(klcl)-thd(2))/(ze(klcl)-ze(2)) &
*(ze(k)-ze(2))
END DO
!
!-----------------------------------------------------------------------
!
! Now we need to weight the downdraft relative to the updraft.
! Assume that the dd weight is zero at the LFS, 1/2 of
! updraft at cloud base, and equal to the updraft at cloud
! base at the ground.
!
!-----------------------------------------------------------------------
!
dzdiv=1E20
DO k=1,kmt
IF(ABS(ze(k)-800.) < dzdiv)THEN
kdiv=k
dzdiv=ABS(ze(k)-800.)
END IF
END DO
kdiv=MAX(MIN(klcl,kdiv),2)
IF(kdiv == klcl) kdiv=klcl-1
DO k=1,nkp
wtd(k)=0.
END DO
wtlfs=0.
wtlcl=.1
wtdiv=.2
wtgnd=1.
DO k=klcl+1,klfs
wtd(k)=wtlcl+(wtlfs-wtlcl)/(ze(klfs)-ze(klcl)) &
*(ze(k)-ze(klcl))
END DO
DO k=kdiv,klcl
wtd(k)=wtdiv+(wtlcl-wtdiv)/(ze(klcl)-ze(kdiv)) &
*(ze(k)-ze(kdiv))
END DO
DO k=2,kdiv-1
wtd(k)=wtgnd+(wtdiv-wtgnd)/(ze(kdiv)-ze(2)) &
*(ze(k)-ze(2))
END DO
ELSE ! for idownd=0 option
DO k=1,nkp ! set the downdraft weights to zero
wtd(k)=0.
END DO
DO k=2,klcl-1 ! below cloud base, no difference
thu(k)=the(k)
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Compute the b parameter. Use Fritsch/Chappell's precipitation
! efficiency. Notice that the unit of (DV/DZ) is 10e-3 s-1.
!
!-----------------------------------------------------------------------
!
envshr=SQRT((upe(kct)-upe(klfc))**2 &
+(vpe(kct)-vpe(klfc))**2) &
/(ze(kct)-ze(klfc))*1E3
IF (envshr > 1.35) THEN
preff=1.591-.639*envshr+.0953*envshr**2-.00496*envshr**3
ELSE
preff=.9 !parameter (1-b)
END IF
preff=MAX(0.0,preff) !keep it positive
bkuo=1.-preff !parameter b
!
!-----------------------------------------------------------------------
!
! Initializing the vertical profiles of convective heating and
! moistening
!
!-----------------------------------------------------------------------
!
DO k=2,kmt
vheat(k)=0.
vmois(k)=0.
vmdry(k)=0.
END DO
!
!-----------------------------------------------------------------------
!
! Find the weighted THETA to use for the convection. Where WTD is the
! weights for downdraft. When WTD(k)=0, THCON(k)=THU(k).
!
!-----------------------------------------------------------------------
!
DO k=2,kct
thcon(k)=wtd(k)*thd(k)+(1.-wtd(k))*thu(k)
END DO
!
!-----------------------------------------------------------------------
!
! Heating profile is difference between convective THETAs and
! environment.
!
!-----------------------------------------------------------------------
!
DO k=2,kct
vheat(k)=thcon(k)-the(k) !delta theta
END DO
!
!-----------------------------------------------------------------------
!
! Moisture profile is difference between vapor's of updraft and
! environment in the cloud layer. Below cloud base, air is
! dried by SUPPLY. Downdrafts are assumed to have no effect
! on this.
!
!-----------------------------------------------------------------------
!
zdetr=.66667*ze(kct)
dzdet=1000000.
DO k=klcl,kct
IF(ABS(ze(k)-zdetr) < dzdet)THEN
dzdet=ABS(ze(k)-zdetr)
kdet=k
END IF
END DO
DO k=kdet,kct
vmois(k)=1.
END DO
DO k=klcl,kdet-1
vmois(k)=rsu(k)-rve(k)
END DO
DO k=2,klcl-1
vmdry(k)=rve(k)
END DO
vhint=0.
vmint=0.
vdint=0.
DO k=2,kmt
vhint=vhint+vheat(k)*(zc(k)-zc(k-1)) ! sum (dtheta*dz)
vmint=vmint+vmois(k)*(zc(k)-zc(k-1)) ! sum(dq*dz) for moist case
vdint=vdint+vmdry(k)*(zc(k)-zc(k-1)) ! sum(dq*dz) for dry case
END DO
!
!-----------------------------------------------------------------------
!
! If VHINT is less than 0, there is more negative area than
! positive area. No convection allowed.
!
!-----------------------------------------------------------------------
!
IF (vhint <= 0.) THEN
igo=0
RETURN
END IF
!
!-----------------------------------------------------------------------
!
! Also require that there is a minimum average
! temperature difference between the updraft and environment
! from the LFC to the ETL. This eliminates the cases where
! VHINT is very small and the heating and cooling rates get
! astronomically large.
!
!-----------------------------------------------------------------------
!
avgmin=.10
avtdiff=0.
DO k=klfc,ketl-1
avtdiff=avtdiff+(thcon(k)-the(k))
END DO
avtdiff=avtdiff/MAX(1,ketl-klfc)
IF (avtdiff < avgmin) THEN
igo=0
RETURN
END IF
!
!-----------------------------------------------------------------------
!
! Heating and moistening rates
!
!-----------------------------------------------------------------------
!
3100 CONTINUE
DO k=2,kmt
ftcon(k)=alvl*preff*supply*vheat(k) &
/(pke(k)*rhoe(k)*vhint)
END DO
DO k=klcl,kct
frcon(k)=bkuo*supply*vmois(k)/(rhoe(k)*vmint)
END DO
DO k=2,klcl-1
! FRCON(K)=-SUPPLY*VMDRY(K)/(RHOE(K)*VDINT)
frcon(k)=0.0
END DO
!
!-----------------------------------------------------------------------
!
DO k=klfc,ketl-1
qvct1(k)=the(k)+confrq*ftcon(k)
END DO
overmax=0.
DO k=klfc,ketl-1
IF(qvct1(k)-thu(k) > overmax)THEN
overmax=(qvct1(k)-thu(k))/(ftcon(k)*confrq)
kover=k
END IF
END DO
IF(overmax > 0.) THEN
factr=1.-overmax
supply=factr*supply
GO TO 3100
END IF
cprecip=preff*supply
RETURN
END SUBROUTINE kuocp
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CP2MOD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE cp2mod(nzz) 1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Convert the cumulus heating and moistening in the convective grid
! to model grid.
!
!-----------------------------------------------------------------------
!
! AUTHOR: V. Wong, L. Zhao and X. Song
! 3/06/95
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nzz
! See the include file cumucst.inc
!
! OUTPUT:
!
! See the include file cumucst.inc
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nzz
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: k
INTEGER :: nz1,l
REAL :: tftm,frres,ftres,tfrm,tftc,s_sum,dzlft,tfrc
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'cumucst.inc'
REAL :: vctr5(nkp),vctr6(nkp)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Compute integrated heating and moistening energy tendencies
!
!-----------------------------------------------------------------------
!
DO k=2,kmt
qvct1(k)=rhoe(k)*ftcon(k)*pke(k)
qvct2(k)=rhoe(k)*alvl*frcon(k)
qvct3(k)=(zc(k)-zc(k-1))*qvct1(k) !delta (ZE(k))
qvct4(k)=(zc(k)-zc(k-1))*qvct2(k)
END DO
tftc=s_sum(kmt-1,qvct3(2),1) !vertical integration
tfrc=s_sum(kmt-1,qvct4(2),1)
!
!-----------------------------------------------------------------------
!
! Transfer energy tendencies to model grids
!
!-----------------------------------------------------------------------
!
nz1=nzz-1
DO k=1,nzz !initializing data arrays
vctr5(k)=0.
vctr6(k)=0.
END DO
dzlft=0.
l=2
DO k=2,nzz-1
IF(dzlft /= 0.) THEN
vctr5(k)=vctr5(k)+qvct1(l)*dzlft
vctr6(k)=vctr6(k)+qvct2(l)*dzlft
l=l+1
END IF
500 CONTINUE
IF(zc(l) <= zzcon(k)) THEN
vctr5(k)=vctr5(k)+qvct1(l)*(zc(l)-zc(l-1))
vctr6(k)=vctr6(k)+qvct2(l)*(zc(l)-zc(l-1))
l=l+1
dzlft=0.
GO TO 500
ELSE
vctr5(k)=vctr5(k)+qvct1(l)*(zzcon(k)-zc(l-1))
vctr6(k)=vctr6(k)+qvct2(l)*(zzcon(k)-zc(l-1))
dzlft=zc(l)-zzcon(k)
END IF
END DO
tftm=s_sum(nz1,vctr5(2),1)
tfrm=s_sum(nz1,vctr6(2),1)
!
!-----------------------------------------------------------------------
!
! Make sure the transfer from the convective grid to the model
! grid happened correctly.
!
!-----------------------------------------------------------------------
!
ftres=tftm-tftc
frres=tfrm-tfrc
IF(ABS(ftres) > .01*ABS(tftc)) THEN
PRINT*,' Energy error in grid tranfser in convective param.'
PRINT*,' TFTM,TFTC ',tftm,tftc
END IF
!
!-----------------------------------------------------------------------
!
! Change energy tendencies to temperature and mixing ratio
! tendencies.
!
!-----------------------------------------------------------------------
!
DO k=2,nzz-1
ftcon(k)=vctr5(k)/((zzcon(k+1)-zzcon(k))*dncon(k)*picon(k))
frcon(k)=vctr6(k)/((zzcon(k+1)-zzcon(k))*dncon(k)*alvl)
END DO
RETURN
END SUBROUTINE cp2mod
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ENVIRON ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE environ(nzz) 1,10
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Stability check and determination of cloud base and top
!
!-----------------------------------------------------------------------
!
! AUTHOR: V. Wong, L. Zhao and X. Song
! 3/06/95
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nzz
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: k,nkmid
REAL :: rlll,zlll,plll,themax,tlll,dzlll,tdu,rdsu,thdu,dzdd, &
abe,znz,wcpmax
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'cumucst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
REAL :: hz(nkp)
!
!-----------------------------------------------------------------------
!
! Basic constants
!
!-----------------------------------------------------------------------
!
cpr=3.4965 ! Cp/Rd
alvl=2.68E6 ! Binbin Zhuo, July 24, 1997
aliv=2.837E6 ! Lf
dzlow=200. ! higher vertical increment
dzhigh=500. ! lower vertical increment
zmid=3000. ! height of cloud base
cdzmin=3000. ! depth of cloud
!
!-----------------------------------------------------------------------
!
! Compute moist static energy profile. Cp*T+G*Z+Lv*q
!
!-----------------------------------------------------------------------
!
DO k=1,nzz
hz(k)=cp*tmpcon(k)+g*zcon(k)+alvl*qvcon(k)
END DO
!
!-----------------------------------------------------------------------
!
! Check for conditional instability. If no instability exist, there
! will be no convection.
!
!-----------------------------------------------------------------------
!
igo=0 !index variable for convection
DO k=2,nzz
IF(hz(k) > hz(k+1))THEN
igo=1
EXIT
END IF
END DO
! 500 CONTINUE
IF(igo == 0) RETURN
!
!-----------------------------------------------------------------------
!
! Check the upward motion. Only it is greater than some value, like
! WCONMIN, under ZMID, the convection is possible.
!
!-----------------------------------------------------------------------
!
igo=0
wcpmax=-1.e10
DO k=2,nzz
IF(zcon(k) > zmid)EXIT
wcpmax=MAX(wcpmax,wcon(k)) !w in cloud base
END DO
! 510 CONTINUE
IF (wcpmax > 0.0.AND.wcpmax > wcldbs ) igo=1 ! WCLDBS is input parm.
IF (igo == 0) RETURN
!
!-----------------------------------------------------------------------
!
! If the above two conditions are satisfied, set the convective analysis
! grids. There are two type of resolutions. Lower resolution for upper
! levels (above ZMID), higher resolution for lower levels (below ZMID)
!
!-----------------------------------------------------------------------
!
nkmid=zmid/dzlow+1 !# of grid points below cloud base
zc(1)=0.0
DO k=2,nkmid
zc(k)=zc(1)+(k-1)*dzlow
END DO
DO k=nkmid+1,nkp
zc(k)=zc(nkmid)+(k-nkmid)*dzhigh
END DO
!
!-----------------------------------------------------------------------
!
! Computer the height at PI grid points.
!
!-----------------------------------------------------------------------
!
ze(1)=0.
DO k=2,nkp
ze(k)=(zc(k)+zc(k-1))*.5
END DO
!
!-----------------------------------------------------------------------
!
! Find the model top on the convective analysis grids.
!
!-----------------------------------------------------------------------
!
znz=zcon(nzz) ! the height of the model top
DO k=nkp,1,-1 ! to find the cloud grid close to the model top
IF(ze(k) < znz)GO TO 520
END DO
CALL arpsstop('arpsstop called from environ',1)
! means cloud is out of the analysis domain
520 CONTINUE
kmt=k ! the level of cloud top in cloud grids
!
!-----------------------------------------------------------------------
!
! Interpolate model(environment) variables to the convective analysis
! grids.
!
!-----------------------------------------------------------------------
!
CALL htint(nzz,ucon,zcon,kmt,upe,ze) ! u
CALL htint(nzz,vcon,zcon,kmt,vpe,ze) ! v
CALL htint(nzz,wcon,zzcon,kmt,wpe,ze) ! w
CALL htint(nzz,thtcon,zcon,kmt,the,ze) ! theta
CALL htint(nzz,qvcon,zcon,kmt,rve,ze) ! q
!
!-----------------------------------------------------------------------
!
! Set the minimum value of qv.
!
!-----------------------------------------------------------------------
!
DO k=1,kmt
rve(k)=MAX(rve(k),1E-8) ! set the limited qv
END DO
!
!-----------------------------------------------------------------------
!
! Computer theta V, theta E, and get PI pressure profile.
!
!-----------------------------------------------------------------------
!
DO k=1,kmt
thve(k)=the(k)*(1.+.61*rve(k)) ! virtual potential temperature
END DO
pke(1)=picon(1)
DO k=2,kmt
pke(k)=pke(k-1)-g*2.*(ze(k)-ze(k-1)) &
/(thve(k)+thve(k-1))
END DO
DO k=1,kmt
te(k)=the(k)*pke(k)/cp ! temp
pe(k)=(pke(k)/cp)**cpr*p0 ! pressure
rhoe(k)=pe(k)/(rd*te(k)*(1.+.61*rve(k))) ! density
END DO
!
!-----------------------------------------------------------------------
!
! Computer theta E.
!
!-----------------------------------------------------------------------
!
DO k=1,kmt
CALL thetae(pe(k),te(k),rve(k),thee(k))
END DO
!
!-----------------------------------------------------------------------
!
! Find the main source level of the updraft. First, test if any
! inversion below 1.2 km.
!
!-----------------------------------------------------------------------
!
DO k=3,nkmid !check the levels below ZMID
IF(te(k) > te(k-1).AND.te(k) > te(k+1) &
.AND.ze(k) <= 1200.)THEN
kcon=k
GO TO 530
END IF
END DO
!
!-----------------------------------------------------------------------
!
! If there isn't an inversion, use the level of highest theta E.
!
!-----------------------------------------------------------------------
!
themax=0.
DO k=2,nkmid
IF(thee(k) > themax)THEN
themax=thee(k)
kcon=k
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Find the LCL of a layer average around the source level.
!
!-----------------------------------------------------------------------
!
530 CONTINUE
tlll=(te(kcon)+te(kcon+1)+te(kcon-1))/3. ! averaged
plll=pe(kcon)
rlll=(rve(kcon)+rve(kcon+1)+rve(kcon-1))/3. ! averaged
zlll=ze(kcon)
!
!-----------------------------------------------------------------------
!
! Find the (p,t) at the condesation level and the height of LCL relative
! to the source level.
!
!-----------------------------------------------------------------------
!
CALL lcl(tlll,plll,rlll,tlcl,plcl,dzlcl) ! condesation level
!
!-----------------------------------------------------------------------
!
! Find the closest level on the convective grid near the LCL.
!
!-----------------------------------------------------------------------
!
dzlll=1E20
DO k=1,kmt
dzdd=ABS(ze(k)-(zlll+dzlcl))
IF(dzdd < dzlll)THEN
dzlll=dzdd
klcl=k
END IF
END DO
!
!-----------------------------------------------------------------------
!
! If there is not upward motion at the LCL, no convection (must be
! greater than wconmin)
!
!-----------------------------------------------------------------------
!
IF(wpe(klcl) < 0.0.OR.wpe(klcl) < wcldbs)THEN
igo=0
RETURN
END IF
!
!-----------------------------------------------------------------------
!
! Locate equilibrium temperature level of an unentrained parcel.
!
!-----------------------------------------------------------------------
!
theu(klcl)=(3.376/tlcl-0.00254)*1000*rve(kcon) &
*(1+0.81*rve(kcon))
theu(klcl)=the(kcon)*EXP(theu(klcl))
DO k=klcl,kmt
IF(theu(klcl) <= thve(k))THEN
ketl=k
GO TO 540
END IF
END DO
CALL arpsstop("arpsstop called from environ at 540",1)
540 CONTINUE
!
!-----------------------------------------------------------------------
!
! If the cloud depth is less than a critical value (CDZMID), then
! convection is not allowed.
!
!-----------------------------------------------------------------------
!
IF(ze(ketl)-ze(klcl) < cdzmin)THEN
igo=0
RETURN
END IF
!
!-----------------------------------------------------------------------
!
! Computer initial ABE (Avaliable Buoyant Energy). Convection is allowed
! only if ABE is positive.
!
!-----------------------------------------------------------------------
!
abe=0.
DO k=klcl,ketl
CALL the2t(theu(klcl),pe(k),thdu,tdu,rdsu)
abe=abe+(thdu*(1.+.61*rdsu)-thve(k))/thve(k)*(zc(k)-zc(k-1))
END DO
IF(abe <= 0.)THEN
igo=0
RETURN
END IF
!
RETURN
END SUBROUTINE environ
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE LCL ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE lcl(t0,pp0,r0,tlcl,plcl,dzlcl) 1,1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Determine the LCL.
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL :: t0,pp0,r0,tlcl,plcl,dzlcl
REAL :: cpg,cpr,p00k
DATA cpg/102.45/,cpr/3.4965/,p00k/26.91535/
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: nitt
REAL :: dz,td,ttth0,ttd,rvs,r_rs,ti,pki,pppi,p0k,pi0i
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
plcl=pp0
tlcl=t0
p0k=pp0**rddcp
pi0i=p0k/p00k*cp
ttth0=t0*p00k/p0k
ttd=td(pp0,r0)
dz=cpg*(t0-ttd)
IF(dz <= 0.)THEN
dzlcl=0.
RETURN
END IF
DO nitt=1,50
pki=pi0i-g*dz/(ttth0*(1.+.61*r0))
pppi=(pki/cp)**cpr*p0
ti=ttth0*pki/cp
rvs=r_rs(pppi,ti)
IF (ABS(rvs-r0) < .00001) GO TO 110
ttd=td(pppi,r0)
dz=dz+cpg*(ti-ttd)
END DO
CALL arpsstop("arpstop called from lcl at 100",1)
110 CONTINUE
plcl=pppi
tlcl=ti
dzlcl=dz
RETURN
END SUBROUTINE lcl
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE THE2T ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE the2t(the,p,th,t,r) 3,1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Convert potential temperature to temperature.
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL :: the,p,th,t,r
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: itter
REAL :: TO,pi,alvl,tn,r_rs,cp
DATA cp/1004./,alvl/2.68E6/ ! new, suggested by G S Bhat
! and modified by Binbin Zhou
pi=(p*1E-5)**.286
TO=the/EXP((3.376/295.-0.00254)*1000*.012*(1+0.81*.012))*pi
DO itter=1,50
r=r_rs(p,TO)
th=the/EXP((3.376/TO-0.00254)*1000*r*(1+0.81*r))
tn=th*pi
IF(ABS(TO-tn) < 0.005)GO TO 12
TO=TO+(tn-TO)*.3
END DO
WRITE(6,1) the,p,TO,tn,th,r
1 FORMAT(' STOP IN ROUTINE THE2T '/' THE,P,TO,TN,TH,R',6E15.6)
CALL arpsstop("arpstop called from the2t",1)
12 CONTINUE
t=tn
RETURN
END SUBROUTINE the2t
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION R_RS ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION r_rs(p,t)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Compute saturation mixing ratio.
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL :: r_rs,p,t
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL :: es
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
es=610.78*EXP(17.269*(t-273.16)/(t-35.86))
r_rs=.622*es/(p-es)
RETURN
END FUNCTION r_rs
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION TD ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION td(p,rs)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL :: td,p,rs
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL :: rr,es,esln
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
rr=rs+1E-8
es=p*rr/(.622+rr)
esln=LOG(es)
td=(35.86*esln-4947.2325)/(esln-23.6837)
RETURN
END FUNCTION td
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE HTINT ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE htint(nzz1,vctra,eleva,nzz2,vctrb,elevb) 5,1
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nzz1,nzz2
REAL :: vctra(nzz1),vctrb(nzz2),eleva(nzz1),elevb(nzz2)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: l,k
REAL :: wt
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! stop SGI pfa compiler from parallelize the code
!*$*NOCONCURRENTIZE
l=1
DO k=1,nzz2
30 CONTINUE
IF (elevb(k) < eleva(1)) GO TO 35
IF (elevb(k) >= eleva(l).AND.elevb(k) <= eleva(l+1)) GO TO 35
IF (elevb(k) > eleva(nzz1)) GO TO 36
l=l+1
IF (l == nzz1) CALL arpsstop('arpstop called from htint',1)
GO TO 30
35 CONTINUE
wt=(elevb(k)-eleva(l))/(eleva(l+1)-eleva(l))
vctrb(k)=vctra(l)+(vctra(l+1)-vctra(l))*wt
CYCLE
36 CONTINUE
wt=(elevb(k)-eleva(nzz1))/(eleva(nzz1-1)-eleva(nzz1))
vctrb(k)=vctra(nzz1)+(vctra(nzz1-1)-vctra(nzz1))*wt
END DO
! SGI pfa compiler directive:
!*$*CONCURRENTIZE
RETURN
END SUBROUTINE htint
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE THETAE ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE thetae(p,t,rv,the) 1,1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Compute equivalent potential temperature.
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL :: p,t,rv,the
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: itter
REAL :: ttd,td,pit,tupo,tupn,tmn,dz,cpg
REAL :: cp,g,r,alvl
DATA cp/1004./,g/9.8/,r/287./,alvl/2.68E6/ ! modified by
! Binbin Zhou
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
cpg=cp/g
pit=p
tupo=t
ttd=td(p,rv)
dz=cpg*(t-ttd)
IF(dz <= 0.)GO TO 20
DO itter=1,50
tupn=t-g/cp*dz
tmn=(tupn+t)*.5*(1.+.61*rv)
pit=p*EXP(-g*dz/(r*tmn))
IF(ABS(tupn-tupo) < 0.001)GO TO 20
ttd=td(pit,rv)
tupo=tupn
dz=dz+cpg*(tupn-ttd)
END DO
CALL arpsstop("arpstop called from thetae",1)
20 CONTINUE
the=(3.376/tupo-0.00254)*1000*rv*(1+0.81*rv)
the=tupo*(1E5/pit)**.286*EXP(the)
RETURN
END SUBROUTINE thetae
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION S_SUM ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION s_sum(nn,vctr,inc)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nn,inc
REAL :: vctr(nn)
REAL :: s_sum
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: k
REAL :: sum
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
sum=0.
DO k=1,nn,inc
sum=sum+vctr(k)
END DO
s_sum=sum
RETURN
END FUNCTION s_sum
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ES_CAL ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE es_cal(t,ess)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! This function computer vapor pressure over water (in Pascal).
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL :: t,ess
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL :: a1,a2,b1,b2,tt,t1,t2,tice
DATA a1/17.2694/a2/21.87456/
DATA b1/35.86/b2/7.66/
DATA tice/260.0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
tt=t
IF (tt <= tice) GO TO 10
t1=a1*(tt-273.16)
t2=tt-b1
GO TO 11
10 t1=a2*(tt-273.16)
t2=tt-b2
11 ess=610.78*EXP(t1/t2)
RETURN
END SUBROUTINE es_cal