!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FRCUVW ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE frcuvw( nx,ny,nz,exbcbufsz,dtbig1, & 1,43
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,pbldpth, &
ubar,vbar,ptbar,pbar,ptbari,pbari,rhostr,qvbar, &
usflx,vsflx, x,y,z,zp, mapfct, &
j1,j2,j3,aj3x,aj3y,aj3z,j3inv, sinlat,ptsfc, &
uforce,vforce,wforce,kmh,kmv,rprntl,lenscl,defsq, &
exbcbuf, bcrlx,rhofct,phydro, &
tem1,tem2, &
tem3,tem4,tem5,tem6,tem7,tem8,tem9, &
tem10,tem11,tem12,tem13,tem14,tem15,tem16 )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the inactive acoustic forcing terms in the momentum
! equations. These forcing terms include the advection, mixing, Coriolis
! force and buoyancy terms, and are accumulated into one array for
! each equation, i.e.,
! uforce = - uadv + umix + ucorio.
! vforce = - vadv + vmix + vcorio.
! wforce = - wadv + wmix + wbuoy + wcorio.
! These forcing terms are held fixed during the small acoustic wave
! integration time steps. The turbulent mixing coefficient for
! momentum km is an output which will be used to calculate the
! turbulent mixing of heat and water substances.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 3/21/92.
!
! MODIFICATION HISTORY:
!
! 5/05/92 (M. Xue)
! Added full documentation.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 9/10/94 (D. Weber & Y. Lu)
! Cleaned up documentation.
!
! 01/28/95 (G. Bassett)
! Option added to turn off buoyancy terms (for buoyopt=0).
!
! 01/23/96 (Donghai Wang, Ming Xue and Yuhe Liu)
! Added the map factor to forcing terms.
!
! 4/1/96 (Donghai Wang, X. Song and M. Xue)
! Added the implicit treatment for the vertical mixing.
!
! 10/15/97 (Donghai Wang)
! Using total density (rho) in the calculation of the pressure
! gradient force terms (if rhofctopt=1).
!
! 11/05/97 (D. Weber)
! Added phydro array for use in the bottom boundary condition for
! pertrubation pressure (hydrostatic).
!
! 11/06/97 (D. Weber)
! Added three additional levels to the mapfct array. The three
! levels (4,5,6) represent the inverse of the first three in order.
! The inverse map factors are computed to improve efficiency.
!
! 11/06/97 (D. Weber)
! Added tem12 for use in turbulent mixing (optimizing tmix3d.f).
!
! 9/18/98 (D. Weber)
! Removed single operation do loops and placed the calculations
! into the subroutines that generate the forcing terms.
! Added arrays storing the average of j3 in the x,y, and z
! directions.
!
! 1999/11/24 (Limin Zhao and Alan Shapiro)
! Added the forward data assimilation package, which renews the
! the wind, temperature and pressure fields before the forcing terms
! are calculated when the assimilation mode is turned on
! (assimopt=1).
!
!-----------------------------------------------------------------------
!
! 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
! dtbig1 Local large time step size.
!
! u x component of velocity (m/s)
! v y component of velocity (m/s)
! w Vertical component of velocity in Cartesian
! coordinates (m/s).
! wcont Contravariant vertical velocity (m/s)
! ptprt Perturbation potential temperature (K)
! pprt Perturbation pressure (Pascal)
! qv Water vapor specific humidity (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rainwater mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
! tke Turbulent Kinetic Energy ((m/s)**2)
! pbldpth Planetary boundary layer depth (m)
!
! ubar Base state x velocity component (m/s)
! vbar Base state y velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! ptbari Inverse Base state potential temperature (K)
! pbari Inverse Base state pressure (Pascal)
! rhostr Base state density rhobar times j3 (kg/m**3)
! qvbar Base state water vapor specific humidity (kg/kg)
!
! usflx Surface flux of u-momentum
! vsflx Surface flux of v-momentum
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
! zp Vertical coordinate of grid points in physical space(m)
!
! mapfct Map factors at scalar, u and v points
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
! aj3x Avgx of the coordinate transformation Jacobian d(zp)/dz
! aj3y Avgy of the coordinate transformation Jacobian d(zp)/dz
! aj3z Avgz of the coordinate transformation Jacobian d(zp)/dz
! j3inv Inverse of the coordinate transformation j3 d(zp)/dz
!
! sinlat Sin of latitude at each grid point
! ptsfc Potential temperature at the ground level (K)
!
! OUTPUT:
!
! uforce Acoustically inactive forcing terms in u-momentum
! equation (kg/(m*s)**2). uforce= -uadv + umix + ucorio
! vforce Acoustically inactive forcing terms in v-momentum
! equation (kg/(m*s)**2). vforce= -vadv + vmix + vcorio
! wforce Acoustically inactive forcing terms in w-momentum
! equation (kg/(m*s)**2).
! wforce= -wadv + wmix + wbuoy + wcorio
!
! phydro Big time step forcing term for use in computing the
! hydrostatic pressure at k=1.
!
! kmh Horizontal turb. mixing coef. for momentum ( m**2/s )
! kmv Vertical turb. mixing coef. for momentum ( m**2/s )
! rprntl Reciprocal of Prandtl number
! lenscl Turbulent mixing length scale (m)
! defsq Deformation squared (1/s**2)
!
! 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.
! tem11 Temporary work array.
! tem12 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'timelvls.inc'
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: dtbig1 ! Local large time step size.
REAL :: u (nx,ny,nz,nt) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity (m/s)
REAL :: w (nx,ny,nz,nt) ! Total w-velocity (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz,nt) ! Perturbation pressure (Pascal)
REAL :: qv (nx,ny,nz,nt) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz,nt) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz,nt) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz,nt) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz,nt) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz,nt) ! Hail mixing ratio (kg/kg)
REAL :: tke (nx,ny,nz,nt) ! Turbulent kinetic energy ((m/s)**2)
REAL :: pbldpth(nx,ny,nt) ! Planetary boundary layer depth (m)
REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal).
REAL :: ptbari(nx,ny,nz) ! Inverse Base state pot. temperature (K)
REAL :: pbari (nx,ny,nz) ! Inverse Base state pressure (Pascal).
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific
! humidity(kg/kg)
REAL :: x (nx) ! x-coord. of the physical and compu-
! tational grid. Defined at u-point(m).
REAL :: y (ny) ! y-coord. of the physical and compu-
! tational grid. Defined at v-point(m).
REAL :: z (nz) ! z-coord. of the computational grid.
! Defined at w-point on the staggered
! grid(m).
REAL :: zp (nx,ny,nz) ! Physical height coordinate defined at
! w-point of the staggered grid(m).
REAL :: mapfct(nx,ny,8) ! Map factors at scalar, u and v points
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/d(x)
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/d(y)
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian
! d(zp)/d(z)
REAL :: aj3x (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
REAL :: aj3y (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
REAL :: aj3z (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
REAL :: j3inv (nx,ny,nz) ! Coordinate transformation Jacobian
! d(zp)/d(z)
REAL :: sinlat(nx,ny) ! Sin of latitude at each grid point
REAL :: ptsfc (nx,ny) ! Potential temperature at the ground level (K)
REAL :: usflx (nx,ny) ! Surface flux of u-momentum
! (kg/(m*s**2))
REAL :: vsflx (nx,ny) ! Surface flux of v-momentum
! (kg/(m*s**2))
REAL :: uforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in u-momentum equation (kg/(m*s)**2)
! uforce= -uadv + umix + ucorio
REAL :: vforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in v-momentum equation (kg/(m*s)**2)
! vforce= -vadv + vmix + vcorio
REAL :: wforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in w-momentum equation (kg/(m*s)**2)
! wforce= -wadv + wmix + wbuoy
REAL :: phydro(nx,ny) ! Big time step forcing for computing
! hydrostatic pprt at k=1.
REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: rprntl(nx,ny,nz) ! Reciprocal of Prandtl number
REAL :: lenscl(nx,ny,nz) ! Turbulent mixing length scale (m)
REAL :: defsq (nx,ny,nz) ! Deformation squared (1/s**2)
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: bcrlx(nx,ny) ! EXBC relaxation coefficients
REAL :: rhofct(nx,ny,nz) ! rho-factor:rhobar/rho
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.
REAL :: tem11 (nx,ny,nz) ! Temporary work array.
REAL :: tem12 (nx,ny,nz) ! Temporary work array.
REAL :: tem13 (nx,ny,nz) ! Temporary work array.
REAL :: tem14 (nx,ny,nz) ! Temporary work array.
REAL :: tem15 (nx,ny,nz) ! Temporary work array.
REAL :: tem16 (nx,ny,nz) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
INTEGER :: tlevel
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global constants that control model
! execution
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
INCLUDE 'assim.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Calculate the total mixing (which includes subgrid scale
! turbulent mixing and computational mixing) terms for u, v, w
! equations as well as the mixing coefficient kmh and kmv, and
! the inverse Prandtl number.
!
! These mixing terms are accumulated in the arrays
! uforce, vforce and wforce.
!
! Since all mixing terms are evaluated at the past time level,
! we pass only the variable fields at time tpast into the routine.
!
!-----------------------------------------------------------------------
!
IF ( assimopt == 1 ) THEN
!
!-----------------------------------------------------------------------
!
! ASSIMCON is a driver for the whole data assimilation system, which
! determines the control parameters and switches in the assimilation
! system.
!
!------------------------------------------------------------------------
!
CALL assimcon
(nx,ny,nz,x,y,z,zp, &
ubar,vbar,pbar,ptbar,rhostr,qvbar, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
j1,j2,j3, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9, &
tem10,tem11)
! CALL assimcon
!
!------------------------------------------------------------------------
!
! ASSIMDRIV is a driver for the variational velocity adjustment,which
! controls the ways to perform the varaitional adjustment and the
! velocity blending.
!
! Note: the arrays uforce, vforce, wforce, defsq are used in ASSIMDRIV
! as working array. This will not effect any model results since
! these array will be specified after the data assimilation.
!
!------------------------------------------------------------------------
!
CALL assimdriv(nx,ny,nz,x,y,z,zp, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
tke,kmh,kmv,mapfct, &
ubar,vbar,ptbar,pbar,rhostr,qvbar,j1,j2,j3, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9, &
tem10,tem11,tem12, &
tem13,tem14,tem15,tem16)
! CALL assimdriv
END IF
!
!------------------------------------------------------------------------
!
! If the assim mode is on, the modified model variables will be used
! to calculate the force terms in the following subroutines.
!
!
! Calculate the total mixing (which includes subgrid scale
! turbulent mixing and computational mixing) terms for u, v, w
! equations as well as the mixing coefficient kmh and kmv, and
! the inverse Prandtl number.
!
! These mixing terms are accumulated in the arrays
! uforce, vforce and wforce.
!
! Since all mixing terms are evaluated at the past time level,
! we pass only the variable fields at time tpast into the routine.
!
!-----------------------------------------------------------------------
!
tlevel = tpast
!call test_dump (u(1,1,1,tlevel),"XXXBmixuvw_u",nx,ny,nz,1,1)
!call test_dump (v(1,1,1,tlevel),"XXXBmixuvw_v",nx,ny,nz,2,1)
!call test_dump (w(1,1,1,tlevel),"XXXBmixuvw_w",nx,ny,nz,3,1)
!call test_dump (ptprt(1,1,1,tlevel),"XXXBmixuvw_ptprt",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,tlevel),"XXXBmixuvw_pprt",nx,ny,nz,0,1)
!call test_dump (qv(1,1,1,tlevel),"XXXBmixuvw_qv",nx,ny,nz,0,1)
!call test_dump (qc(1,1,1,tlevel),"XXXBmixuvw_qc",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tlevel),"XXXBmixuvw_qr",nx,ny,nz,0,1)
!call test_dump (qi(1,1,1,tlevel),"XXXBmixuvw_qi",nx,ny,nz,0,1)
!call test_dump (qs(1,1,1,tlevel),"XXXBmixuvw_qs",nx,ny,nz,0,1)
!call test_dump (qh(1,1,1,tlevel),"XXXBmixuvw_qh",nx,ny,nz,0,1)
!call test_dump (tke,"XXXBmixuvwtke",nx,ny,nz,0,0)
!call test_dump (uforce,"XXXBmixuvw_uforce",nx,ny,nz,1,0)
CALL mixuvw(nx,ny,nz, exbcbufsz, &
u (1,1,1,tlevel),v (1,1,1,tlevel), &
w (1,1,1,tlevel), &
ptprt(1,1,1,tlevel),pprt(1,1,1,tlevel), &
qv (1,1,1,tlevel),qc (1,1,1,tlevel), &
qr (1,1,1,tlevel), &
qi (1,1,1,tlevel),qs (1,1,1,tlevel), &
qh (1,1,1,tlevel), &
tke, pbldpth(1,1,tpresent), &
ubar,vbar,ptbar, &
pbar,rhostr,qvbar, &
usflx,vsflx, x,y,z,zp,mapfct, &
j1,j2,j3,aj3x,aj3y,aj3z,j3inv, ptsfc, &
uforce,vforce,wforce,kmh,kmv,rprntl,lenscl,defsq, &
exbcbuf, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9,tem10, &
tem11,tem12)
!call test_dump (uforce,"XXXAmixuvw_uforce",nx,ny,nz,1,0)
IF( (tmixopt /= 0 .OR. cmix2nd /= 0 .OR. cmix4th /= 0 ) &
.AND. lvldbg >= 4 ) THEN
CALL set_acct(misc_acct)
CALL checkuhx(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'umixx', tem1)
CALL checkuhy(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'umixy', tem1)
CALL checkvhx(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vmixx', tem1)
CALL checkvhy(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vmixy', tem1)
CALL checkwhx(wforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wmixx', tem1)
CALL checkwhy(wforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wmixy', tem1)
END IF
!
!
!-----------------------------------------------------------------------
!
! Calculate the advection terms in the momentum equations, using
! the equivalent advective formulation.
!
! On exit of advuvw, the advection terms are included in the
! forcing arrays which contain the turbulent mixing terms.
!
!-----------------------------------------------------------------------
!
!
CALL set_acct(advuvw_acct)
CALL advuvw(nx,ny,nz,u,v,w,wcont,rhostr,ubar,vbar, mapfct, &
uforce,vforce,wforce, &
tem1,tem2,tem3, &
tem4,tem5,tem6,tem7,tem8,tem9)
!call test_dump (uforce,"XXXAadvuvw_uforce",nx,ny,nz,1,0)
IF( lvldbg >= 4 ) THEN
CALL checkuhx(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'uforce after advu', tem9)
CALL checkuhy(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'uforce after advu', tem9)
CALL checkvhx(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vforce after advv', tem9)
CALL checkvhy(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vforce after advv', tem9)
CALL checkwhx(wforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wforce after advw', tem9)
CALL checkwhy(wforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wforce after advw', tem9)
END IF
!-----------------------------------------------------------------------
!
! Calculate the advection terms from inhomogeneous map factor for
! u and v, and store them into tem1 and tem2.
!
! Ustr and vstr were calculated by advuvw and stored in tem4 and
! tem5.
!
!-----------------------------------------------------------------------
IF ( mptrmopt /= 0 ) THEN
CALL maptrm( nx,ny,nz, u,v, tem4,tem5, mapfct, &
uforce, vforce, tem7,tem8,tem9 )
!call test_dump (uforce,"XXXAmaptrm_uforce",nx,ny,nz,1,0)
IF ( lvldbg >= 5 ) THEN
CALL checkuhx(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'uforce after maptrm', tem9)
CALL checkuhy(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'uforce after maptrm', tem9)
CALL checkvhx(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vforce after maptrm', tem9)
CALL checkvhy(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vforce after maptrm', tem9)
END IF
END IF
!-----------------------------------------------------------------------
!
! Calculate the Coriolis terms in u, v and w equations and add
! to the forcing terms INSIDE the call to coriol.
!
!-----------------------------------------------------------------------
IF( coriopt /= 0 ) THEN
CALL set_acct(coriol_acct)
tlevel = tpresent
CALL coriol(nx,ny,nz, &
u(1,1,1,tlevel),v(1,1,1,tlevel),w(1,1,1,tlevel), &
ubar,vbar,rhostr, sinlat, &
uforce, vforce, wforce, tem4,tem5,tem6,tem7)
!call test_dump (uforce,"XXXAcoriol_uforce",nx,ny,nz,1,0)
IF( lvldbg >= 4 ) THEN
CALL checkuhx(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'uforce after ucorx', tem9)
CALL checkuhy(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'uforce after ucory', tem9)
CALL checkvhx(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vforce after vcorx', tem9)
CALL checkvhy(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vforce after vcory', tem9)
CALL checkwhx(wforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wforce after wcorx', tem9)
CALL checkwhy(wforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wforce after wcory', tem9)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! The temperature and pressure fields are modified by Gal-Chen's
! thermodynamic recovery method if the assim mode is turned on.
!
!-----------------------------------------------------------------------
!
IF ( assimopt == 1 ) THEN
CALL set_acct(tinteg_acct)
CALL assimptpr(nx,ny,nz,x,y,z,zp, &
u,v,w,ptprt,pprt,qv, &
qc,qr,qi,qs,qh, &
tke,kmh,kmv, &
ubar,vbar,ptbar,pbar,qvbar,rhostr, &
uforce,vforce,wforce,j1,j2,j3)
! CALL assimptpr
WRITE(6,*)'back from assimptpr'
END IF
!
!-----------------------------------------------------------------------
!
! Calculate the rho-factor: rhobar/rho, to correct the calculation
! for the pressure gradient terms
!
!-----------------------------------------------------------------------
!
CALL set_acct(buoy_acct)
IF (rhofctopt /= 0) THEN
tlevel = tpresent
!call test_dump2(qr(1,1,1,tlevel),"XXXBrhofactor_qr",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
CALL rhofactor(nx,ny,nz,ptprt(1,1,1,tlevel),pprt(1,1,1,tlevel), &
qv(1,1,1,tlevel),qc(1,1,1,tlevel),qr(1,1,1,tlevel), &
qi(1,1,1,tlevel),qs(1,1,1,tlevel),qh(1,1,1,tlevel), &
ptbar,pbar,ptbari,pbari,qvbar,rhofct,tem1)
!call test_dump2(rhofct,"XXXArhofactor_rhofct",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
ELSE
DO k=1,nz
DO j=1,ny
DO i=1,nx
rhofct(i,j,k) = 1.0
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Average rhofct at w points, stored in tem3
!
!-----------------------------------------------------------------------
!
CALL avgsw(rhofct,nx,ny,nz, 1,nx-1, 1,ny-1,tem3)
!call test_dump2(rhofct,"XXXAavgsw_rhofct",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!
!-----------------------------------------------------------------------
!
! Calculate the buoyancy term for the w-equation.
! The buoyancy term is stored in tem1 then added to wforce.
!
! If buoyopt = 0 then buoyancy is turned off.
!
!-----------------------------------------------------------------------
!
tlevel = tpresent
IF ( buoyopt /= 0 ) THEN
CALL buoycy(nx,ny,nz, ptprt(1,1,1,tlevel),pprt(1,1,1,tlevel), &
qv(1,1,1,tlevel),qc(1,1,1,tlevel),qr(1,1,1,tlevel), &
qi(1,1,1,tlevel),qs(1,1,1,tlevel),qh(1,1,1,tlevel), &
ptbar,pbar,ptbari,pbari,rhostr,qvbar, &
tem1, & ! wbuoy
tem2)
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
wforce(i,j,k)=wforce(i,j,k)+tem3(i,j,k)*tem1(i,j,k)
END DO
END DO
END DO
DO j=1,ny-1 ! store buoyancy at k=2 into phydro.
DO i=1,nx-1 ! for use in pprt(i,j,1) boundary condition
phydro(i,j)=tem1(i,j,2)
END DO
END DO
END IF
IF( lvldbg >= 4 ) THEN
CALL checkwhx(tem1, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wbuox', tem9)
CALL checkwhy(tem1, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wbuoy', tem9)
END IF
!
!-----------------------------------------------------------------------
!
! To calculate additional boundary relaxation and mixing terms for
! the wind fields in the case of externally forced boundary condition.
!
!-----------------------------------------------------------------------
!
IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN
CALL set_acct(cmix_acct)
CALL brlxuvw( nx,ny,nz, dtbig1, &
u(1,1,1,1),v(1,1,1,1),w(1,1,1,1),rhostr, &
uforce,vforce,wforce, &
exbcbuf(nu0exb), exbcbuf(nv0exb), &
exbcbuf(nw0exb), exbcbuf(nudtexb), &
exbcbuf(nvdtexb),exbcbuf(nwdtexb),bcrlx, &
tem1,tem2,tem3,tem4,tem5,tem6 )
!call test_dump (uforce,"XXXAbrlxuvw_uforce",nx,ny,nz,1,0)
END IF
!
!-----------------------------------------------------------------------
!
! To calculate the vertically implicit mixing terms,
! and add to uforce, vforce, and wforce.
!
!-----------------------------------------------------------------------
!
IF (trbvimp == 1) THEN ! Vertical implicit application
CALL set_acct(tmix_acct)
CALL vmiximpuvw(nx,ny,nz,dtbig1,u(1,1,1,1),v(1,1,1,1),w(1,1,1,1), &
rhostr,kmv,j1,j2,j3inv, &
uforce,vforce,wforce, &
tem1,tem2,tem3,tem4, &
tem5,tem6,tem7,tem8,tem9,tem10,tem11)
!call test_dump (uforce,"XXXAvmiximp_uforce",nx,ny,nz,1,0)
END IF
RETURN
END SUBROUTINE frcuvw
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FRCP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE frcp( nx,ny,nz,exbcbufsz, dtbig1, & 1,6
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
ptbar,pbar,rhostr,qvbar,mapfct,j1,j2,j3,aj3x,aj3y,aj3z, &
pforce, &
exbcbuf, bcrlx, &
padv,tem1,tem2,tem3,tem4,tem5,tem6,tem7,mp_tem)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the total acoustically inactive forcing terms in
! the pressure equation (advection and other source/sink terms).
!
! These terms are invariant during the small time-step integration.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 3/21/92.
!
! MODIFICATION HISTORY:
!
! 5/05/92 (M. Xue)
! Added full documentation.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 9/10/94 (D. Weber & Y. Lu)
! Cleaned up documentation.
!
! 01/23/96 (Donghai Wang, Yuhe Liu and Ming Xue)
! Added the map factor to forcing terms.
!
! 9/18/98 (D. Weber)
! Added arrays aj3x,aj3y,aj3z.
!
!-----------------------------------------------------------------------
!
! 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 direction.
! dtbig1 Local large time step size.
!
! u x component of velocity (m/s)
! v y component of velocity (m/s)
! w Vertical component of velocity in Cartesian
! coordinates (m/s).
! wcont Contravariant vertical velocity (m/s)
! ptprt Perturbation potential temperature (K)
! pprt Perturbation pressure (Pascal)
! qv Water vapor specific humidity (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rainwater mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
!
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! rhostr Base state density rhobar times j3 (kg/m**3)
! qvbar Base state water vapor specific humidity (kg/kg)
!
! mapfct Map factors at scalar points
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
! aj3x Avgx of the coordinate transformation Jacobian d(zp)/dz
! aj3y Avgy of the coordinate transformation Jacobian d(zp)/dz
! aj3z Avgz of the coordinate transformation Jacobian d(zp)/dz
!
! OUTPUT:
!
! pforce Acoustically inactive forcing terms in pressure
! equation (Pascal/s). pforce= -padv
!
! WORK ARRAYS:
!
! padv Advection of pressure
! 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.
! mp_tem Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'timelvls.inc'
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: dtbig1 ! Local large time step size.
REAL :: u (nx,ny,nz,nt) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity (m/s)
REAL :: w (nx,ny,nz,nt) ! Total w-velocity (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature
! from that of base state atmosphere
! (Kelvin)
REAL :: pprt (nx,ny,nz,nt) ! Perturbation pressure from that
! of base state atmosphere (Pascal)
REAL :: qv (nx,ny,nz,nt) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz,nt) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz,nt) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz,nt) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz,nt) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz,nt) ! Hail mixing ratio (kg/kg)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal).
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific
! humidity(kg/kg)
REAL :: mapfct(nx,ny) ! Map factors at scalar points
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dx.
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dy.
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian
! d(zp)/dz.
REAL :: aj3x (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
REAL :: aj3y (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
REAL :: aj3z (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
REAL :: pforce(nx,ny,nz) ! Sound-wave independent forcing terms
! in pressure equation (Pascal/s).
! pforce = -padv
REAL :: padv (nx,ny,nz) ! The advection term of pressure eq.
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: bcrlx(nx,ny) ! EXBC relaxation coefficients
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 :: mp_tem(MAX(nx+1,ny+1)*(nz+1)) ! Temporary message passing array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k ! local varaibles.
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global constants that control model
! execution
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Calculate the advection term in the pressure equation
! using the equivalent advective formulation
!
!-----------------------------------------------------------------------
!
!
CALL set_acct(advs_acct)
IF( peqopt == 1 ) THEN
CALL advp(nx,ny,nz,pprt,u,v,w,wcont,rhostr,mapfct, &
j3,aj3x,aj3y,aj3z, &
padv, tem1,tem2,tem3,tem4,tem5,tem6,tem7,mp_tem)
IF( lvldbg >= 4) THEN
CALL checkshx(padv, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2, &
'padvx', tem3)
CALL checkshy(padv, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2, &
'padvy', tem3)
END IF
!
!-----------------------------------------------------------------------
!
! Store the total forcing in array pforce.
!
!-----------------------------------------------------------------------
!
!
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
pforce(i,j,k)= -padv(i,j,k)
END DO
END DO
END DO
!call test_dump (pforce,"XXXIfrcp_padv_pforce",nx,ny,nz,0,0)
ELSE
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
pforce(i,j,k)= 0.0
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! To calculate additional boundary relaxation and mixing terms for
! pressure in the case of externally forced boundary condition.
!
!-----------------------------------------------------------------------
!
IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN
CALL set_acct(bc_acct)
CALL brlxp(nx,ny,nz, dtbig1, pprt(1,1,1,1),rhostr, pforce, &
exbcbuf(npr0exb),exbcbuf(nprdtexb), bcrlx, &
tem1,tem2,tem3,tem4)
!call test_dump (pforce,"XXXAbrlxp_pforce",nx,ny,nz,0,0)
END IF
RETURN
END SUBROUTINE frcp
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FRCPT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE frcpt(nx,ny,nz, exbcbufsz, dtbig1,ptprt,u,v,w,wcont, & 1,9
ptbar,rhostr,rhostri,kmh,kmv,rprntl, &
usflx,vsflx,ptsflx,pbldpth, &
x,y,z,zp,mapfct,j1,j2,j3,aj3x,aj3y,j3inv,ptsfc, &
ptforce, &
exbcbuf, bcrlx, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
tem8,tem9,tem10,tem11, &
tem1_0,tem2_0,tem3_0,mp_tem)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate gravity wave or inactive acoustic wave terms in the
! potential temperature equation.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 9/17/94
!
! MODIFICATION HISTORY:
!
! 8/15/95 (Ming Xue)
! Corrected a bug related the calculation of addtional boundary zone
! relaxation. The ptmix term was effectively added twice to
! ptforce when lbcopt=2.
!
! 01/23/96 (Donghai Wang, Yuhe Liu and Ming Xue)
! Added the map factor to forcing terms.
!
! 4/1/96 (Donghai Wang, X. Song and M. Xue)
! Added the implicit treatment for the vertical mixing.
!
! 7/10/1997 (Fanyou Kong -- CMRP)
! Added the positive definite advection option (sadvopt = 5)
!
! 7/17/1998 (M. Xue)
! Changed call to ADVPTFCT.
!
! 9/18/1998 (D. Weber)
! Added arrays aj3x,y.
!
!-----------------------------------------------------------------------
!
! 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
!
! dtbig1 The big time step size (s)
!
! ptprt Perturbation potential temperature at times tpast and
! tpresent (K)
!
! u x component of velocity at all time levels (m/s)
! v y component of velocity at all time levels (m/s)
! w Vertical component of Cartesian velocity at times
! tpast and tpresent (m/s)
! wcont Contravariant vertical velocity (m/s)
!
! ptbar Base state potential temperature (K)
! rhostr Base state density rhobar times j3 (kg/m**3)
! rhostri Inv. base state density rhobar times j3 (kg/m**3)
!
! kmh Horizontal turb. mixing coef. for momentum ( m**2/s )
! kmv Vertical turb. mixing coef. for momentum ( m**2/s )
! rprntl Reciprocal of Prandtl number
!
! ptsflx Surface flux of heat (K*kg/(m**2*s)).
! pbldpth Planetary boundary layer depth (m)
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
! zp Vertical coordinate of grid points in physical space(m)
!
! mapfct Map factors at scalar points
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
! aj3x Avgx of the coordinate transformation Jacobian d(zp)/dz
! aj3y Avgy of the coordinate transformation Jacobian d(zp)/dz
!
! OUTPUT:
!
! ptforce Gravity wave inactive forcing terms in pt-eq.
! (K*kg/(m**3*s))
!
! 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.
! tem11 Temporary work array.
!
! tem1_0 Temporary work array.
! tem2_0 Temporary work array.
! tem3_0 Temporary work array.
!
! mp_tem Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INCLUDE 'timelvls.inc'
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: dtbig1 ! Local big time step size
REAL :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: u (nx,ny,nz,nt) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity (m/s)
REAL :: w (nx,ny,nz,nt) ! Total w-velocity (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: rhostri(nx,ny,nz) ! Inv. base state density rhobar times j3.
REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: rprntl(nx,ny,nz) ! Reciprocal of Prandtl number
REAL :: usflx (nx,ny) ! Surface flux of u-momentum
! (kg/(m*s**2))
REAL :: vsflx (nx,ny) ! Surface flux of v-momentum
! (kg/(m*s**2))
REAL :: ptsflx(nx,ny) ! Surface flux of heat (K*kg/(m**2*s))
REAL :: pbldpth(nx,ny,nt) ! Planetary boundary layer depth (m)
REAL :: x (nx) ! x-coord. of the physical and compu-
! tational grid. Defined at u-point.
REAL :: y (ny) ! y-coord. of the physical and compu-
! tational grid. Defined at v-point.
REAL :: z (nz) ! z-coord. of the computational grid.
! Defined at w-point on the staggered
! grid.
REAL :: zp (nx,ny,nz) ! Physical height coordinate defined at
! w-point of the staggered grid.
REAL :: mapfct(nx,ny,8) ! Map factors at scalar points
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian
! defined as - d( zp )/d( x ).
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian
! defined as - d( zp )/d( y ).
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian
! defined as d( zp )/d( z ).
REAL :: aj3x (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
REAL :: aj3y (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
REAL :: j3inv (nx,ny,nz) ! Coordinate transformation Jacobian
! defined as d( zp )/d( z ).
REAL :: ptsfc (nx,ny) ! Ground surface potential temperature (K)
REAL :: ptforce(nx,ny,nz) ! Gravity wave inactive forcing terms
! in pressure equation (Pascal/s)
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: bcrlx(nx,ny) ! EXBC relaxation coefficients
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
REAL :: tem11 (nx,ny,nz) ! Temporary work array
REAL :: tem1_0(0:nx,0:ny,0:nz) ! Temporary work array.
REAL :: tem2_0(0:nx,0:ny,0:nz) ! Temporary work array.
REAL :: tem3_0(0:nx,0:ny,0:nz) ! Temporary work array.
REAL :: mp_tem(MAX(nx+1,ny+1)*(nz+1)) ! Temporary message passing array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k, tstrtlvl
REAL :: deltat
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Compute the advection term of the potential temperature
! equation and store it in array ptadv.
!
!-----------------------------------------------------------------------
!
CALL set_acct(advs_acct)
IF( sadvopt == 4 ) THEN ! FCT advection
CALL advptfct(nx,ny,nz,dtbig1,ptprt,u,v,w,wcont, &
rhostr,rhostri,ptbar,mapfct,j3,j3inv, &
ptforce, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem9,tem10,tem11,tem1_0,tem2_0,tem3_0,mp_tem)
deltat = dtbig1
tstrtlvl = tpresent
ELSE
CALL advpt(nx,ny,nz,ptprt,u,v,w,wcont, rhostr,ptbar,mapfct, &
j3,j3inv,ptforce, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)
deltat = dtbig1*2
tstrtlvl = tpast
END IF
!
!-----------------------------------------------------------------------
!
! Compute the mixing terms in the potential temperature equation.
! This includes both physical and computational mixing.
! Store in array ptmix.
!
!-----------------------------------------------------------------------
!
CALL set_acct(cmix_acct)
CALL mixpt(nx,ny,nz, exbcbufsz, &
ptprt(1,1,1,tstrtlvl),ptbar,rhostr, &
kmh,kmv,rprntl, &
usflx,vsflx,ptsflx,pbldpth(1,1,tpresent), &
x,y,z,zp,mapfct,j1,j2,j3,aj3x,aj3y,j3inv,ptsfc, &
tem7, exbcbuf, &
tem1,tem2,tem3,tem4,tem5,tem6)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
ptforce(i,j,k)=-ptforce(i,j,k)+tem7(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate additional relaxation and mixing on ptprt for the
! base grid when external boundary forcing is used.
!
!-----------------------------------------------------------------------
!
IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN
CALL set_acct(bc_acct)
CALL brlxpt(nx,ny,nz,deltat*0.5,ptprt(1,1,1,tstrtlvl), &
rhostr,ptforce, &
exbcbuf(npt0exb),exbcbuf(nptdtexb),bcrlx, &
tem1,tem2,tem3,tem4)
END IF
!
!-----------------------------------------------------------------------
!
! Treat the vertically implicit mixing
!
!-----------------------------------------------------------------------
!
IF (trbvimp == 1) THEN ! Vertical implicit application
CALL set_acct(tmix_acct)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=ptbar(i,j,k)+ptprt(i,j,k,tstrtlvl)
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)=rhostr(i,j,k)*kmv(i,j,k)*rprntl(i,j,k) &
*j3inv(i,j,k)*j3inv(i,j,k)
END DO
END DO
END DO
CALL vmiximps(nx,ny,nz,deltat*0.5,tem1,rhostr,tem2, &
ptforce,tem3,tem4,tem5,tem6)
END IF
RETURN
END SUBROUTINE frcpt
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CORIOL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE coriol(nx,ny,nz, & 1
u,v,w,ubar,vbar,rhostr, sinlat, uforce, vforce, wforce, &
fcorio1, fcorio2, tem1,tem2)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the Coriolis force terms in the u, v and w equations.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 4/20/92.
!
! MODIFICATION HISTORY:
!
! 5/05/92 (M. Xue)
! Added full documentation.
!
! 8/28/92 (M. Xue)
! Included the domain translation effect into the Coriolis force
! terms.
!
! 9/10/94 (D. Weber & Y. Lu)
! Cleaned up documentation.
!
! 9/14/98 (D. Weber)
! Removed operators and merged loops and added coriolis terms
! to the forcing terms inside this subroutine.
!
!-----------------------------------------------------------------------
!
! 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 direction.
!
! 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)
! ubar Base state x velocity component (m/s)
! vbar Base state y velocity component (m/s)
! rhostr Base state density rhobar times j3 (kg/m**3)
! sinlat Sin of latitude at each grid point
!
! OUTPUT:
!
! uforce Acoustically inactive forcing terms in u-momentum
! equation (kg/(m*s)**2). uforce= -uadv + umix + ucorio
! vforce Acoustically inactive forcing terms in v-momentum
! equation (kg/(m*s)**2). vforce= -vadv + vmix + vcorio
! wforce Acoustically inactive forcing terms in w-momentum
! equation (kg/(m*s)**2).
!
! WORK ARRAYS:
!
! fcorio1 Temporary work array.
! fcorio2 Temporary work array.
! tem1 Temporary work array.
! tem2 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: u (nx,ny,nz) ! Total u-velocity at a time level (m/s)
REAL :: v (nx,ny,nz) ! Total v-velocity at a time level (m/s)
REAL :: w (nx,ny,nz) ! Total w-velocity at a time level (m/s)
REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s)
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: sinlat(nx,ny) ! Sin of latitude at each grid point
REAL :: uforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in u-momentum equation (kg/(m*s)**2)
! uforce= -uadv + umix + ucorio
REAL :: vforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in v-momentum equation (kg/(m*s)**2)
! vforce= -vadv + vmix + vcorio
REAL :: wforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in w-momentum equation (kg/(m*s)**2)
! wforce= -wadv + wmix + wbuoy
REAL :: fcorio1 (nx,ny) ! Temporary work array.
REAL :: fcorio2 (nx,ny) ! Temporary work array.
REAL :: tem1 (nx,ny,nz) ! Temporary work array.
REAL :: tem2 (nx,ny,nz) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
INTEGER :: onvf
REAL :: omega2,sinclat,cosclat
REAL :: gumove,gvmove
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF ( grdtrns == 0 ) THEN
gumove = 0.0
gvmove = 0.0
ELSE
gumove = umove
gvmove = vmove
END IF
omega2 = 2.0* omega
IF( coriopt == 1 ) THEN
sinclat = SIN( ATAN(1.0)/45.0 * ctrlat )
DO j=1,ny
DO i=1,nx
fcorio1(i,j) = omega2* sinclat
fcorio2(i,j) = 0.0
END DO
END DO
ELSE IF( coriopt == 2 ) THEN
sinclat = SIN( ATAN(1.0)/45.0 * ctrlat )
cosclat = SQRT( 1.0 - sinclat*sinclat )
DO j=1,ny
DO i=1,nx
fcorio1(i,j) = omega2* sinclat
fcorio2(i,j) = omega2* cosclat
END DO
END DO
ELSE IF( coriopt == 3 ) THEN
DO j=1,ny-1
DO i=1,nx-1
fcorio1(i,j) = omega2*sinlat(i,j)
fcorio2(i,j) = 0.0
END DO
END DO
ELSE IF( coriopt == 4 ) THEN
DO j=1,ny-1
DO i=1,nx-1
fcorio1(i,j) = omega2* sinlat(i,j)
fcorio2(i,j) = omega2* SQRT( 1.0 - sinlat(i,j)**2 )
END DO
END DO
END IF
!-----------------------------------------------------------------------
!
! Coriolis terms in u-eq. if coriotrm=1
!
! ucorio = avgx( rhostr * (fcorio1*avgy(v)-fcorio2*avgz(w)) )
!
! or if coriotrm=2,
!
! ucorio = avgx( rhostr * (fcorio1*avgy(v-vbar)-fcorio2*avgz(w)) )
!
! where fcorio1 = 2*omega*sinlat.
! and fcorio2 = 2*omega* sqrt(1-sinlat**2).
!
!-----------------------------------------------------------------------
IF( coriotrm == 1 ) THEN
DO k=1,nz-1
DO j=1,ny
DO i=1,nx-1
tem2(i,j,k) = v(i,j,k) + gvmove
END DO
END DO
END DO
ELSE
DO k=1,nz-1
DO j=1,ny
DO i=1,nx-1
tem2(i,j,k) = v(i,j,k) + gvmove - vbar(i,j,k)
END DO
END DO
END DO
END IF
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = 0.5*rhostr(i,j,k)* &
(fcorio1(i,j)*(tem2(i,j+1,k)+tem2(i,j,k)) &
-fcorio2(i,j)*(w(i,j,k+1)+w(i,j,k)))
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=2,nx-1
uforce(i,j,k) = uforce(i,j,k) + 0.5*(tem1(i-1,j,k)+tem1(i,j,k))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Coriolis terms in v-eq. and w-eq. if coriotrm=1:
!
! vcorio = - avgy( rhostr * fcorio1 * avgx(u) )
! wcorio = avgz( rhostr * fcorio2 * avgx(u) )
!
! or if coriotrm=2:
!
! vcorio = - avgy( rhostr * fcorio1 * avgx(u-ubar) )
! wcorio = avgz( rhostr * fcorio2 * avgx(u) )
!
!
!-----------------------------------------------------------------------
IF( coriotrm == 1 ) THEN
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx
tem2(i,j,k) = u(i,j,k) + gumove
END DO
END DO
END DO
ELSE
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx
tem2(i,j,k) = u(i,j,k) + gumove -ubar(i,j,k)
END DO
END DO
END DO
END IF
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = - rhostr(i,j,k)*fcorio1(i,j)* &
0.5*(tem2(i+1,j,k)+tem2(i,j,k))
END DO
END DO
END DO
DO k=1,nz-1
DO j=2,ny-1
DO i=1,nx-1
vforce(i,j,k) = vforce(i,j,k)+0.5*(tem1(i,j-1,k)+tem1(i,j,k))
END DO
END DO
END DO
IF( coriopt == 2 .OR. coriopt == 4 ) THEN ! compute the w term.
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = fcorio2(i,j) * rhostr(i,j,k)* &
( 0.5*(u(i+1,j,k)+u(i,j,k)) + gumove )
END DO
END DO
END DO
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
wforce(i,j,k) = wforce(i,j,k)+0.5*(tem1(i,j,k-1)+tem1(i,j,k))
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE coriol
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE BUOYCY ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE buoycy(nx,ny,nz,ptprt,pprt,qv,qc,qr,qi,qs,qh, & 1
ptbar,pbar,ptbari,pbari,rhostr,qvbar, wbuoy, tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the total buoyancy including liquid and solid water
! loading.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/91.
!
! MODIFICATION HISTORY:
!
! 5/05/92 (M. Xue)
! Added full documentation.
!
! 3/10/93 (M. Xue)
! The buoyancy term is reformulated. The previous formula was
! in error. The water loading was calculated wrong, resulting in
! a value of the water loading that is typically an order of
! magnitude too small.
!
! 3/25/94 (G. Bassett & M. Xue)
! The buoyancy terms are reformulated for better numerical accuracy.
! Instead of storing numbers which had the form (1+eps)*(1+eps1)
! (eps << 1 and eps1 <<1), terms were expanded out, and most of the
! high order terms neglected, except for the second order terms
! in ptprt, pprt and qvbar.
!
! 9/10/94 (D. Weber & Y. Lu)
! Cleaned up documentation.
!
! 6/21/95 (Alan Shapiro)
! Fixed bug involving missing qvpert term in buoyancy formulation.
!
! 10/15/97 (Donghai Wang)
! Added a new option for including the second order terms.
!
! 11/05/97 (D. Weber)
! Changed lower loop bounds in DO LOOP 400 for computing the
! buoyancy term from k=3,nz-2 to k=2,nz-1. Level k=2 data will be
! used in the hydrostatic pprt lower boundary condition (removed
! DO LOOP 410 used to set wbuoy = 0.0 for k= 2 and nz-1).
!
!-----------------------------------------------------------------------
!
! 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 direction.
!
! ptprt Perturbation potential temperature at a time level (K)
! pprt Perturbation pressure at a given time level (Pascal)
! qv Water vapor specific humidity at a given time level
! (kg/kg)
! qc Cloud water mixing ratio at a given time level (kg/kg)
! qr Rainwater mixing ratio at a given time level (kg/kg)
! qi Cloud ice mixing ratio at a given time level (kg/kg)
! qs Snow mixing ratio at a given time level (kg/kg)
! qh Hail mixing ratio at a given time level (kg/kg)
!
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! ptbari Inverse base state potontial temperature (K)
! pbari Inverse base state pressure (Pascal)
! rhostr Base state density rhobar times j3 (kg/m**3)
! qvbar Base state water vapor specific humidity (kg/kg)
!
! OUTPUT:
!
! wbuoy The total buoyancy force (kg/(m*s)**2)
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature
! at a given time level (K)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure at a given time
! level (Pascal)
REAL :: qv (nx,ny,nz) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal).
REAL :: ptbari(nx,ny,nz) ! Inverse base state pot. temperature (K)
REAL :: pbari (nx,ny,nz) ! Inverse base state pressure (Pascal).
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific
! humidity(kg/kg)
REAL :: wbuoy(nx,ny,nz) ! Total buoyancy in w-eq. (kg/(m*s)**2)
REAL :: tem1 (nx,ny,nz) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
REAL :: g5
REAL :: pttem,tema
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global model control parameters
INCLUDE 'phycst.inc' ! Physical constants
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! The total buoyancy
!
! wbuoy = rhostr*g ( ptprt/ptbar-pprt/(rhobar*csndsq)+
! qvprt/(rddrv+qvbar)-(qvprt+qc+qr+qs+qi+qh)/(1+qvbar)
! -(ptprt*ptprt)/(ptbar*ptbar) !2nd-order
! +0.5*(ptprt*pprt)/(cpdcv*ptbar*pbar)) !2nd-order
!
! and rddrv=rd/rv, cp, cv, rd and rv are defined in phycst.inc.
!
! Here, the contribution from pprt (i.e., term pprt/(rhobar*csndsq))
! is evaluated inside the small time steps, therefore wbuoy
! does not include this part.
!
! The contribution from ptprt is calculated inside the small time
! steps if the potential temperature equation is solved inside
! small time steps, i.e., if ptsmlstp=1.
!
!-----------------------------------------------------------------------
!
IF( ptsmlstp == 1 ) THEN
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = 0.0
END DO
END DO
END DO
ELSE
IF (buoy2nd == 0) THEN !1st-order
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = ptprt(i,j,k)*ptbari(i,j,k)
END DO
END DO
END DO
ELSE !2nd-order
IF ( bsnesq == 1 ) THEN
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
pttem = ptprt(i,j,k)*ptbari(i,j,k)
tem1(i,j,k) = pttem-pttem*pttem
END DO
END DO
END DO
ELSE
tema = 1.0/cpdcv
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
pttem = ptprt(i,j,k)*ptbari(i,j,k)
tem1(i,j,k) = pttem* &
(1.0-pttem+0.5*pprt(i,j,k)*(tema*pbari(i,j,k)))
END DO
END DO
END DO
END IF
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Add on the contributions to the buoyancy from the water vapor
! content and the liquid and ice water loading.
!
!-----------------------------------------------------------------------
!
IF( moist == 1 .AND. ice == 0 ) THEN ! Moist case, no ice.
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = tem1(i,j,k) &
+ (qv(i,j,k) - qvbar(i,j,k))/(rddrv + qvbar(i,j,k)) &
- (qv(i,j,k) - qvbar(i,j,k) + qc(i,j,k) + qr(i,j,k)) &
/(1 + qvbar(i,j,k))
END DO
END DO
END DO
ELSE IF(moist == 1 .AND. ice == 1) THEN
! Full microphysics case, loading
! of liquid and ice water.
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = tem1(i,j,k) &
+ (qv(i,j,k) - qvbar(i,j,k))/(rddrv + qvbar(i,j,k)) &
- (qv(i,j,k) - qvbar(i,j,k) + qc(i,j,k) + qr(i,j,k) + &
qs(i,j,k) + qi(i,j,k) + qh(i,j,k))/(1 + qvbar(i,j,k))
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Then the total buoyancy:
!
! wbuoy = tem1 * rhostr * g
!
! averged to the w-point on the staggered grid.
!
!-----------------------------------------------------------------------
!
g5 = g*0.5
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
wbuoy(i,j,k)= (tem1(i,j, k )*rhostr(i,j, k ) &
+tem1(i,j,k-1)*rhostr(i,j,k-1))*g5
END DO
END DO
END DO
RETURN
END SUBROUTINE buoycy
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE PGRAD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE pgrad(nx,ny,nz, pprt, div, & 1,4
j1,j2,j3, upgrad,vpgrad,wpgrad,tem1,tem2)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the pressure gradient terms in the momentum equations.
! These terms are evaluated every small time step.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/91.
!
! MODIFICATION HISTORY:
!
! 5/05/92 (M. Xue)
! Added full documentation.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 5/25/93 (M. Xue & K. Brewster)
! Fixed and error in the vertical pressure gradient force term.
! The error is present in the finite differenced w equation
! of ARPS 3.0 User's guide. The equation in continuous form is
! correct.
!
! 9/10/94 (D. Weber & Y. Lu)
! Cleaned up documentation.
!
! 01/23/96 (Donghai Wang, Yuhe Liu and Ming Xue)
! Added the map factor.
!
! 07/31/96 (Ming Xue and Yuhe Liu)
! Added the isotropic option for divergence damping
!
! 9/14/98 (D. Weber)
! Removed operators and merged loops.
!
!-----------------------------------------------------------------------
!
! 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 direction.
!
! pprt Perturbation pressure at a given time level (Pascal)
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
!
! OUTPUT:
!
! upgrad Pressure gradient force in u-eq. (kg/(m*s)**2)
! vpgrad Pressure gradient force in v-eq. (kg/(m*s)**2)
! wpgrad Pressure gradient force in w-eq. (kg/(m*s)**2)
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
! tem2 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: pprt (nx,ny,nz) ! Perturbation pressure at a given time
! level (Pascal).
REAL :: div (nx,ny,nz) ! Mass weighted divergence
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dx.
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dy.
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian
! d(zp)/dz.
REAL :: upgrad(nx,ny,nz) ! Pressure gradient force in u-eq.
! (kg/(m*s)**2)
REAL :: vpgrad(nx,ny,nz) ! Pressure gradient force in v-eq.
! (kg/(m*s)**2)
REAL :: wpgrad(nx,ny,nz) ! Pressure gradient force in w-eq.
REAL :: tem1 (nx,ny,nz) ! Temporary work array.
REAL :: tem2 (nx,ny,nz) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: onvf
INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid parameters
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! pprt - cdvdmpv*div(i,j,k) in vertical direction
!
!-----------------------------------------------------------------------
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=pprt(i,j,k)-cdvdmpv*div(i,j,k)
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! d(pprt-cdvdmpv*div)/dz at w point - p grad. force in z-dir
!
!-----------------------------------------------------------------------
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
wpgrad(i,j,k)=dzinv*(tem1(i,j,k)-tem1(i,j,k-1))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Array div may be used only for horizontal now
!
! div = pprt - cdvdmph*div(i,j,k) in horizontal direction
!
!-----------------------------------------------------------------------
!call test_dump (div,"XXX1f_div",nx,ny,nz,1,0)
!call test_dump (pprt,"XXX1f_pprt",nx,ny,nz,1,0)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
div(i,j,k)=pprt(i,j,k)-cdvdmph*div(i,j,k)
tem1(i,j,k)=div(i,j,k)*j3(i,j,k)
END DO
END DO
END DO
!call test_dump (pprt,"XXXf_pprt",nx,ny,nz,1,0)
!call test_dump (div,"XXXf_div",nx,ny,nz,1,0)
!call test_dump (tem1,"XXXf_tem1",nx,ny,nz,1,0)
!-----------------------------------------------------------------------
!
! d(j3*pprt)/dx at u point - 1st component of p grad. force in x-dir
!
!-----------------------------------------------------------------------
DO k=1,nz-1
DO j=1,ny-1
DO i=2,nx-1
upgrad(i,j,k)=dxinv*(tem1(i,j,k)-tem1(i-1,j,k))
END DO
END DO
END DO
!call test_dump (upgrad,"XXXf_upgrad",nx,ny,nz,1,0)
!-----------------------------------------------------------------------
!
! d(j3*pprt)/dy at v point - 1st component of p grad. force in y-dir
!
!-----------------------------------------------------------------------
DO k=1,nz-1
DO j=2,ny-1
DO i=1,nx-1
vpgrad(i,j,k)=dyinv*(tem1(i,j,k)-tem1(i,j-1,k))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! If there is no terrain, i.e. the ground is flat, skip the
! following calculations.
!
!-----------------------------------------------------------------------
IF( ternopt /= 0 ) THEN
!-----------------------------------------------------------------------
!
! d(j1*pprt)/dz at u point - 2nd component of p grad. force in x-dir
! due to terrain
!
!-----------------------------------------------------------------------
DO k=2,nz-1
DO j=1,ny-1
DO i=2,nx-1
tem2(i,j,k)= j1(i,j,k)* &
0.25*((div(i-1,j,k-1)+div(i,j,k-1)) &
+(div(i-1,j,k) +div(i,j,k)))
END DO
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
DO i=2,nx-1
upgrad(i,j,k)=upgrad(i,j,k)+ &
dzinv*(tem2(i,j,k+1)-tem2(i,j,k))
END DO
END DO
END DO
!call test_dump (upgrad,"XXX1f_upgrad",nx,ny,nz,1,0)
!-----------------------------------------------------------------------
!
! d(j2*pprt)/dz at v point - 2nd component of p grad. force in y-dir
! due to terrain
!
!-----------------------------------------------------------------------
DO k=2,nz-1
DO j=2,ny-1
DO i=1,nx-1
tem2(i,j,k)= j2(i,j,k)* &
0.25*((div(i,j-1,k-1)+div(i,j,k-1)) &
+(div(i,j-1,k) +div(i,j,k)))
END DO
END DO
END DO
DO k=2,nz-2
DO j=2,ny-1
DO i=1,nx-1
vpgrad(i,j,k)=vpgrad(i,j,k)+ &
dzinv*(tem2(i,j,k+1)-tem2(i,j,k))
END DO
END DO
END DO
END IF
IF( lvldbg >= 5) THEN
CALL checkuhx(upgrad, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'upgrdx', tem2)
CALL checkuhy(upgrad, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'upgrdy', tem2)
CALL checkvhx(vpgrad, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vpgrdx', tem2)
CALL checkvhy(vpgrad, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vpgrdy', tem2)
END IF
RETURN
END SUBROUTINE pgrad
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE UVWRHO ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE uvwrho(nx,ny,nz,u,v,wcont,rhostr, & 3,1
ustr,vstr,wstr)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Compute ustr=u*rhostr, vstr=v*rhostr, wstr=wcont*rhostr.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/91.
!
! MODIFICATION HISTORY:
!
! 5/05/92 (M. Xue)
! Added full documentation.
!
! 9/10/94 (D. Weber & Y. Lu)
! Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
! 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)
! wcont Contravariant vertical velocity (m/s)
! rhostr Base state density rhobar times j3 (kg/m**3)
!
! OUTPUT:
!
! ustr u * rhostr at u-point
! vstr v * rhostr at v-point
! wstr wcont * rhostr at w-point
!
! WORK ARRAYS:
!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of 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 :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: ustr (nx,ny,nz) ! u * rhostr
REAL :: vstr (nx,ny,nz) ! v * rhostr
REAL :: wstr (nx,ny,nz) ! wcont * rhostr
INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Calculate ustr=rhostr*u, vstr=rhostr*v, wstr=rhostr*wcont
!
!-----------------------------------------------------------------------
!
CALL rhouvw(nx,ny,nz,rhostr,ustr,vstr,wstr)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx
ustr(i,j,k)=u(i,j,k)*ustr(i,j,k)
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny
DO i=1,nx-1
vstr(i,j,k)=v(i,j,k)*vstr(i,j,k)
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
wstr(i,j,k)=wcont(i,j,k)*wstr(i,j,k)
END DO
END DO
END DO
RETURN
END SUBROUTINE uvwrho
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RHOUVW ######
!###### ######
!###### Developed by ######
!###### Center for the Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rhouvw(nx,ny,nz,rhostr,rhostru,rhostrv,rhostrw) 17,3
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate rhostr averaged to u, v, and w points.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue & Hao Jin
! 3/8/1993.
!
! MODIFICATION HISTORY:
!
! 9/10/94 (D. Weber & Y. Lu)
! Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
! 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
!
! rhostr j3 times base state density rhobar(kg/m**3).
!
! OUTPUT:
!
! rhostru Average rhostr at u points (kg/m**3).
! rhostrv Average rhostr at v points (kg/m**3).
! rhostrw Average rhostr at w points (kg/m**3).
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! The number of grid points in 3
! directions
REAL :: rhostr(nx,ny,nz) ! j3 times base state density rhobar
! (kg/m**3).
REAL :: rhostru(nx,ny,nz) ! Average rhostr at u points (kg/m**3).
REAL :: rhostrv(nx,ny,nz) ! Average rhostr at v points (kg/m**3).
REAL :: rhostrw(nx,ny,nz) ! Average rhostr at w points (kg/m**3).
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
CALL avgsu(rhostr,nx,ny,nz, 1,ny-1, 1,nz-1, rhostru, rhostrw)
! rhostrw used here as a temporary array
!call test_dump (rhostr,"XXXbavgsv_rhostr",nx,ny,nz,0,1)
CALL avgsv(rhostr,nx,ny,nz, 1,nx-1, 1,nz-1, rhostrv, rhostrw)
! rhostrw used here as a temporary array
CALL avgsw(rhostr,nx,ny,nz, 1,nx-1, 1,ny-1, rhostrw)
RETURN
END SUBROUTINE rhouvw
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MAPTRM ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE maptrm( nx,ny,nz, u, v, ustr,vstr, mapfct, & 1
uforce, vforce, tem1,tem2,tem3 )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the term caused by inhomogeneous map factor in the
! advection terms of u and v.
!
! mpuadv = avgx( avgy(v*)
! * (avgy(v)*difx(mapfct_u)-avgx(u)*dify(mapfct_v)) )
!
! mpvadv = avgy( avgx(u*)
! * (avgy(v)*difx(mapfct_u)-avgx(u)*dify(mapfct_v)) )
!
! where mapfct_u is the mapfct at u points and mapfct_v at v points.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Donghai Wang, Yuhe Liu and Ming Xue
! 01/30/96
!
! MODIFICATION HISTORY:
!
! 11/06/97 (D. Weber)
! Added three additional levels to the mapfct array. The three
! levels (4,5,6) represent the inverse of the first three in order.
! The inverse map factors are computed to improve efficiency.
!
! 9/14/98 (D. Weber)
! Removed operators and merged loops and included the results into
! the forcing arrays.
!
!c#######################################################################
!
! 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
!
! mapfct Map factors at scalar, u and v points
!
! u x component of velocity (m/s)
! v y component of velocity (m/s)
! ustr u*rhostr
! vstr v*rhostr
!
! OUTPUT:
!
! uforce Acoustically inactive forcing terms in u-momentum
! equation (kg/(m*s)**2).
! vforce Acoustically inactive forcing terms in v-momentum
! equation (kg/(m*s)**2).
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: mapfct(nx,ny,8) ! Map factors at scalar, u and v points
REAL :: u (nx,ny,nz) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz) ! Total v-velocity (m/s)
REAL :: ustr (nx,ny,nz) ! u * rhostr
REAL :: vstr (nx,ny,nz) ! v * rhostr
REAL :: uforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in u-momentum equation (kg/(m*s)**2)
REAL :: vforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in v-momentum equation (kg/(m*s)**2)
REAL :: tem1 (nx,ny,nz) ! Temporary work array
REAL :: tem2 (nx,ny,nz) ! Temporary work array
REAL :: tem3 (nx,ny,nz) ! Temporary work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
INTEGER :: onvf
REAL :: tema
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid parameters
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!-----------------------------------------------------------------------
!
! Calculate avgy(v)*difx(mapfct_u) and store it to tem3
!
!-----------------------------------------------------------------------
tema = 0.5*dxinv
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem3(i,j,k) = tema*(v(i,j+1,k)+v(i,j,k))* &
(mapfct(i+1,j,2)-mapfct(i,j,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! At this point tem3 contains avgy(v)*difx(mapfct_u)
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! Calculate avgx(u)*dify(mapfct_v) and substract it from tem3
!
!-----------------------------------------------------------------------
tema = 0.5*dyinv
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem3(i,j,k) = tem3(i,j,k) - tema*(u(i+1,j,k)+u(i,j,k))* &
(mapfct(i,j+1,3)-mapfct(i,j,3))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! At this point tem3 contains
!
! avgy(v)*difx(mapfct_u) - avgx(u)*dify(mapfct_v)
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! Calculate avgsu(avgy(v*))*tem3 and subtract from uforce.
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k) = tem3(i,j,k)*0.5*(vstr(i,j+1,k)+vstr(i,j,k))
END DO
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
DO i=2,nx-1
uforce(i,j,k) = uforce(i,j,k) - 0.5*(tem2(i,j,k)+tem2(i-1,j,k))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Calculate avgsv(avgx(u*))*tem3 and add to vforce.
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k) = tem3(i,j,k)*0.5*(ustr(i+1,j,k)+ustr(i,j,k))
END DO
END DO
END DO
DO k=2,nz-2
DO j=2,ny-1
DO i=1,nx-1
vforce(i,j,k) = vforce(i,j,k)+0.5*(tem2(i,j-1,k)+tem2(i,j,k))
END DO
END DO
END DO
RETURN
END SUBROUTINE maptrm
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RHOFACTOR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rhofactor(nx,ny,nz,ptprt,pprt,qv,qc,qr,qi,qs,qh, & 1
ptbar,pbar,ptbari,pbari,qvbar,rhofct,tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the rho-factor term: rhobar/rho to correct the
! calculation of the pressure gradient force terms
!
!-----------------------------------------------------------------------
!
! AUTHOR: Donghai Wang
! 05/26/97.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! 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 direction.
!
! ptprt Perturbation potential temperature at a time level (K)
! pprt Perturbation pressure at a given time level (Pascal)
! qv Water vapor specific humidity at a given time level
! (kg/kg)
! qc Cloud water mixing ratio at a given time level (kg/kg)
! qr Rainwater mixing ratio at a given time level (kg/kg)
! qi Cloud ice mixing ratio at a given time level (kg/kg)
! qs Snow mixing ratio at a given time level (kg/kg)
! qh Hail mixing ratio at a given time level (kg/kg)
!
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! ptbari Inverse Base state pot. temperature (K)
! pbari Inverse Base state pressure (Pascal)
! qvbar Base state water vapor specific humidity (kg/kg)
!
! OUTPUT:
!
! rhofct the factor at scalar point
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
! rhofct Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature
! at a given time level (K)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure at a given time
! level (Pascal)
REAL :: qv (nx,ny,nz) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal).
REAL :: ptbari(nx,ny,nz) ! Inverse Base state pot. temperature (K)
REAL :: pbari (nx,ny,nz) ! Inverse Base state pressure (Pascal).
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific
! humidity(kg/kg)
REAL :: rhofct(nx,ny,nz) !
REAL :: tem1 (nx,ny,nz) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL :: pttem,ptem,tema
INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global model control parameters
INCLUDE 'phycst.inc' ! Physical constants
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!
! rhofct = rhobar/rho
! = 1.0/(1.0-B/g)
! = 1.0/(1.0-tem1)
!
! where B/g = ptprt/ptbar-pprt/(cpdcv*pbar)+ !1st-order
! qvprt/(rddrv+qvbar)-(qvprt+qc+qr+qs+qi+qh)/(1+qvbar(i,j,k)) !1st-order
! -(ptprt*ptprt)/(ptbar*ptbar)-1.0/(2.0*cpdcv)*(1.0/cpdcv-1.0)!2nd-order
! *(pprt*pprt)/(pbar*pbar)+(ptprt*pprt)/(2.0*cpdcv*ptbar*pbar)!2nd-order
!
! B is the buoyancy term,
! and rddrv=rd/rv,cpdcv=cp/cv, cp, cv, rd and rv are defined
! in phycst.inc.
!
!-----------------------------------------------------------------------
!
tema = 1.0/cpdcv
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
pttem = ptprt(i,j,k)*ptbari(i,j,k)
ptem = pprt(i,j,k)*tema*pbari(i,j,k)
tem1(i,j,k) = pttem-pttem*pttem &
-ptem-0.5*(1.0-cpdcv)*ptem*ptem &
+0.5*pttem*ptem
END DO
END DO
END DO
!call test_dump2(tem1,"XXX2rhofactor_tem1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(ptprt,"XXX2rhofactor_ptprt",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(ptbari,"XXX2rhofactor_ptbari",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(pprt,"XXX2rhofactor_pprt",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(pbari,"XXX2rhofactor_pbari",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!
!-----------------------------------------------------------------------
!
! Add on the contributions to the term from the water vapor
! content and the liquid and ice water loading.
!
!-----------------------------------------------------------------------
!
IF( moist == 1 .AND. ice == 0 ) THEN ! Moist case, no ice.
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = tem1(i,j,k) &
+ (qv(i,j,k) - qvbar(i,j,k))/(rddrv + qvbar(i,j,k)) &
- (qv(i,j,k) - qvbar(i,j,k) + qc(i,j,k) + qr(i,j,k)) &
/(1 + qvbar(i,j,k))
END DO
END DO
END DO
!call test_dump2(tem1,"XXX1rhofactor_tem1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
ELSE IF(moist == 1 .AND. ice == 1) THEN
! Full microphysics case, loading
! of liquid and ice water.
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = tem1(i,j,k) &
+ (qv(i,j,k) - qvbar(i,j,k))/(rddrv + qvbar(i,j,k)) &
- (qv(i,j,k) - qvbar(i,j,k) + qc(i,j,k) + qr(i,j,k) + &
qs(i,j,k) + qi(i,j,k) + qh(i,j,k))/(1 + qvbar(i,j,k))
END DO
END DO
END DO
!call test_dump2(tem1,"XXXrhofactor_tem1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qv,"XXXrhofactor_qv",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qvbar,"XXXrhofactor_qvbar",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qc,"XXXrhofactor_qc",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr,"XXXrhofactor_qr",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qs,"XXXrhofactor_qs",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qi,"XXXrhofactor_qi",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qh,"XXXrhofactor_qh",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
END IF
!
!-----------------------------------------------------------------------
!
!
! rhofct = 1.0/(1.0-B/g), B is the buoyancy, tem1=B/g
!
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
rhofct(i,j,k)=1.0/(1.0-tem1(i,j,k))
END DO
END DO
END DO
RETURN
END SUBROUTINE rhofactor