!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SOLVUV ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE solvuv(mptr, nx,ny,nz, exbcbufsz, dtsml1, curtsml, & 1,43
u,v,wcont,pprt, &
udteb,udtwb,udtnb,udtsb, &
vdteb,vdtwb,vdtnb,vdtsb, &
rhostr, uforce,vforce,ubar,vbar, &
x,y,z,zp, mapfct, j1,j2,j3,j3inv, &
exbcbuf,rhofct, &
rhostru,rhostrv,rhostrw,dtmrxrsu,dtmryrsv,wpgrad, &
upgrad,vpgrad,tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinate the time stepping of u and v momentum equations and the
! setting of boundary conditions.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Ming Xue
! 10/10/92
!
! MODIFICATION HISTORY:
!
! 5/20/92 (K. Droegemeier and M. Xue)
! Added full documentation.
!
! 6/07/92 ( M. Xue)
! Now only tfuture fields of u, v, w and pprt are passed in.
!
! 2/10/93 (K. Droegemeier)
! Cleaned up documentation.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 9/1/94 (Y. Lu)
! Cleaned up documentation.
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 10/16/97 (Donghai Wang)
! Using total density (rho) for the calculation of the pressure
! gradient force terms.
!
! 11/06/97 (Dan 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/28/98 (Dan Weber)
! Reformulated do-loops to improve efficiency, brought in
! pre-computed variables. Note u,vforce contains dtsml and
! 1/rhostru,v.
!
! 07/23/2001 (K. Brewster)
! Added mptr to argument list.
! Added optional writing of total divergence as a diagnostic noise
! parameter.
!
!-----------------------------------------------------------------------
!
! 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
! dtsml1 Local value of small time step size
! curtsml Local curttim at a small time step
!
! u x component of velocity at tfuture (m/s)
! v y component of velocity at tfuture (m/s)
! wcont Contravariant vertical velocity (m/s)
!
! pprt Perturbation pressure at tfuture (Pascal)
!
! udteb Time tendency of the u field at the east boundary
! udtwb Time tendency of the u field at the west boundary
! udtnb Time tendency of the u field at the north boundary
! udtsb Time tendency of the u field at the south boundary
!
! vdteb Time tendency of the v field at the east boundary
! vdtwb Time tendency of the v field at the west boundary
! vdtnb Time tendency of the v field at the north boundary
! vdtsb Time tendency of the v field at the south boundary
!
! rhostr Base state density rhobar times j3 (kg/m**3)
!
! uforce Acoustically inactive forcing terms in u-eq.
! (kg/(m*s)**2)
! vforce Acoustically inactive forcing terms in v-eq.
! (kg/(m*s)**2)
!
! ubar Base state x velocity component (m/s)
! vbar Base state y velocity component (m/s)
!
! 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
! j3inv Inverse of j3
!
! 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).
!
! OUTPUT:
!
! u x component of velocity at time tfuture updated
! in time by dtsml1 (m/s)
! v y component of velocity at time tfuture updated
! in time by dtsml1 (m/s)
! w Vertical component of Cartesian velocity at tfuture
! updated in time by dtsml1 (m/s)
! wpgrad Pressure gradient term in w-eq. (kg/(m*s)**2)
!
! WORK ARRAYS:
!
! upgrad Pressure gradient term in u-eq. (kg/(m*s)**2)
! vpgrad Pressure gradient term in v-eq. (kg/(m*s)**2)
! tem1 Temporary work array
! tem2 Temporary work array
! tem3 Temporary work array
! dtmrxrsu dtsml*mapfct*avgx(rhofct)/rhostru
! dtmryrsv dtsml*mapfct*avgy(rhofct)/rhostrv
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INTEGER :: mptr
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: u (nx,ny,nz) ! u-velocity at tfuture (m/s)
REAL :: v (nx,ny,nz) ! v-velocity at tfuture (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure at tfuture
! (Pascal)
REAL :: udteb (ny,nz) ! Time tendency of u at east boundary
REAL :: udtwb (ny,nz) ! Time tendency of u at west boundary
REAL :: udtnb (nx,nz) ! Time tendency of u at north boundary
REAL :: udtsb (nx,nz) ! Time tendency of u at south boundary
REAL :: vdteb (ny,nz) ! Time tendency of v at east boundary
REAL :: vdtwb (ny,nz) ! Time tendency of v at west boundary
REAL :: vdtnb (nx,nz) ! Time tendency of v at north boundary
REAL :: vdtsb (nx,nz) ! Time tendency of v at south boundary
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
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 :: ubar(nx,ny,nz)
REAL :: vbar(nx,ny,nz)
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 compui-
! 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, u and v 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 :: j3inv (nx,ny,nz) ! Inverse of j3
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: rhofct(nx,ny,nz) ! rho-factor: rhobar/rho
REAL :: rhostru(nx,ny,nz) ! Averaged rhostr at u points (kg/m**3).
REAL :: rhostrv(nx,ny,nz) ! Averaged rhostr at v points (kg/m**3).
REAL :: rhostrw(nx,ny,nz) ! Averaged rhostr at w points (kg/m**3).
REAL :: wpgrad(nx,ny,nz) ! Pressure gradient term in w-eq.
!
!-----------------------------------------------------------------------
!
! Temporary WORK ARRAYS:
!
!-----------------------------------------------------------------------
!
REAL :: upgrad(nx,ny,nz) ! Pressure gradient term in u-eq.
REAL :: vpgrad(nx,ny,nz) ! Pressure gradient term in v-eq.
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 :: dtmrxrsu(nx,ny,nz) ! dtsml*map*avgx(rhofct)/rhostru
REAL :: dtmryrsv(nx,ny,nz) ! dtsml*map*avgy(rhofct)/rhostrv
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL :: curtsml ! Local curttim at a small time step
REAL :: dtsml1 ! Local small time step
REAL :: sumdiv
INTEGER :: i, j, k
INTEGER :: istart,iend,jstart,jend
INTEGER :: ebc1,wbc1,nbc1,sbc1
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
INTEGER :: noisewrt
INTEGER :: ncalls(mgrdmax), nchdiv1(mgrdmax)
CHARACTER (LEN=128 ) :: divfn
INTEGER :: ldivfn,istat
INTEGER :: mptag
SAVE ncalls,nchdiv1
DATA ncalls /mgrdmax*0/
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
noisewrt=0
ncalls(mptr)=ncalls(mptr) + 1
IF(noisewrt == 1 .AND. ncalls(mptr) == 1 ) THEN
divfn = runname(1:lfnkey)//'.absdiv'
ldivfn = 7 + lfnkey
IF(nestgrd == 1) THEN
WRITE(divfn((ldivfn+1):(ldivfn+4)),'(a,i2.2)')'.g',mptr
ldivfn = ldivfn + 4
END IF
WRITE(6,'(1x,a,a,a/,1x,a)') &
'Check to see if file ',divfn(1:ldivfn),' already exists.', &
'If so, append a version number to the filename.'
CALL fnversn
( divfn, ldivfn )
CALL getunit ( nchdiv1(mptr) )
OPEN(nchdiv1(mptr),FORM='formatted',STATUS='new', &
FILE=divfn(1:ldivfn),IOSTAT=istat)
IF( istat /= 0) THEN
WRITE(6,'(/a,i2,/a/)') &
' Error occured when opening file '//divfn(1:ldivfn)// &
' using FORTRAN unit ',nchdiv1(mptr), &
' Program stopped in MAXMIN.'
CALL arpsstop(' ',1)
END IF
WRITE(nchdiv1(mptr),'(a)') runname
WRITE(nchdiv1(mptr),'(t2,a,t15,a)') 'Time_(s)','MnAbsDiv'
END IF
!-----------------------------------------------------------------------
!
! Divergence damping is activated if divdmp = 1.
! Otherwise, this portion of the code is skipped to save
! computations.
!
!-----------------------------------------------------------------------
IF( divdmp == 1 .OR. divdmp == 2 ) THEN
!
!-----------------------------------------------------------------------
!
! The divergence damping terms are defined as
!
! d(cdvdmph*div)/dx for u eq.,
! d(cdvdmph*div)/dy for v eq.,
! d(cdvdmpv*div)/dz for w eq.,
!
! where
!
! div = (difx(u*rhostr) + dify(v*rhostr) + difz(wcont*rhostr))/j3
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! Compute the momentum divergence term, defined as
!
! div = d(u*rhostr)/dx + d(v*rhostr)/dy + d(wcont*rhostr)/dz.
!
! and combine the pressure and divergence damping into one
! array so that the pressure gradient force and divergence damping
! can be calculated in a single step.
!
!-----------------------------------------------------------------------
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx
upgrad(i,j,k)=u(i,j,k)*rhostru(i,j,k)*mapfct(i,j,5)
END DO
END DO
END DO
!call test_dump (upgrad,"XXX0upgrad",nx,ny,nz,1,0)
DO k=1,nz-1
DO j=1,ny
DO i=1,nx-1
vpgrad(i,j,k)=v(i,j,k)*rhostrv(i,j,k)*mapfct(i,j,6)
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=wcont(i,j,k)*rhostrw(i,j,k)
END DO
END DO
END DO
IF(noisewrt == 1 .AND. MOD(ncalls(mptr),2) == 0) THEN
IF(lbcopt == 2) THEN
istart=MAX(2,(ngbrz+1))
iend=MIN((nx-2),(nx-ngbrz-1))
jstart=MAX(2,(ngbrz+1))
jend=MIN((ny-2),(ny-ngbrz-1))
ELSE
istart=2
iend=nx-2
jstart=2
jend=ny-2
END IF
sumdiv=0.
DO k=2,nz-2
DO j=jstart,jend
DO i=istart,iend
sumdiv = sumdiv+abs(j3inv(i,j,k) &
* ( mapfct(i,j,7) &
* ((upgrad(i+1,j,k)-upgrad(i,j,k))*dxinv &
+(vpgrad(i,j+1,k)-vpgrad(i,j,k))*dyinv) &
+ (tem1(i,j,k+1)-tem1(i,j,k))*dzinv ))
END DO
END DO
END DO
sumdiv=1.0E06*sumdiv/ &
float((iend-istart+1)*(jend-jstart+1)*(nz-3))
write(nchdiv1(mptr),'(f10.2,f12.4)') (curtim+2*dtsml1),sumdiv
END IF
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem3(i,j,k) = j3inv(i,j,k) &
* ( mapfct(i,j,7) &
* ((upgrad(i+1,j,k)-upgrad(i,j,k))*dxinv &
+(vpgrad(i,j+1,k)-vpgrad(i,j,k))*dyinv) &
+ (tem1(i,j,k+1)-tem1(i,j,k))*dzinv )
END DO
END DO
END DO
!call test_dump (upgrad,"XXXtem3_upgrad",nx,ny,nz,1,0)
!call test_dump (vpgrad,"XXXtem3_vpgrad",nx,ny,nz,2,0)
!call test_dump (vpgrad,"XXXtem3_vpgrad2",nx,ny,nz,1,0)
!call test_dump (tem1,"XXXtem3_tem1",nx,ny,nz,3,0)
!call test_dump (tem1,"XXXtem3_tem12",nx,ny,nz,1,0)
ELSE
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem3(i,j,k)= 0.0
END DO
END DO
END DO
END IF
!-----------------------------------------------------------------------
!
! Compute the pressure gradient force terms for the three
! momentum equations and store the values in upgrad, vpgrad,
! and wpgrad.
!
! When divergence damping is activated (divdmp=1), these arrays also
! contain the divergence damping terms.
!
!-----------------------------------------------------------------------
!call test_dump (tem3,"XXXBpgrad_tem3",nx,ny,nz,1,0)
CALL pgrad(nx,ny,nz, pprt, tem3, &
j1,j2,j3, upgrad,vpgrad,wpgrad,tem1,tem2)
!-----------------------------------------------------------------------
!
! Integrate the u, and v equations forward by one small
! timestep dtsml1 using a forward-in-time integration scheme
! (that is, forward relative to the pressure gradient terms).
!
!-----------------------------------------------------------------------
!call test_dump (dtmrxrsu,"XXX3dtmrxrsu",nx,ny,nz,1,0)
IF( nestgrd == 1 .AND. mgrid /= 1 ) THEN
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-1
u(i,j,k)=u(i,j,k) + uforce(i,j,k) - &
dtmrxrsu(i,j,k)*upgrad(i,j,k)
END DO
END DO
END DO
DO k=2,nz-2
DO j=2,ny-1
DO i=2,nx-2
v(i,j,k)=v(i,j,k) + vforce(i,j,k) - &
dtmryrsv(i,j,k)*vpgrad(i,j,k)
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Set the lateral boundary conditions for u and v given the
! time tendencies in the case of nested inner grid.
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=1,ny-1
u( 1,j,k)=u( 1,j,k)+dtsml1 * udtwb(j,k)
u(nx,j,k)=u(nx,j,k)+dtsml1 * udteb(j,k)
END DO
END DO
DO k=2,nz-2
DO i=2,nx-1
u(i, 1,k)=u(i, 1,k)+dtsml1 * udtsb(i,k)
u(i,ny-1,k)=u(i,ny-1,k)+dtsml1 * udtnb(i,k)
END DO
END DO
DO k=2,nz-2
DO j=1,ny
v( 1,j,k)=v( 1,j,k)+dtsml1 * vdtwb(j,k)
v(nx-1,j,k)=v(nx-1,j,k)+dtsml1 * vdteb(j,k)
END DO
END DO
DO k=2,nz-2
DO i=2,nx-2
v(i, 1,k)=v(i, 1,k)+dtsml1 * vdtsb(i,k)
v(i,ny,k)=v(i,ny,k)+dtsml1 * vdtnb(i,k)
END DO
END DO
!-----------------------------------------------------------------------
!
! Set the top and bottom boundary conditions for u and v.
!
!-----------------------------------------------------------------------
!call test_dump (u,"XXXsolve_u",nx,ny,nz,1,0)
! bc's for mp not implemented for nested grids
CALL acct_interrupt(bc_acct)
CALL bcu(nx,ny,nz, dtsml1, u, udteb,udtwb,udtnb,udtsb, &
0,0,0,0 ,tbc,bbc)
CALL acct_stop_inter
!call test_dump (u,"XXXAsolve_u",nx,ny,nz,1,1)
!call test_dump (v,"XXXsolve_v",nx,ny,nz,2,0)
! bc's for mp not implemented for nested grids
CALL acct_interrupt(bc_acct)
CALL bcv(nx,ny,nz, dtsml1, v, vdteb,vdtwb,vdtnb,vdtsb, &
0,0,0,0 ,tbc,bbc)
CALL acct_stop_inter
!call test_dump (v,"XXXAsolve_v",nx,ny,nz,2,1)
ELSE
!call test_dump (u,"XXXuBuforce",nx,ny,nz,1,0)
!call test_dump (uforce,"XXXuforce",nx,ny,nz,1,0)
!call test_dump (dtmrxrsu,"XXXdtmrxrsu",nx,ny,nz,1,0)
!call test_dump (upgrad,"XXXupgrad",nx,ny,nz,1,0)
DO k=2,nz-2
DO j=1,ny-1
DO i=2,nx-1
u(i,j,k)=u(i,j,k) + uforce(i,j,k) - &
dtmrxrsu(i,j,k)*upgrad(i,j,k)
END DO
END DO
END DO
!call test_dump (dtmrxrsu,"XXXdtmrxrsu_tst",nx,ny,nz,1,0)
DO k=2,nz-2
DO j=2,ny-1
DO i=1,nx-1
v(i,j,k)=v(i,j,k) + vforce(i,j,k) - &
dtmryrsv(i,j,k)*vpgrad(i,j,k)
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Set the boundary conditions for u and v.
!
!-----------------------------------------------------------------------
IF ( lbcopt == 1 ) THEN
ebc1=ebc
wbc1=wbc
sbc1=sbc
nbc1=nbc
IF( ebc == 4 ) ebc1=0
IF( wbc == 4 ) wbc1=0
IF( nbc == 4 ) nbc1=0
IF( sbc == 4 ) sbc1=0
!call test_dump (u,"XXX1solve_u",nx,ny,nz,1,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(u,nx,ny,nz,ebc,wbc,1,mptag,tem1)
CALL mprecv2dew(u,nx,ny,nz,ebc,wbc,1,mptag,tem1)
CALL mpsend2dns(u,nx,ny,nz,nbc1,sbc1,1,mptag,tem1)
CALL mprecv2dns(u,nx,ny,nz,nbc1,sbc1,1,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
CALL bcu(nx,ny,nz, dtsml1, u, udteb,udtwb,udtnb,udtsb, &
ebc,wbc,nbc1,sbc1,tbc,bbc)
CALL acct_stop_inter
!call test_dump (u,"XXX1Asolve_u",nx,ny,nz,1,1)
!call test_dump (v,"XXX1solve_v",nx,ny,nz,2,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(v,nx,ny,nz,ebc1,wbc1,2,mptag,tem1)
CALL mprecv2dew(v,nx,ny,nz,ebc1,wbc1,2,mptag,tem1)
CALL mpsend2dns(v,nx,ny,nz,nbc,sbc,2,mptag,tem1)
CALL mprecv2dns(v,nx,ny,nz,nbc,sbc,2,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
CALL bcv(nx,ny,nz, dtsml1, v, vdteb,vdtwb,vdtnb,vdtsb, &
ebc1,wbc1,nbc,sbc,tbc,bbc)
CALL acct_stop_inter
!call test_dump (v,"XXX1Asolve_v",nx,ny,nz,2,1)
ELSE
!call test_dump (u,"XXX2solve_u",nx,ny,nz,1,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(u,nx,ny,nz,0,0,1,mptag,tem1)
CALL mprecv2dew(u,nx,ny,nz,0,0,1,mptag,tem1)
CALL mpsend2dns(u,nx,ny,nz,0,0,1,mptag,tem1)
CALL mprecv2dns(u,nx,ny,nz,0,0,1,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
CALL bcu(nx,ny,nz, dtsml1, u, udteb,udtwb,udtnb,udtsb, &
0,0,0,0,tbc,bbc)
CALL acct_stop_inter
!call test_dump (u,"XXX2Asolve_u",nx,ny,nz,1,1)
!call test_dump (v,"XXX2solve_v",nx,ny,nz,2,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(v,nx,ny,nz,0,0,2,mptag,tem1)
CALL mprecv2dew(v,nx,ny,nz,0,0,2,mptag,tem1)
CALL mpsend2dns(v,nx,ny,nz,0,0,2,mptag,tem1)
CALL mprecv2dns(v,nx,ny,nz,0,0,2,mptag,tem1)
END IF
CALL acct_interrupt(mp_acct)
CALL bcv(nx,ny,nz, dtsml1, v, vdteb,vdtwb,vdtnb,vdtsb, &
0,0,0,0,tbc,bbc)
!call test_dump (v,"XXX2Asolve_v",nx,ny,nz,2,1)
CALL exbcuv(nx,ny,nz, curtsml, u,v, &
exbcbuf(nu0exb), exbcbuf(nv0exb), &
exbcbuf(nudtexb),exbcbuf(nvdtexb))
CALL acct_stop_inter
!call test_dump (u,"XXX2AAsolve_u",nx,ny,nz,1,1)
!call test_dump (v,"XXX2AAsolve_v",nx,ny,nz,2,1)
END IF
END IF
RETURN
END SUBROUTINE solvuv
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SOLVWPEX ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE solvwpex(mptr, nx,ny,nz,exbcbufsz, dtsml1, curtsml, & 1,51
u,v,w,wcont,ptprt,pprt,phydro, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
rhostr,ptbar,ptbari,pbari,csndsq, &
wforce,wpgrad,pforce, &
x,y,z,zp, mapfct, j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
rhostru,rhostrv,rhostrw, &
exbcbuf,rhofct, &
div,pdiv,tem1,tem2,tem3,rstpbi,rstptbi, &
dtrzrstw,dxj3xm,dyj3ym,rstj3i)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinate the time stepping of the w momentum and pressure equations
! using explicit time integration scheme. Divergence damping and the
! acoustically inactive forcing terms in w and pressure equations
! have been evaluated prior to this subroutine and are passed into
! this routine.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/91
!
! MODIFICATION HISTORY:
!
! 5/20/92 (K. Droegemeier and M. Xue)
! Added full documentation.
!
! 6/07/92 ( M. Xue)
! Now only tfuture fields of pprt, u, v and w are passed in.
!
! 2/10/93 (K. Droegemeier)
! Cleaned up documentation.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 9/6/94 (M.Xue)
! Pressure perturbation buoyancy and base state pressure
! advection terms moved into the small time steps.
!
! 9/1/94 (Y. Lu)
! Cleaned up documentation.
!
! 01/28/95 (G. Bassett)
! Option added to turn off buoyancy terms (for buoyopt=0).
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 4/27/1996 (M. Xue)
! Added code for peqopt=2 case.
!
! 5/6/96 (M. Xue)
! Replaced csndsq in pressure buoyancy term by cpdcv*pbar/rhobar.
!
! 10/16/97 (Donghai Wang)
! Using total density (rho) for the calculation of the pressure
! gradient force terms.
!
! 10/16/97 (Donghai Wang)
! Added the second order terms in the linerized buoyancy terms.
!
! 11/05/97 (D. Weber)
! Added phydro array for use in the bottom boundary condition for
! perturbation 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.
!
! 9/18/98 (D. Weber)
! Added arrays aj3x,y,z,pbari,ptbari,rstpbi,rstptbi,dtrzrstw,dxj3xm,
! dyj3ym,rstj3i.
! do not alter arrays!!!
!
! 07/23/2001 (K. Brewster)
! Added mptr to argument list.
! Added optional writing of mean abs(dp/dt) as a diagnostic noise
! parameter.
!
!-----------------------------------------------------------------------
!
! 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
!
! dtsml1 Local value of small time step size
! curtsml Local curttim at a small time step
!
! u x component of velocity at tfuture (m/s)
! v y component of velocity at tfuture (m/s)
! w Vertical component of Cartesian velocity at tfuture
! (m/s)
! wcont Contravariant vertical velocity (m/s)
! ptprt Perturbation potential temperature at all time levels
! (K)
! pprt Perturbation pressure at tfuture (Pascal)
!
! phydro Big time step forcing term for use in computing the
! hydrostatic pressure at k=1.
!
! wdteb Time tendency of the w field at the east boundary
! wdtwb Time tendency of the w field at the west boundary
! wdtnb Time tendency of the w field at the north boundary
! wdtsb Time tendency of the w field at the south boundary
!
! pdteb Time tendency of the pressure field at the east
! boundary
! pdtwb Time tendency of the pressure field at the west
! boundary
! pdtnb Time tendency of the pressure field at the north
! boundary
! pdtsb Time tendency of the pressure field at the south
! boundary
!
! rhostr Base state density rhobar times j3 (kg/m**3)
! ptbar Base state potential temperature (K)
! ptbari Inverse base state potential temperature (K)
! pbari Inverse base state pressure (Pascal)
! csndsq Sound wave speed squared.
!
! wforce Acoustically inactive forcing terms in w-eq.
! (kg/(m*s)**2)
! pforce Acoustically inactive terms in pressure eq. (Pascal/s)
!
! 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
!
! 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).
!
! OUTPUT:
!
! pprt Perturbation pressure at tfuture updated in time
! by stsml1 (Pascal)
! w Vertical component of Cartesian velocity at tfuture
! updated in time by dtsml1 (m/s)
!
! WORK ARRAYS:
!
! wpgrad Pressure gradient term in w-eq. (kg/(m*s)**2)
! div Velocity divergence (1/s), a local array
! pdiv Divergence term in pressure eq.
! (Pascal/s), a local array
! tem1 Temporary work array
! tem2 Temporary work array
! tem3 Temporary work array
! rstpbi rhostr*pbari
! rstptbi rhostr*ptbari
! dtrzrstw dtsml*avgz(rhofct)/rhostrw
! dxj3xm dxinv*aj3x*mapfct(i,j,5)
! dyj3ym dyinv*aj3y*mapfct(i,j,6)
! rstj3i rhostr*j3inv
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INTEGER :: mptr
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: curtsml ! Local curttim at a small time step
REAL :: dtsml1 ! Local small time step size (s)
REAL :: u (nx,ny,nz) ! u-velocity at tfuture (m/s)
REAL :: v (nx,ny,nz) ! v-velocity at tfuture (m/s)
REAL :: w (nx,ny,nz) ! w-velocity at tfuture (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure at tfuture
! (Pascal)
REAL :: phydro(nx,ny) ! Big time step forcing for computing
! hydrostatic pprt at k=1.
REAL :: wdteb (ny,nz) ! Time tendency of w at east boundary
REAL :: wdtwb (ny,nz) ! Time tendency of w at west boundary
REAL :: wdtnb (nx,nz) ! Time tendency of w at north boundary
REAL :: wdtsb (nx,nz) ! Time tendency of w at south boundary
REAL :: pdteb (ny,nz) ! Time tendency of pressure at east
! boundary
REAL :: pdtwb (ny,nz) ! Time tendency of pressure at west
! boundary
REAL :: pdtnb (nx,nz) ! Time tendency of pressure at north
! boundary
REAL :: pdtsb (nx,nz) ! Time tendency of pressure at south
! boundary
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: ptbari(nx,ny,nz) ! Inverse base state pot. temperature (K)
REAL :: pbari (nx,ny,nz) ! Inverse base state pressure (Pascal).
REAL :: csndsq(nx,ny,nz) ! Sound wave speed squared.
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: wforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in w-momentum equation (kg/(m*s)**2)
REAL :: pforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in the pressure equation (Pascal/s)
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, u and v 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 :: 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
! defined as d( zp )/d( z ).
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: rhofct(nx,ny,nz) ! rho-factor:rhobar/rho
REAL :: rhostru(nx,ny,nz) ! Averaged rhostr at u points (kg/m**3).
REAL :: rhostrv(nx,ny,nz) ! Averaged rhostr at v points (kg/m**3).
REAL :: rhostrw(nx,ny,nz) ! Averaged rhostr at w points (kg/m**3).
!
!-----------------------------------------------------------------------
!
! Temporary work arrays:
!
!-----------------------------------------------------------------------
!
REAL :: wpgrad(nx,ny,nz) ! Pressure gradient term in w-eq.
REAL :: div (nx,ny,nz) ! Velocity divergence (1/s)
REAL :: pdiv (nx,ny,nz) ! Divergence term in pressure eq.
! (Pascal/s)
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 :: rstpbi(nx,ny,nz) ! rhostr*pbari
REAL :: rstptbi(nx,ny,nz) ! rhostr*ptbari
REAL :: dtrzrstw(nx,ny,nz) ! dtsml1*avgz(rhofct)/rhostrw
REAL :: dxj3xm(nx,ny,nz) ! dxinv*aj3x*mapfct(i,j,5)
REAL :: dyj3ym(nx,ny,nz) ! dyinv*aj3y*mapfct(i,j,6)
REAL :: rstj3i(nx,ny,nz) ! rhostr*j3inv
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'phycst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
INTEGER :: istart,iend,jstart,jend
INTEGER :: ebc1,wbc1,nbc1,sbc1
INTEGER :: wpprt ! Switch for pprt/(rhobar*csndsq) term in w-eq.
INTEGER :: prgw ! Switch for rhobar*g*w term in w-eq.
INTEGER :: noisewrt
REAL :: prgwg5, g05
REAL :: pttem,ptem,tem,tema,temb,temc
REAL :: sumdpdt
CHARACTER (LEN=128 ) :: dpdtfn
INTEGER :: ldpdtfn,istat
INTEGER :: ncalls(mgrdmax), nchdpdt1(mgrdmax)
SAVE ncalls,nchdpdt1
DATA ncalls /mgrdmax*0/
INTEGER :: mptag
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
noisewrt=0
ncalls(mptr)=ncalls(mptr) + 1
IF (noisewrt == 1 .AND. ncalls(mptr) == 1) THEN
dpdtfn = runname(1:lfnkey)//'.dpdt'
ldpdtfn = 5 + lfnkey
IF(nestgrd == 1) THEN
WRITE(dpdtfn((ldpdtfn+1):(ldpdtfn+4)),'(a,i2.2)')'.g',mptr
ldpdtfn = ldpdtfn + 4
END IF
WRITE(6,'(1x,a,a,a/,1x,a)') &
'Check to see if file ',dpdtfn(1:ldpdtfn),' already exists.', &
'If so, append a version number to the filename.'
CALL fnversn( dpdtfn, ldpdtfn )
CALL getunit ( nchdpdt1(mptr) )
OPEN(nchdpdt1(mptr),FORM='formatted',STATUS='new', &
FILE=dpdtfn(1:ldpdtfn),IOSTAT=istat)
IF( istat /= 0) THEN
WRITE(6,'(/a,i2,/a/)') &
' Error occured when opening file '//dpdtfn(1:ldpdtfn)// &
' using FORTRAN unit ',nchdpdt1(mptr), &
' Program stopped in SOLVWPEX.'
CALL arpsstop('arpsstop called from SOLVWPEX error opening '// &
'file',1)
END IF
WRITE(nchdpdt1(mptr),'(a)') runname
WRITE(nchdpdt1(mptr),'(t2,a,t15,a)') 'Time (s)','MnAbsdpdt'
END IF
g05 = g*0.5
IF( bsnesq == 1 ) THEN
wpprt = 0 ! Switch for pprt/(rhobar*csndsq) term in w-eq.
prgw = 0 ! Switch for rhobar*g*w term in w-eq.
ELSE
wpprt = 1 ! Switch for pprt/(rhobar*csndsq) term in w-eq.
prgw = 1 ! Switch for rhobar*g*w term in w-eq.
END IF
!
!-----------------------------------------------------------------------
!
! The contribution from ptprt to the buoyancy term is calculated
! inside the small time steps if the potential temperature equation
! is solved inside small time steps, i.e., if ptsmlstp=1.
!
! The contribution from pressure perturbation to the buoyancy is
! always calculated inside the small time steps for better
! computational stability
!
! If buoyopt = 0 then turn off all buoyancy.
!
!-----------------------------------------------------------------------
!
tem = 0.5*(1.0-cpdcv)
tema = 1.0/cpdcv
IF ( buoyopt /= 0 ) THEN
CALL acct_interrupt(buoy_acct)
IF ( ptsmlstp == 1 ) THEN
IF ( buoy2nd == 0) THEN !1st-order
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)=rstptbi(i,j,k)*ptprt(i,j,k) &
-wpprt*pprt(i,j,k)*tema*rstpbi(i,j,k)
END DO
END DO
END DO
ELSE !2nd-order
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 = wpprt*pprt(i,j,k)*(tema*pbari(i,j,k))
tem2(i,j,k)=rhostr(i,j,k)*( &
pttem-pttem*pttem-ptem-tem*ptem*ptem &
+ 0.5*pttem*ptem)
END DO
END DO
END DO
END IF
ELSE
IF (buoy2nd == 0) THEN !1st-order
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)= &
-wpprt*pprt(i,j,k)*rstpbi(i,j,k)*tema
END DO
END DO
END DO
ELSE !2nd-order
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
ptem = wpprt*pprt(i,j,k)*(tema*pbari(i,j,k))
tem2(i,j,k)=-rhostr(i,j,k)*(ptem+tem*ptem*ptem )
END DO
END DO
END DO
END IF
END IF
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=(tem2(i,j,k)+tem2(i,j,k-1))*g05
END DO
END DO
END DO
CALL acct_stop_inter
ELSE
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
END IF
!
!-----------------------------------------------------------------------
!
! To assume mirror symmetry about the top and bottom boundary,
! the density/temperature has zero gradient across the boundary
! but g is assumed to change sign, so that the buoyancy is zero
! at the boundary. c.f. subroutine BUOYCY.
!
!-----------------------------------------------------------------------
!
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,2)=0.0
tem1(i,j,nz-1)=0.0
END DO
END DO
!-----------------------------------------------------------------------
!
! Integrate w equations forward by one small timestep dtsml1 using
! a forward-in-time integration scheme (that is, forward relative
! to the pressure gradient terms).
!
!-----------------------------------------------------------------------
IF( nestgrd == 1 .AND. mgrid /= 1 ) THEN
DO k=2,nz-1
DO j=2,ny-2
DO i=2,nx-2
w(i,j,k)=w(i,j,k) + wforce(i,j,k) - &
(wpgrad(i,j,k)-tem1(i,j,k))*dtrzrstw(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set the lateral boundary conditions for w given the time tendencies
! in the case of nested inner grid.
!
!-----------------------------------------------------------------------
!
DO k=3,nz-2
DO j=1,ny-1
w( 1,j,k)=w( 1,j,k)+dtsml1 * wdtwb(j,k)
w(nx-1,j,k)=w(nx-1,j,k)+dtsml1 * wdteb(j,k)
END DO
END DO
DO k=3,nz-2
DO i=2,nx-2
w(i, 1,k)=w(i, 1,k)+dtsml1 * wdtsb(i,k)
w(i,ny-1,k)=w(i,ny-1,k)+dtsml1 * wdtnb(i,k)
END DO
END DO
ELSE
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
w(i,j,k)=w(i,j,k) + wforce(i,j,k) - &
(wpgrad(i,j,k)-tem1(i,j,k))*dtrzrstw(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Apply the lateral boundary conditions for w.
!
! For the open boundary case, w at the lateral boundarues are solved
! from w equation directly therefore does not need to be reset.
!
!-----------------------------------------------------------------------
!
!call test_dump (w,"XXXsolve_w",nx,ny,nz,3,0)
IF( lbcopt == 1 ) THEN ! Internal boundary conditions
ebc1=ebc
wbc1=wbc
sbc1=sbc
nbc1=nbc
IF( ebc == 4 ) ebc1=0
IF( wbc == 4 ) wbc1=0
IF( nbc == 4 ) nbc1=0
IF( sbc == 4 ) sbc1=0
!call test_dump (w,"XXX4solve_w",nx,ny,nz,3,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mprecv2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mpsend2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
CALL mprecv2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
END IF
CALL acct_interrupt(mp_acct)
CALL lbcw(nx,ny,nz,dtsml,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, &
ebc1,wbc1,nbc1,sbc1)
CALL acct_stop_inter
!call test_dump (w,"XXX4Asolve_w",nx,ny,nz,3,1)
ELSE ! Apply zero gradient condition
ebc1 = 3
wbc1 = 3
sbc1 = 3
nbc1 = 3
IF( ebc == 0 ) ebc1 = 0
IF( wbc == 0 ) wbc1 = 0
IF( nbc == 0 ) nbc1 = 0
IF( sbc == 0 ) sbc1 = 0
!call test_dump (w,"XXX5solve_w",nx,ny,nz,3,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mprecv2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mpsend2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
CALL mprecv2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
END IF
CALL acct_interrupt(mp_acct)
CALL lbcw(nx,ny,nz,dtsml,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, &
ebc1,wbc1,nbc1,sbc1)
CALL acct_stop_inter
!call test_dump (w,"XXX5Asolve_w",nx,ny,nz,3,1)
END IF
END IF
!call test_dump (w,"XXXAsolve_w",nx,ny,nz,3,1)
!
!-----------------------------------------------------------------------
!
! Calculate wcont at time tfuture, including the boundary
! points. Wcont at the lateral boundaries is calculated
! from boundary u, v and w. Wcont at the top and bottom
! depends on the boundary condition option.
!
!-----------------------------------------------------------------------
!
!call test_dump (w,"XXX1sBwc_w1",nx,ny,nz,1,0)
CALL wcontra(nx,ny,nz,u,v,w,mapfct,j1,j2,j3,aj3z, &
rhostr,rhostru,rhostrv,rhostrw,wcont,tem1,tem2)
!
!-----------------------------------------------------------------------
!
! Set the top and bottom boundary conditions for w based on u, v and
! wcont at the top and bottom boundaries.
!
!-----------------------------------------------------------------------
!
CALL acct_interrupt(bc_acct)
CALL vbcw(nx,ny,nz,w,wcont,tbc,bbc,u,v, &
rhostr,rhostru,rhostrv,rhostrw, &
j1,j2,j3)
CALL acct_stop_inter
IF( peqopt == 1) THEN
!-----------------------------------------------------------------------
!
! Calculate the velocity divergence using newly updated velocity.
!
!-----------------------------------------------------------------------
tema = dzinv*dtsml1
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
div(i,j,k) = mapfct(i,j,7) &
*( u(i+1,j,k)*dxj3xm(i+1,j,k)- u(i ,j,k)*dxj3xm(i ,j,k) &
+ v(i,j+1,k)*dyj3ym(i,j+1,k)- v(i,j ,k)*dyj3ym(i,j ,k)) &
+( wcont(i,j,k+1)*aj3z(i,j,k+1) &
- wcont(i,j,k)*aj3z(i,j,k) ) * tema
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Compute the divergence term and the base state pressure advection
! in the pressure equation.
!
! The pressure divergence term is: -rhostr*csndsq*div/j3
! The base state pressure advection is -j3*w*d(pbar)/dzp = w*rhostr*g.
!
! These two terms are saved in pdiv.
!
!-----------------------------------------------------------------------
prgwg5=prgw * g05
tema = dtsml1*prgw * g05
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
pdiv(i,j,k)= &
tema*(w(i,j,k)+w(i,j,k+1))*rhostr(i,j,k) &
-csndsq(i,j,k)*rstj3i(i,j,k)*div(i,j,k)
END DO
END DO
END DO
ELSE
!-----------------------------------------------------------------------
!
! Calculate the density weighted mass divergence using newly
! updated velocity. Then calculate the pressure divergence term:
! -csndsq*div
!
!-----------------------------------------------------------------------
!
DO k= 2,nz-2
DO j= 1,ny-1
DO i= 1,nx
tem1(i,j,k)=u(i,j,k)*rhostru(i,j,k)*mapfct(i,j,5)
END DO
END DO
END DO
DO k= 2,nz-2
DO j= 1,ny
DO i= 1,nx-1
tem2(i,j,k)=v(i,j,k)*rhostrv(i,j,k)*mapfct(i,j,6)
END DO
END DO
END DO
DO k= 2,nz-1
DO j= 1,ny-1
DO i= 1,nx-1
tem3(i,j,k)=wcont(i,j,k)*rhostrw(i,j,k)
END DO
END DO
END DO
tema = dxinv*dtsml1
temb = dyinv*dtsml1
temc = dzinv*dtsml1
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
pdiv(i,j,k)= -csndsq(i,j,k)*( mapfct(i,j,7) &
*((tem1(i+1,j,k)-tem1(i,j,k))*tema &
+ (tem2(i,j+1,k)-tem2(i,j,k))*temb ) &
+ (tem3(i,j,k+1)-tem3(i,j,k))*temc )
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Integrate forward the pressure equation by one small timestep,
! using a backward (relative to the divergence term) integration
! scheme.
!
! And set lateral boundary conditions for the pressure equation.
!
!-----------------------------------------------------------------------
!
IF( nestgrd == 1 .AND. mgrid /= 1 ) THEN
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-2
pprt(i,j,k)=pprt(i,j,k) &
+pforce(i,j,k)+pdiv(i,j,k)*j3inv(i,j,k)
END DO
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
pprt( 1,j,k)=pprt( 1,j,k)+dtsml1 * pdtwb(j,k)
pprt(nx-1,j,k)=pprt(nx-1,j,k)+dtsml1 * pdteb(j,k)
END DO
END DO
DO k=2,nz-2
DO i=2,nx-2
pprt(i, 1,k)=pprt(i, 1,k)+dtsml1 * pdtsb(i,k)
pprt(i,ny-1,k)=pprt(i,ny-1,k)+dtsml1 * pdtnb(i,k)
END DO
END DO
!
! Call the pprt bottom boundary condition subroutine to
! compute the hydrostatic pprt at k=1.
!
CALL acct_interrupt(bc_acct)
CALL pprtbbc(nx,ny,nz,g05,buoy2nd,rhostr,pprt,ptprt, &
pbari,ptbari,phydro, &
tem1,tem2) ! tem1 = new pprt at k=1.
CALL acct_stop_inter
!
!-----------------------------------------------------------------------
!
! And top and bottom boundary conditions for the pressure equation.
!
!-----------------------------------------------------------------------
!
! bc's for mp not implemented for nested grids
CALL acct_interrupt(mp_acct)
CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, &
tem1(1,1,1),0,0,0,0 ,tbc,bbc)
CALL acct_stop_inter
ELSE ! Not nested grid
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
pprt(i,j,k)=pprt(i,j,k) &
+pforce(i,j,k)+pdiv(i,j,k)*j3inv(i,j,k)
END DO
END DO
END DO
IF(noisewrt == 1 .AND. MOD(ncalls(mptr),2) == 0 ) THEN
IF(lbcopt == 2) THEN
istart=MAX(2,(ngbrz+1))
iend=MIN((nx-2),(nx-ngbrz-1))
jstart=MAX(2,(ngbrz+1))
jend=MIN((ny-2),(ny-ngbrz-1))
ELSE
istart=2
iend=nx-2
jstart=2
jend=ny-2
END IF
sumdpdt=0.
DO k=2,nz-2
DO j=jstart,jend
DO i=istart,iend
sumdpdt=sumdpdt+abs(pforce(i,j,k)+pdiv(i,j,k)*j3inv(i,j,k))
END DO
END DO
END DO
sumdpdt=sumdpdt/ &
(float((iend-istart+1)*(jend-jstart+1)*(nz-3))*dtsml1)
write(nchdpdt1(mptr),'(f10.2,f12.6)') (curtim+2*dtsml1),sumdpdt
END IF
!
! Call the pprt bottom boundary condition subroutine to
! compute the hydrostatic pprt at k=1.
!
!call test_dump (pprt,"XXXsolve_pprt",nx,ny,nz,0,0)
CALL acct_interrupt(bc_acct)
CALL pprtbbc(nx,ny,nz,g05,buoy2nd,rhostr,pprt,ptprt, &
pbari,ptbari,phydro, &
tem1,tem2) ! tem1 = new pprt at k=1.
CALL acct_stop_inter
!call test_dump (pprt,"XXXAsolve_pprt",nx,ny,nz,0,1)
!
!-----------------------------------------------------------------------
!
! Apply the boundary conditions for the pressure equation.
!
! For the open boundary case, the boundary value is solved from the
! equation.
!
!-----------------------------------------------------------------------
!
!call test_dump (pprt,"XXX1solve_pprt",nx,ny,nz,0,0)
IF( lbcopt == 1 ) THEN ! Internal boundary conditions
ebc1=ebc
wbc1=wbc
sbc1=sbc
nbc1=nbc
IF( ebc == 4 ) ebc1=0
IF( wbc == 4 ) wbc1=0
IF( nbc == 4 ) nbc1=0
IF( sbc == 4 ) sbc1=0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(pprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mprecv2dew(pprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mpsend2dns(pprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
CALL mprecv2dns(pprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
END IF
CALL acct_interrupt(mp_acct)
CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, &
tem1(1,1,1), ebc1,wbc1,nbc1,sbc1,tbc,bbc)
CALL acct_stop_inter
ELSE ! External boundary condition
ebc1 = 3
wbc1 = 3
sbc1 = 3
nbc1 = 3
IF( ebc == 0 ) ebc1 = 0
IF( wbc == 0 ) wbc1 = 0
IF( nbc == 0 ) nbc1 = 0
IF( sbc == 0 ) sbc1 = 0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(pprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mprecv2dew(pprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mpsend2dns(pprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
CALL mprecv2dns(pprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
END IF
CALL acct_interrupt(mp_acct)
CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, &
tem1(1,1,1), 0,0,0,0,tbc,bbc)
CALL exbcp(nx,ny,nz, curtsml, pprt, &
exbcbuf(npr0exb),exbcbuf(nprdtexb))
CALL acct_stop_inter
END IF
END IF
!call test_dump (pprt,"XXX1Asolve_pprt",nx,ny,nz,0,1)
RETURN
END SUBROUTINE solvwpex
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SOLVWPIM ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE solvwpim(mptr, nx,ny,nz,exbcbufsz, dtsml1, curtsml, & 1,52
u,v,w,wcont,ptprt,pprt,phydro, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
rhostr,ptbar,ptbari,pbari,csndsq, &
wforce,wpgrad,pforce, &
x,y,z,zp, mapfct, j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
rhostru,rhostrv,rhostrw, &
exbcbuf,rhofct, &
fw,fp,wcontuv,tem1,tem2,tem3,tem4,j3irst, &
csj32irst,rstpbi,rstwi, &
dxij3xm,dyij3ym,pbzi)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Perform the time integration of w and pressure equations using
! implicit schem in vertical direction. The acoustically inactive
! forcing terms in these equations have been evaluated prior to this
! subroutine and are stored in wforce and pforce.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: M. Xue & H. Jin
! 8/20/93.
!
! 9/6/94 (M.Xue)
! A new version of vertically implicit small time step solver,
! The perturbation pressure contribution to buoyancy and the
! base state pressure advection terms are also treated implicitly
! in the vertical direction inside the small time steps.
!
! 9/1/94 (Y. Lu)
! Cleaned up documentation.
!
! 3/2/1995 (M. Xue and A. Shapiro)
! Significant bug fix in loop 600. rhostru there was mistakenly
! written as rhostrw.
!
! 10/31/95 (D. Weber)
! Added linear hydrostatic upper radiation condition for w-p.
! References are Klemp and Durran (JAS, 1983) and Chen (1991).
! Includes the use of trigs1,trigs2,ifax1,ifax2.
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 4/27/1996 (M. Xue)
! Corrected the formulations of the coefficients of the tridiagonal
! equations. The oringal formulation was slightly inacurate.
!
! 5/6/96 (M. Xue)
! Replaced csndsq in pressure buoyancy term by cpdcv*pbar/rhobar.
!
! 3/11/1997 (M. Xue and D. Weber)
! Corrected a minor bug related to the radiation top boundary
! condition.
!
! 07/22/97 (D. Weber)
! Added even fft for linear hydrostatic w-p upper radiation condition.
! References are Klemp and Durran (JAS, 1983) and Chen (1991).
! Includes the use of wsave1,wsave2,vwork1,vwork2 (fftopt=2).
!
! 10/21/97 (Donghai Wang)
! Using total density (rho) in the calculation of the pressure
! gradient force terms.
!
! 10/21/97 (Donghai Wang)
! Added the second order terms in the linerized buoyancy terms.
!
! 11/05/97 (D. Weber)
! Added phydro array for use in the bottom boundary condition for
! perturbation 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.
! Computed constants outside loops in selected loops.
!
! 9/18/98 (D. Weber)
! Added arrays aj3x,y,z,ptbari, and pbari to improve the speed
! of the code. Also added tem5-12 for pre-computing groups of
! commonly used constants.
!
! 3/18/99 (D. Weber)
! Bug fix to second order buoyancy term (do loop 350).
!
! 07/23/2001 (K. Brewster)
! Added mptr to argument list.
! Added optional writing of mean abs(dp/dt) as a diagnostic noise
! parameter.
!
!-----------------------------------------------------------------------
!
! 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
!
! dtsml1 Local value of small time step size
! curtsml Local curttim at a small time step
!
! u x component of velocity at tfuture (m/s)
! v y component of velocity at tfuture (m/s)
! w Vertical component of Cartesian velocity at tfuture
! (m/s)
! wcont Contravariant vertical velocity (m/s)
! ptprt Perturbation potential temperature at time tpresent (K)
! pprt Perturbation pressure at time tpresent (Pascal)
!
! wdteb Time tendency of the w field at the east boundary
! wdtwb Time tendency of the w field at the west boundary
! wdtnb Time tendency of the w field at the north boundary
! wdtsb Time tendency of the w field at the south boundary
!
! pdteb Time tendency of the pressure field at the east
! boundary
! pdtwb Time tendency of the pressure field at the west
! boundary
! pdtnb Time tendency of the pressure field at the north
! boundary
! pdtsb Time tendency of the pressure field at the south
! boundary
!
! phydro Big time step forcing term for use in computing the
! hydrostatic pressure at k=1.
!
! rhostr Base state density rhobar times j3 (kg/m**3)
! ptbar Base state potential temperature (K)
! ptbari Inverse base state potential temperature (K)
! pbari Inverse base state pressure (Pascal)
! csndsq Sound wave speed squared.
!
! wforce Acoustically inactive forcing terms in w-eq.
! (kg/(m*s)**2)
! wpgrad Pressure gradient term in w-eq. (kg/(m*s)**2)
! pforce Acoustically inactive terms in pressure eq. (Pascal/s)
!
! 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
! trigs1 Array containing pre-computed trig function for
! fftopt=1.
! trigs2 Array containing pre-computed trig function for
! fftopt=1.
! ifax1 Array containing the factors of nx for fftopt=1.
! ifax2 Array containing the factors of ny for fftopt=1.
!
! vwork1 2-D work array for fftopt = 2.
! vwork2 2-D work array for fftopt = 2.
! wsave1 Work array for fftopt = 2.
! wsave2 Work array for fftopt = 2.
!
! 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).
!
! OUTPUT:
!
! pprt Perturbation pressure at tfuture updated in time
! by dtsml1 (Pascal)
! w Vertical component of Cartesian velocity at tfuture
! updated in time by dtsml1 (m/s)
!
! WORK ARRAYS:
!
! fw Work array to carry force terms in w-eqation.
! fp Work array to carry force terms in p-eqation.
! wcontuv Work array to carry the contributions of u & v to wcont
! tem1 Work array
! tem2 Work array
! tem3 Work array
! tem4 Work array
! j3irst j3inv*rhostr
! csj32irst csndsq*j3inv*j3inv*rhostr
! rstpbi rhostr*pbari
! rstwi 1/rhostrw
! dxij3xm dxinv*aj3x*mapfct(i,j,5)
! dyij3ym dyinv*aj3y*mapfct(i,j,6)
! pbzi 1/avgz(pbar)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INTEGER :: mptr
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: dtsml1 ! Local small time step size (s)
REAL :: curtsml ! Local curttim at a small time step
REAL :: u (nx,ny,nz) ! u-velocity at tfuture (m/s)
REAL :: v (nx,ny,nz) ! v-velocity at tfuture (m/s)
REAL :: w (nx,ny,nz) ! w-velocity at tfuture (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure at tfuture
! (Pascal)
REAL :: wdteb (ny,nz) ! Time tendency of w at east boundary
REAL :: wdtwb (ny,nz) ! Time tendency of w at west boundary
REAL :: wdtnb (nx,nz) ! Time tendency of w at north boundary
REAL :: wdtsb (nx,nz) ! Time tendency of w at south boundary
REAL :: pdteb (ny,nz) ! Time tendency of pressure at east
! boundary
REAL :: pdtwb (ny,nz) ! Time tendency of pressure at west
! boundary
REAL :: pdtnb (nx,nz) ! Time tendency of pressure at north
! boundary
REAL :: pdtsb (nx,nz) ! Time tendency of pressure at south
! boundary
REAL :: phydro(nx,ny) ! Big time step forcing for computing
! hydrostatic pprt at k=1.
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: ptbari(nx,ny,nz) ! Inverse base state pot. temperature (K)
REAL :: pbari (nx,ny,nz) ! Inverse base state pressure (Pascal).
REAL :: csndsq(nx,ny,nz) ! Sound wave speed squared.
REAL :: wforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in w-momentum equation (kg/(m*s)**2)
REAL :: wpgrad(nx,ny,nz) ! Pressure gradient term in w-eq.
! (kg/(m*s)**2)
REAL :: pforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in the pressure equation (Pascal/s)
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, u and v 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 :: 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
! defined as d( zp )/d( z ).
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: rhofct(nx,ny,nz) ! rho-factor: rhobar/rho
REAL :: trigs1(3*(nx-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
REAL :: trigs2(3*(ny-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
INTEGER :: ifax1(13) ! Array containing the factors of nx for
! fftopt=1.
INTEGER :: ifax2(13) ! Array containing the factors of ny for
! fftopt=1.
REAL :: vwork1 (nx+1,ny+1) ! 2-D work array for fftopt = 2.
REAL :: vwork2 (ny,nx+1) ! 2-D work array for fftopt = 2.
REAL :: wsave1 (3*(ny-1)+15) ! Work array for fftopt = 2.
REAL :: wsave2 (3*(nx-1)+15) ! Work array for fftopt = 2.
REAL :: rhostru(nx,ny,nz) ! Averaged rhostr at u points (kg/m**3).
REAL :: rhostrv(nx,ny,nz) ! Averaged rhostr at v points (kg/m**3).
REAL :: rhostrw(nx,ny,nz) ! Averaged rhostr at w points (kg/m**3).
!
!-----------------------------------------------------------------------
!
! Temporary WORK ARRAYS:
!
!-----------------------------------------------------------------------
!
REAL :: fp (nx,ny,nz) ! Pressure gradient term in w-eq.
REAL :: fw (nx,ny,nz) ! Force in w equation.(kg/(m*s)**2)
REAL :: wcontuv(nx,ny,nz) ! Contributions of u & v to wcont
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 :: j3irst(nx,ny,nz) ! j3inv*rhostr
REAL :: csj32irst (nx,ny,nz) ! csndsq*j3inv*j3inv*rhostr
REAL :: rstpbi(nx,ny,nz) ! rhostr*pbari
REAL :: rstwi (nx,ny,nz) ! 1/rhostrw
REAL :: dxij3xm(nx,ny,nz) ! dxinv*aj3x*mapfct(i,j,5)
REAL :: dyij3ym(nx,ny,nz) ! dyinv*aj3y*mapfct(i,j,6)
REAL :: pbzi (nx,ny,nz) ! 1/avgz(pbar)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
INTEGER :: istart,iend,jstart,jend
INTEGER :: ebc1,wbc1,nbc1,sbc1
REAL :: g05,pk,nk,g05wp,g05pr
INTEGER :: wpprt, prgw, itema
REAL :: tema,temb,w1,w2,nrho
INTEGER :: buoy2swt !Switch for 1st-order or 2nd-order in buoyancy
REAL :: ptemk,ptemk1,pttemk,pttemk1
INTEGER :: noisewrt
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
INCLUDE 'phycst.inc' ! Physical constants
INCLUDE 'mp.inc' ! Message passing parameters.
REAL :: sumdpdt
CHARACTER (LEN=128 ) :: dpdtfn
INTEGER :: ldpdtfn,istat
INTEGER :: ncalls(mgrdmax), nchdpdt1(mgrdmax)
SAVE ncalls,nchdpdt1
DATA ncalls /mgrdmax*0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF( bsnesq == 1 ) THEN
wpprt = 0 ! Switch for pprt/(rhobar*csndsq) term in w-eq.
prgw = 0 ! Switch for rhobar*g*w term in w-eq.
ELSE
wpprt = 1 ! Switch for pprt/(rhobar*csndsq) term in w-eq.
prgw = 1 ! Switch for rhobar*g*w term in w-eq.
END IF
g05 = g*0.5
g05wp = g05*wpprt
g05pr = g05*prgw
noisewrt = 0
ncalls(mptr) = ncalls(mptr) + 1
IF (noisewrt == 1 .AND. ncalls(mptr) == 1) THEN
dpdtfn = runname(1:lfnkey)//'.dpdt'
ldpdtfn = 5 + lfnkey
IF(nestgrd == 1) THEN
WRITE(dpdtfn((ldpdtfn+1):(ldpdtfn+4)),'(a,i2.2)')'.g',mptr
ldpdtfn = ldpdtfn + 4
END IF
WRITE(6,'(1x,a,a,a/,1x,a)') &
'Check to see if file ',dpdtfn(1:ldpdtfn),' already exists.', &
'If so, append a version number to the filename.'
CALL fnversn( dpdtfn, ldpdtfn )
CALL getunit ( nchdpdt1(mptr) )
OPEN(nchdpdt1(mptr),FORM='formatted',STATUS='new', &
FILE=dpdtfn(1:ldpdtfn),IOSTAT=istat)
IF( istat /= 0) THEN
WRITE(6,'(/a,i2,/a/)') &
' Error occured when opening file '//dpdtfn(1:ldpdtfn)// &
' using FORTRAN unit ',nchdpdt1(mptr), &
' Program stopped in SOLVWPIM.'
CALL arpsstop('arpsstop called from SOLVWPIM error opening '// &
'file',1)
END IF
WRITE(nchdpdt1(mptr),'(a)') runname
WRITE(nchdpdt1(mptr),'(t2,a,t15,a)') 'Time (s)','MnAbsdpdt'
END IF
!-----------------------------------------------------------------------
!
! Calculate the horizontal velocity divergence using newly updated
! u and v velocity plus half vertical divergence from wcont, and
! store the result in tem3.
!
! Namely, tem3 = difx(u*avgx(j3))+dify(v*avgy(j3))+
! (1-tacoef)*difz(wcont*avgz(j3))
!
! note tem9=dtsml*aj3x*mapfct(5)*dxinv
! note tem10=dtsml*aj3y*mapfct(6)*dyinv
!-----------------------------------------------------------------------
tema = dtsml1*(1.0-tacoef)*dzinv
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem3(i,j,k) = mapfct(i,j,7) &
* ( (u(i+1,j,k)*dxij3xm(i+1,j,k) &
-u(i,j,k)*dxij3xm(i,j,k)) &
+ (v(i,j+1,k)*dyij3ym(i,j+1,k) &
-v(i,j,k)*dyij3ym(i,j,k)) ) &
+(wcont(i,j,k+1)*aj3z(i,j,k+1)-wcont(i,j,k)*aj3z(i,j,k)) &
*tema
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Compute the forcing terms in pressure equation
!
! fp = dtsml/j3 *(pforce-rhobar*csndsq*tem3+(1-tacoef)*rhostr*g*w)
!
! note pforce includes dtsml1 and j3inv...
!
! rhostr*g*w is the base state pressure advection term.
!
!-----------------------------------------------------------------------
!
tema = dtsml1*g05pr*(1.0-tacoef)
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
fp(i,j,k)=pforce(i,j,k) &
-csj32irst(i,j,k)*tem3(i,j,k) &
+tema*(w(i,j,k)+w(i,j,k+1))*j3irst(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Compute the first two terms in the contravariant vertical
! velocity formulation:
!
! wcontuv = (avgx(avgz(ustr)*j1)+avgy(avgz(vstr)*j2))/rhostrw
!
!-----------------------------------------------------------------------
!
IF( ternopt == 0 ) THEN
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
wcontuv(i,j,k)=0.0
END DO
END DO
END DO
ELSE
DO k= 1,nz-1
DO j= 1,ny-1
DO i= 1,nx
tem1(i,j,k)=u(i,j,k)*rhostru(i,j,k)
END DO
END DO
END DO
DO k= 1,nz-1
DO j= 1,ny
DO i= 1,nx-1
tem2(i,j,k)=v(i,j,k)*rhostrv(i,j,k)
END DO
END DO
END DO
DO k= 2,nz-1
DO j= 1,ny-1
DO i= 1,nx-1
wcontuv(i,j,k)=((tem1(i ,j,k)+tem1(i ,j,k-1))*j1(i ,j,k) &
+(tem1(i+1,j,k)+tem1(i+1,j,k-1))*j1(i+1,j,k) &
+(tem2(i ,j,k)+tem2(i ,j,k-1))*j2(i ,j,k) &
+(tem2(i,j+1,k)+tem2(i,j+1,k-1))*j2(i,j+1,k)) &
*mapfct(i,j,8) * rstwi(i,j,k)
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
wcontuv(i,j,nz-1)=0.0
wcontuv(i,j,2) = ((u(i ,j,2)+u(i ,j,1))*j1(i,j,2) &
+(u(i+1,j,2)+u(i+1,j,1))*j1(i+1,j,2) &
+(v(i,j ,2)+v(i,j ,1))*j2(i,j,2) &
+(v(i,j+1,2)+v(i,j+1,1))*j2(i,j+1,2)) &
*mapfct(i,j,8)
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Average rhofct at w points, stored in tem4
!
!-----------------------------------------------------------------------
!
CALL avgsw(rhofct,nx,ny,nz, 1,nx-1, 1,ny-1,tem4)
!
!-----------------------------------------------------------------------
!
! Compute the right-hand-side forcing term in tridiagonal linear
! equation. Array w is used to store this forcing term.
!
! First, we add the contribution of pressure perturbation to the
! buoyancy to wforce and wpgrad.
!
! This term has a form of - rhostr*g*pprt/(cpdcv*pbar),
!
!-----------------------------------------------------------------------
!
tema = g05wp/cpdcv
IF ( buoy2nd == 0) THEN !1st-order
buoy2swt = 0 !Switch for 1st-order or 2nd-order
ELSE !2nd-order
buoy2swt = 1
END IF
temb = 0.5*buoy2swt*(1.0-cpdcv)/cpdcv
DO k=3,nz-2
DO j=1,ny-1
DO i=1,nx-1
ptemk = pprt(i,j,k )*pbari(i,j,k ) ! new code...
ptemk1= pprt(i,j,k-1)*pbari(i,j,k-1)
fw(i,j,k) = wforce(i,j,k)-tem4(i,j,k)*(wpgrad(i,j,k) &
+rhostrw(i,j,k)*(ptemk+ptemk1+ &
temb*(ptemk*ptemk+ptemk1*ptemk1)) &
*tema)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! When potential temperature equation is solved within small time
! steps, the contribution from ptprt to buoyancy term is calculated
! here.
!
!-----------------------------------------------------------------------
!
IF( ptsmlstp == 1 ) THEN
tema = 1.0/(2.0*cpdcv)
DO k=3,nz-2
DO j=1,ny-1
DO i=1,nx-1
ptemk = pprt(i,j,k )*pbari(i,j,k )
ptemk1= pprt(i,j,k-1)*pbari(i,j,k-1)
pttemk = ptprt(i,j,k )*ptbari(i,j,k )
pttemk1= ptprt(i,j,k-1)*ptbari(i,j,k-1)
fw(i,j,k)=fw(i,j,k)+ &
rhostrw(i,j,k)*(pttemk+pttemk1 &
-buoy2swt*(pttemk*pttemk+pttemk1*pttemk1 &
+tema*(ptemk*pttemk + ptemk1*pttemk1))) &
*g05*tem4(i,j,k)
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! tem3=fp-dtsml*rhostr*csndsq*tacoef* d(wcontuv)/dz /(j3*j3)
!
!-----------------------------------------------------------------------
!
tema = tacoef*dtsml1*dzinv
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem3(i,j,k)=fp(i,j,k)-tema*csj32irst(i,j,k) &
*(wcontuv(i,j,k+1)-wcontuv(i,j,k))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! fw=fw-tacoef*d(tem3)/dz-g*tacoef*avgz(rhostr*tem3/(cpdcv*pbar)))
!
!-----------------------------------------------------------------------
!
tema = g05wp/cpdcv
DO k=3,nz-2
DO j=1,ny-1
DO i=1,nx-1
fw(i,j,k)=fw(i,j,k)-tacoef*tem4(i,j,k)*( &
(tem3(i,j,k)-tem3(i,j,k-1))*dzinv + &
(rstpbi(i,j,k )*tem3(i,j,k ) &
+rstpbi(i,j,k-1)*tem3(i,j,k-1))*tema)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! fw=fw*dtsml/rhostr + w
!
!-----------------------------------------------------------------------
!
DO k=3,nz-2
DO j=1,ny-1
DO i=1,nx-1
fw(i,j,k)=dtsml1*(fw(i,j,k)*rstwi(i,j,k)) + w(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Prepare the left-hand-side coefficents for tridiagonal
! equation set:
!
! A(k)*w(k-1)+B(k)*w(k)+C(k)*w(k+1)=D(k) for k=3,nz-2
!
! w(2) and w(nz-1) are used as the boundary conditions.
!
! Due to the lack of work arrays, we are storing the coefficients
! A in rhostru, B in rhostrv, and C in rhostrw. D is stored in fw.
!
! rhostru, rhostrv and rhostrw are re-calculated from rhostr after
! they have been used.
!
!-----------------------------------------------------------------------
!
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)= g05pr*j3irst(i,j,k)
tem2(i,j,k)= dzinv*csj32irst(i,j,k)
END DO
END DO
END DO
tema = (dtsml1*tacoef)**2 * wpprt*g/cpdcv
temb = (dtsml1*tacoef)**2 * dzinv
DO k=3,nz-2
DO j=1,ny-1
DO i=1,nx-1
pk = tem4(i,j,k)*tema*pbzi(i,j,k)
nk = tem4(i,j,k)*temb*rstwi(i,j,k)
rhostru(i,j,k)= (-nk+pk)*(tem1(i,j,k-1)+tem2(i,j,k-1))
rhostrw(i,j,k)= ( nk+pk)*(tem1(i,j,k )-tem2(i,j,k ))
rhostrv(i,j,k)= 1 &
+nk*(tem1(i,j,k)+tem2(i,j,k)-tem1(i,j,k-1)+tem2(i,j,k-1)) &
+pk*(tem1(i,j,k)+tem2(i,j,k)+tem1(i,j,k-1)-tem2(i,j,k-1))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Reset fw on the boundaries using the top and bottom boundary
! conditions of w.
!
! w=0.0 at k=nz-1. At the lower boundary, wcont=0.0, therefore
! w=-wcontuv, where wcontuv was calculated in loop 240.
!
!-----------------------------------------------------------------------
!
DO j=1,ny-1
DO i=1,nx-1
fw(i,j,3) =fw(i,j,3)+wcontuv(i,j,2)*rhostru(i,j,3)
fw(i,j,nz-2)=fw(i,j,nz-2)+0.0*rhostrw(i,j,nz-2)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Call the tridiagonal solver for either a rigid (tbc=1) or open upper
! boundary condition for w (tbc=4).
!
! For tbc=4, we have a choice of fft's for the upper boundary:
! fftopt =1, periodic fft for linearized hydrostatic radiation condition.
! fftopt =2, even fft for linearized hydrostatic radiation condition.
!
! The references are Klemp and Durran Jas (1983) and Chen (MWR) 1991.
!
! w1 is the coef for w(nz-1) and w2 is the coef. for w(nz-2) in the
! pressure equation. The pressure equation is solved at p(nz-2)
! for w(nz-1). tem2(i,j,nz-2) is the summation of the known terms
! in the pressure equation. The only unknowns are p(nz-2), w(nz-1)
! and w(nz-1) at the new time level. In the tridiagonal solver
! the elimination step is performed and w(nz-2) is given in termsx
! of w(nz-1) and known terms. THe next step is to use the rad.
! condition (in fourier space):
!
!
! p = N * rhobar * w(nz-1) / abs(kx+ky)
!
! The end result is a relation for w(nz-1) given known quantities.
! Trigs1 and trigs2 are the predetermined trig functions used in the
! fft program in tridiag.
!
!-----------------------------------------------------------------------
IF(tbc == 4)THEN ! apply linear radiation condition to w at nz-1.
tema = rhostr(nx/2,ny/2,nz-2)*csndsq(nx/2,ny/2,nz-2)* &
j3inv(nx/2,ny/2,nz-2)
temb = dtsml1*j3inv(nx/2,ny/2,nz-2)
w2=temb*tacoef*(tema*dzinv+g05*prgw*rhostr(nx/2,ny/2,nz-2))
w1=temb*tacoef*(-tema*dzinv+g05*prgw*rhostr(nx/2,ny/2,nz-2))
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,2)=pprt(i,j,nz-2)+tem3(i,j,nz-2)
END DO
END DO
nrho=SQRT(g/ptbar(nx/2,ny/2,nz-2)*(ptbar(nx/2,ny/2,nz-1) &
-ptbar(nx/2,ny/2,nz-3))*dzinv*0.5) &
*rhostr(nx/2,ny/2,nz-2)*j3inv(nx/2,ny/2,nz-2)
END IF
!
!-----------------------------------------------------------------------
!
! NOTE: tem2(1,1,4) is passed into subroutine tridiag and becomes
! a 2-D array of size (nx+1,ny+1).
!
!-----------------------------------------------------------------------
!
CALL tridiag(nx,ny,nz,rhostru,rhostrv,rhostrw,fw,tem1, &
tem2(1,1,1),tem2(1,1,2),w1,w2,nrho,tem2(1,1,3),trigs1,trigs2, &
ifax1,ifax2,wsave1,wsave2,vwork1,vwork2,tem2(1,1,4))
!
!-----------------------------------------------------------------------
!
! Restore rhostru, rhostrv and rhostrw after they have been used
! are work arrays.
!
!-----------------------------------------------------------------------
!
CALL rhouvw(nx,ny,nz,rhostr,rhostru,rhostrv,rhostrw)
!
!-----------------------------------------------------------------------
!
! On exit of tridiag, the interior solution of w is stored in fw.
!
!-----------------------------------------------------------------------
IF(tbc == 4)THEN ! set itema = nz-1
itema = nz-1
ELSE ! set itema = nz-2
itema = nz-2
END IF
IF( nestgrd /= 1 .OR. mgrid == 1 ) THEN
! For non-nesting case or the
! base-grid of the nested case
DO k=3,itema
DO j=1,ny-1
DO i=1,nx-1
w(i,j,k) = fw(i,j,k)
END DO
END DO
END DO
!call test_dump (pprt,"XXX2solve_pprt",nx,ny,nz,0,0)
IF ( lbcopt == 1 ) THEN ! Internal boundary condition
ebc1=ebc
wbc1=wbc
sbc1=sbc
nbc1=nbc
IF( ebc == 4 ) ebc1=0
IF( wbc == 4 ) wbc1=0
IF( nbc == 4 ) nbc1=0
IF( sbc == 4 ) sbc1=0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mprecv2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mpsend2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
CALL mprecv2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
END IF
CALL acct_interrupt(mp_acct)
CALL lbcw(nx,ny,nz,dtsml1,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, &
ebc1,wbc1,nbc1,sbc1)
CALL acct_stop_inter
ELSE ! External boundary condition
ebc1 = 3
wbc1 = 3
sbc1 = 3
nbc1 = 3
IF( ebc == 0 ) ebc1 = 0
IF( wbc == 0 ) wbc1 = 0
IF( nbc == 0 ) nbc1 = 0
IF( sbc == 0 ) sbc1 = 0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mprecv2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mpsend2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
CALL mprecv2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
END IF
CALL acct_interrupt(mp_acct)
CALL lbcw(nx,ny,nz,dtsml1,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, &
ebc1,wbc1,nbc1,sbc1)
CALL acct_stop_inter
END IF
ELSE ! For nested interior grid
DO k=3,itema
DO j=2,ny-2
DO i=2,nx-2
w(i,j,k) = fw(i,j,k)
END DO
END DO
END DO
DO k=3,nz-2
DO j=1,ny-1
w( 1,j,k)=w( 1,j,k)+dtsml1 * wdtwb(j,k)
w(nx-1,j,k)=w(nx-1,j,k)+dtsml1 * wdteb(j,k)
END DO
END DO
DO k=3,nz-2
DO i=2,nx-2
w(i, 1,k)=w(i, 1,k)+dtsml1 * wdtsb(i,k)
w(i,ny-1,k)=w(i,ny-1,k)+dtsml1 * wdtnb(i,k)
END DO
END DO
END IF
!call test_dump (pprt,"XXX2Asolve_pprt",nx,ny,nz,0,1)
!
!-----------------------------------------------------------------------
!
! Calculate wcont at time tfuture, including the boundary
! points. Wcont at the lateral boundaries is calculated
! from boundary u, v and w. Wcont at the top and bottom
! depends on the boundary condition option.
!
!-----------------------------------------------------------------------
!
!call test_dump (w,"XXX2sBwc_w1",nx,ny,nz,1,0)
CALL wcontra(nx,ny,nz,u,v,w,mapfct,j1,j2,j3,aj3z, &
rhostr,rhostru,rhostrv,rhostrw,wcont,tem1,tem2)
!
!-----------------------------------------------------------------------
!
! Set the top and bottom boundary conditions for w based on u, v and
! wcont at the top and bottom boundaries.
!
!-----------------------------------------------------------------------
!
CALL acct_interrupt(bc_acct)
CALL vbcw(nx,ny,nz,w,wcont,tbc,bbc,u,v, &
rhostr,rhostru,rhostrv,rhostrw, &
j1,j2,j3)
CALL acct_stop_inter
!
!-----------------------------------------------------------------------
!
! Calculate the new pressure
!
! pprt(new) = pprt(old)+fp+dtsml*tacoef/j3*
! (rhostr*g*avg(w)-rhostr/j3*csndsq*difz(j3*wcont))
!
!-----------------------------------------------------------------------
tema = tacoef*dtsml1*g05pr
temb = tacoef*dtsml1*dzinv
IF( nestgrd == 1 .AND. mgrid /= 1 ) THEN
! For nested interior grid
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-2
pprt(i,j,k)=pprt(i,j,k) + fp(i,j,k) &
+tema*j3irst(i,j,k)*(w(i,j,k+1)+w(i,j,k)) &
-temb*csj32irst(i,j,k)* &
(aj3z(i,j,k+1)*wcont(i,j,k+1)-aj3z(i,j,k)*wcont(i,j,k))
END DO
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
pprt( 1,j,k)=pprt( 1,j,k)+dtsml1 * pdtwb(j,k)
pprt(nx-1,j,k)=pprt(nx-1,j,k)+dtsml1 * pdteb(j,k)
END DO
END DO
DO k=2,nz-2
DO i=2,nx-2
pprt(i, 1,k)=pprt(i, 1,k)+dtsml1 * pdtsb(i,k)
pprt(i,ny-1,k)=pprt(i,ny-1,k)+dtsml1 * pdtnb(i,k)
END DO
END DO
!
! Call the pprt bottom boundary condition subroutine to
! compute the hydrostatic pprt at k=1.
!
CALL acct_interrupt(bc_acct)
CALL pprtbbc(nx,ny,nz,g05,buoy2swt,rhostr,pprt,ptprt, &
pbari,ptbari,phydro, &
tem1,tem2) ! tem1 = new pprt at k=1.
CALL acct_stop_inter
! bc's for mp not implemented for nested grids
CALL acct_interrupt(mp_acct)
CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, &
tem1(1,1,1), 0,0,0,0,tbc,bbc)
CALL acct_stop_inter
ELSE ! For non-nesting case or the
! base-grid of the nested case
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
pprt(i,j,k)=pprt(i,j,k) + fp(i,j,k) &
+tema*j3irst(i,j,k)*(w(i,j,k+1)+w(i,j,k)) &
-temb*csj32irst(i,j,k)* &
(aj3z(i,j,k+1)*wcont(i,j,k+1)-aj3z(i,j,k)*wcont(i,j,k))
END DO
END DO
END DO
IF (noisewrt == 1 .AND. MOD(ncalls(mptr),2) == 0 ) THEN
IF(lbcopt == 2) THEN
istart=MAX(2,(ngbrz+1))
iend=MIN((nx-2),(nx-ngbrz-1))
jstart=MAX(2,(ngbrz+1))
jend=MIN((ny-2),(ny-ngbrz-1))
ELSE
istart=2
iend=nx-2
jstart=2
jend=ny-2
END IF
sumdpdt=0.
DO k=2,nz-2
DO j=jstart,jend
DO i=istart,iend
sumdpdt=sumdpdt+abs(fp(i,j,k) &
+tema*j3irst(i,j,k)*(w(i,j,k+1)+w(i,j,k)) &
-temb*csj32irst(i,j,k)* &
(aj3z(i,j,k+1)*wcont(i,j,k+1)-aj3z(i,j,k)*wcont(i,j,k)))
END DO
END DO
END DO
sumdpdt=sumdpdt/ &
(float((iend-istart+1)*(jend-jstart+1)*(nz-3))*dtsml1)
write(nchdpdt1(mptr),'(f10.2,f12.6)') (curtim+2*dtsml1),sumdpdt
END IF
!
! Call the pprt bottom boundary condition subroutine to
! compute the hydrostatic pprt at k=1.
!
CALL acct_interrupt(bc_acct)
CALL pprtbbc(nx,ny,nz,g05,buoy2swt,rhostr,pprt,ptprt, &
pbari,ptbari,phydro, &
tem1,tem2) ! tem1 = new pprt at k=1.
CALL acct_stop_inter
!
!-----------------------------------------------------------------------
!
! Apply the boundary conditions for the pressure equation.
!
! For the open boundary case, the boundary value of pprt is
! predicted by the pressure equation.
!
!-----------------------------------------------------------------------
!
!call test_dump (pprt,"XXX3solve_pprt",nx,ny,nz,0,0)
IF( lbcopt == 1 ) THEN ! Internal boundary conditions
ebc1=ebc
wbc1=wbc
sbc1=sbc
nbc1=nbc
IF( ebc == 4 ) ebc1=0
IF( wbc == 4 ) wbc1=0
IF( nbc == 4 ) nbc1=0
IF( sbc == 4 ) sbc1=0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(pprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mprecv2dew(pprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mpsend2dns(pprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
CALL mprecv2dns(pprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
END IF
CALL acct_interrupt(mp_acct)
CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, &
tem1(1,1,1), ebc1,wbc1,nbc1,sbc1,tbc,bbc)
CALL acct_stop_inter
ELSE ! External boundary condition
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(pprt,nx,ny,nz,0,0,0,mptag,tem2)
CALL mprecv2dew(pprt,nx,ny,nz,0,0,0,mptag,tem2)
CALL mpsend2dns(pprt,nx,ny,nz,0,0,0,mptag,tem2)
CALL mprecv2dns(pprt,nx,ny,nz,0,0,0,mptag,tem2)
END IF
CALL acct_interrupt(mp_acct)
CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, &
tem1(1,1,1), 0,0,0,0,tbc,bbc)
CALL exbcp(nx,ny,nz, curtsml, pprt, &
exbcbuf(npr0exb),exbcbuf(nprdtexb))
CALL acct_stop_inter
END IF
END IF
!call test_dump (pprt,"XXX3Asolve_pprt",nx,ny,nz,0,1)
RETURN
END SUBROUTINE solvwpim
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SOLVWPIM1 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE solvwpim1(nx,ny,nz,exbcbufsz, dtsml1, curtsml, & 1,49
u,v,w,wcont,ptprt,pprt,phydro, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
rhostr,ptbar,ptbari,pbari,csndsq, &
wforce,wpgrad,pforce, &
x,y,z,zp, mapfct, j1,j2,j3,aj3z,j3inv, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
rhostru,rhostrv,rhostrw, &
exbcbuf,rhofct, &
fw,fp,tem1,tem2,tem3,tem4,tem5, &
rstpbi,rstptbi,rstwi,pbzi,rstj3i)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Perform the time integration of w and pressure equations using
! implicit scheme in vertical direction. The acoustically inactive
! forcing terms in these equations have been evaluated prior to this
! subroutine and are stored in wforce and pforce.
!
! This routine is for peqopt=2.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: M. Xue
! 4/26/1996.
! Written for peqopt=2 based on SOLVWPIM.
!
! MODIFICATION HISTORY:
!
! 07/22/97 (D. Weber)
! Added wsave1,wsave2,vwork1,vwork2 for use in the even fft version
! of the upper w-p radiation condition (fftopt=2).
!
! 10/21/97 (Donghai Wang)
! Using total density (rho) in the calculation of the pressure
! gradient force terms.
!
! 10/21/97 (Donghai Wang)
! Added the second order terms in the linerized buoyancy terms.
!
! 11/05/97 (D. Weber)
! Added phydro array for use in the bottom boundary condition for
! perturbation 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.
! Removed multiplication of constants from loops.
!
! 9/18/98 (D. Weber)
! Added arrays aj3z ptbari,pbari,rstpbi,rstptbi,rstwi,pbzi,
! rstj3i to improve code efficiency.
!
!-----------------------------------------------------------------------
!
! 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
!
! dtsml1 Local value of small time step size
! curtsml Local curttim at a small time step
!
! u x component of velocity at tfuture (m/s)
! v y component of velocity at tfuture (m/s)
! w Vertical component of Cartesian velocity at tfuture
! (m/s)
! wcont Contravariant vertical velocity (m/s)
! ptprt Perturbation potential temperature at time tpresent (K)
! pprt Perturbation pressure at time tpresent (Pascal)
!
! wdteb Time tendency of the w field at the east boundary
! wdtwb Time tendency of the w field at the west boundary
! wdtnb Time tendency of the w field at the north boundary
! wdtsb Time tendency of the w field at the south boundary
!
! pdteb Time tendency of the pressure field at the east
! boundary
! pdtwb Time tendency of the pressure field at the west
! boundary
! pdtnb Time tendency of the pressure field at the north
! boundary
! pdtsb Time tendency of the pressure field at the south
! boundary
!
! phydro Big time step forcing term for use in computing the
! hydrostatic pressure at k=1.
!
! rhostr Base state density rhobar times j3 (kg/m**3)
! ptbar Base state potential temperature (K)
! ptbari Inverse base state potential temperature (K)
! pbari Inverse base state pressure (Pascal)
! csndsq Sound wave speed squared.
!
! wforce Acoustically inactive forcing terms in w-eq.
! (kg/(m*s)**2)
! wpgrad Pressure gradient term in w-eq. (kg/(m*s)**2)
! pforce Acoustically inactive terms in pressure eq. (Pascal/s)
!
! 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
! aj3z Avgz of the coordinate transformation Jacobian d(zp)/dz
! trigs1 Array containing pre-computed trig function for
! fftopt=1.
! trigs2 Array containing pre-computed trig function for
! fftopt=1.
! ifax1 Array containing the factors of nx for fftopt=1.
! ifax2 Array containing the factors of ny for fftopt=1.
!
! vwork1 2-D work array for fftopt=2.
! vwork2 2-D work array for fftopt=2.
! wsave1 Work array for fftopt=2.
! wsave2 Work array for fftopt=2.
!
! 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).
!
! OUTPUT:
!
! pprt Perturbation pressure at tfuture updated in time
! by dtsml1 (Pascal)
! w Vertical component of Cartesian velocity at tfuture
! updated in time by dtsml1 (m/s)
!
! WORK ARRAYS:
!
! fw Work array to carry force terms in w-eqation.
! fp Work array to carry force terms in p-eqation.
! tem1 Work array
! tem2 Work array
! tem3 Work array
! tem4 Work array
! tem5 Work array
! rstpbi pbari*rhostr
! rstptbi ptbari*rhostr
! rstwi 1/rhostrw
! pbzi 2/avgz(pbar)
! rstj3i rhostr*j3inv
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: dtsml1 ! Local small time step size (s)
REAL :: curtsml ! Local curttim at a small time step
REAL :: u (nx,ny,nz) ! u-velocity at tfuture (m/s)
REAL :: v (nx,ny,nz) ! v-velocity at tfuture (m/s)
REAL :: w (nx,ny,nz) ! w-velocity at tfuture (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure at tfuture
! (Pascal)
REAL :: wdteb (ny,nz) ! Time tendency of w at east boundary
REAL :: wdtwb (ny,nz) ! Time tendency of w at west boundary
REAL :: wdtnb (nx,nz) ! Time tendency of w at north boundary
REAL :: wdtsb (nx,nz) ! Time tendency of w at south boundary
REAL :: pdteb (ny,nz) ! Time tendency of pressure at east
! boundary
REAL :: pdtwb (ny,nz) ! Time tendency of pressure at west
! boundary
REAL :: pdtnb (nx,nz) ! Time tendency of pressure at north
! boundary
REAL :: pdtsb (nx,nz) ! Time tendency of pressure at south
! boundary
REAL :: phydro(nx,ny) ! Big time step forcing for computing
! hydrostatic pprt at k=1.
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: ptbari(nx,ny,nz) ! Inverse base state pot. temperature (K)
REAL :: pbari (nx,ny,nz) ! Inverse base state pressure (Pascal).
REAL :: csndsq(nx,ny,nz) ! Sound wave speed squared.
REAL :: wforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in w-momentum equation (kg/(m*s)**2)
REAL :: wpgrad(nx,ny,nz) ! Pressure gradient term in w-eq.
! (kg/(m*s)**2)
REAL :: pforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in the pressure equation (Pascal/s)
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, u and v 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 :: 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
! defined as d( zp )/d( z ).
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: rhofct(nx,ny,nz) ! rho-factor: rhobar/rho
REAL :: trigs1(3*(nx-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
REAL :: trigs2(3*(ny-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
INTEGER :: ifax1(13) ! Array containing the factors of nx
! for fftopt=1.
INTEGER :: ifax2(13) ! Array containing the factors of ny
! for fftopt=1.
REAL :: vwork1 (nx+1,ny+1) ! 2-D work array for fftopt=2 option.
REAL :: vwork2 (ny,nx+1) ! 2-D work array for fftopt=2 option.
REAL :: wsave1 (3*(ny-1)+15) ! Work array for fftopt=2 option.
REAL :: wsave2 (3*(nx-1)+15) ! Work array for fftopt=2 option.
REAL :: rhostru(nx,ny,nz) ! Averaged rhostr at u points (kg/m**3).
REAL :: rhostrv(nx,ny,nz) ! Averaged rhostr at v points (kg/m**3).
REAL :: rhostrw(nx,ny,nz) ! Averaged rhostr at w points (kg/m**3).
!
!-----------------------------------------------------------------------
!
! Temporary WORK ARRAYS:
!
!-----------------------------------------------------------------------
!
REAL :: fp (nx,ny,nz) ! Pressure gradient term in w-eq.
REAL :: fw (nx,ny,nz) ! Force in w 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
REAL :: tem4 (nx,ny,nz) ! Temporary work array
REAL :: tem5 (nx,ny,nz) ! Temporary work array
REAL :: rstpbi(nx,ny,nz) ! pbari*rhostr
REAL :: rstptbi(nx,ny,nz) ! ptbari*rhostr
REAL :: rstwi (nx,ny,nz) ! 1/rhostrw
REAL :: pbzi (nx,ny,nz) ! 2/avgz(pbar)
REAL :: rstj3i(nx,ny,nz) ! rhostr*j3inv
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
INTEGER :: ebc1,wbc1,nbc1,sbc1
REAL :: g05,pk,nk,g05wp
INTEGER :: wpprt, itema
REAL :: tema,temb,w1,w2,nrho
INTEGER :: prgw
INTEGER :: buoy2swt !Switch for 1st-order or 2nd-order in buoyancy
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
INCLUDE 'phycst.inc' ! Physical constants
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF( bsnesq == 1 ) THEN
wpprt = 0 ! Switch for pprt/(rhobar*csndsq) term in w-eq.
ELSE
wpprt = 1 ! Switch for pprt/(rhobar*csndsq) term in w-eq.
END IF
prgw = 0
g05 = g*0.5
g05wp = g05*wpprt
!
!-----------------------------------------------------------------------
!
! Calculate the horizontal velocity divergence using newly updated
! u and v velocity plus half vertical divergence from wcont, and
! store the result in tem3.
!
! Namely, tem3 = difx(u*avgx(j3))+dify(v*avgy(j3))+
! (1-tacoef)*difz(wcont*avgz(j3))
!
!-----------------------------------------------------------------------
!
DO k= 2,nz-2
DO j= 1,ny-1
DO i= 1,nx
tem1(i,j,k)=u(i,j,k)*rhostru(i,j,k)*mapfct(i,j,5)
END DO
END DO
END DO
DO k= 2,nz-2
DO j= 1,ny
DO i= 1,nx-1
tem2(i,j,k)=v(i,j,k)*rhostrv(i,j,k)*mapfct(i,j,6)
END DO
END DO
END DO
DO k= 2,nz-1
DO j= 1,ny-1
DO i= 1,nx-1
tem3(i,j,k)=wcont(i,j,k)*rhostrw(i,j,k)
END DO
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem4(i,j,k)=mapfct(i,j,7) &
*((tem1(i+1,j,k)-tem1(i,j,k))*dxinv &
+(tem2(i,j+1,k)-tem2(i,j,k))*dyinv) &
+(tem3(i,j,k+1)-tem3(i,j,k))*dzinv*(1.0-tacoef)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Compute the forcing terms in pressure equation
!
! fp = dtsml/j3 *(pforce-csndsq*tem4)
!
! rhostr*g*w is the base state pressure advection term.
!
!-----------------------------------------------------------------------
!
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
! fp(i,j,k)=dtsml1*j3inv(i,j,k) *
! : (pforce(i,j,k)-csndsq(i,j,k)*tem4(i,j,k))
fp(i,j,k)=pforce(i,j,k) -dtsml1*j3inv(i,j,k)* &
csndsq(i,j,k)*tem4(i,j,k)
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Compute two more terms in fp related to coordinate transformation:
!
! tem3 = (avgx(avgz(ustr)*j1)+avgy(avgz(vstr)*j2))/avgz(j3)
!
!-----------------------------------------------------------------------
!
IF( ternopt /= 0 ) THEN
DO k= 1,nz-1
DO j= 1,ny-1
DO i= 1,nx
tem1(i,j,k)=u(i,j,k)*rhostru(i,j,k)
END DO
END DO
END DO
DO k= 1,nz-1
DO j= 1,ny
DO i= 1,nx-1
tem2(i,j,k)=v(i,j,k)*rhostrv(i,j,k)
END DO
END DO
END DO
DO k= 2,nz-1
DO j= 1,ny-1
DO i= 1,nx-1
tem3(i,j,k)=((tem1(i ,j,k)+tem1(i ,j,k-1))*j1(i ,j,k) &
+(tem1(i+1,j,k)+tem1(i+1,j,k-1))*j1(i+1,j,k) &
+(tem2(i ,j,k)+tem2(i ,j,k-1))*j2(i ,j,k) &
+(tem2(i,j+1,k)+tem2(i,j+1,k-1))*j2(i,j+1,k)) &
*mapfct(i,j,8)/aj3z(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! fp=fp-(dtsml/j3)*csndsq*tacoef* d(tem3)/dz
!
!-----------------------------------------------------------------------
!
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
fp(i,j,k)=fp(i,j,k) &
-dtsml1*j3inv(i,j,k)*csndsq(i,j,k) &
*tacoef*(tem3(i,j,k+1)-tem3(i,j,k))*dzinv
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Average rhofct at w points, stored in tem5
!
!-----------------------------------------------------------------------
!
CALL avgsw(rhofct,nx,ny,nz, 1,nx-1, 1,ny-1,tem5)
!
!-----------------------------------------------------------------------
!
! Compute the right-hand-side forcing term in tridiagonal linear
! equation. Array w is used to store this forcing term.
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! fw=wforce-wpgrad-tacoef*d(fp)/dz
! -g*avgz(rhostr*(pprt+taceof*pfrc)/(gamma*pbar))
!
!-----------------------------------------------------------------------
!
IF ( buoy2nd == 0) THEN !1st-order
buoy2swt = 0 !Switch for 1st-order or 2nd-order
ELSE !2nd-order
buoy2swt = 1
END IF
tema = 0.5*buoy2swt*(1.0-cpdcv)/cpdcv
temb = 1.0/cpdcv
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
! tem1(i,j,k)=tem6(i,j,k)*((pprt(i,j,k)+tacoef*fp(i,j,k))
tem1(i,j,k)=rstpbi(i,j,k)*((pprt(i,j,k)+tacoef*fp(i,j,k)) &
+tema*pprt(i,j,k)*pprt(i,j,k) &
*pbari(i,j,k))*temb
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! When potential temperature equation is solved within small time
! steps, the contribution from ptprt to buoyancy term is calculated
! here.
!
!-----------------------------------------------------------------------
!
IF( ptsmlstp == 1 ) THEN
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
! tem2(i,j,k)= ptprt(i,j,k)*tem7(i,j,k)
tem2(i,j,k)= ptprt(i,j,k)*rstptbi(i,j,k) &
*(1.0-buoy2swt*ptprt(i,j,k)*ptbari(i,j,k))
END DO
END DO
END DO
ELSE
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)= 0.0
END DO
END DO
END DO
END IF
tema = tacoef*dzinv
DO k=3,nz-2
DO j=1,ny-1
DO i=1,nx-1
fw(i,j,k) = wforce(i,j,k)-tem5(i,j,k)*(wpgrad(i,j,k) &
+tema*(fp(i,j,k)-fp(i,j,k-1)) &
+g05wp*(tem1(i,j,k)+tem1(i,j,k-1)) &
-g05 *(tem2(i,j,k)+tem2(i,j,k-1)))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! fw=fw*dtsml/rhostr + w
!
!-----------------------------------------------------------------------
!
DO k=3,nz-2
DO j=1,ny-1
DO i=1,nx-1
! fw(i,j,k)=fw(i,j,k)*dtsml1*tem8(i,j,k) + w(i,j,k)
fw(i,j,k)=fw(i,j,k)*dtsml1*rstwi(i,j,k) + w(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Prepare the left-hand-side coefficents for tridiagonal
! equation set:
!
! A(k)*w(k-1)+B(k)*w(k)+C(k)*w(k+1)=D(k) for k=3,nz-2
!
! w(2) and w(nz-1) are used as the boundary conditions.
!
! Due to the lack of work arrays, we are storing the coefficients
! A in rhostru, B in rhostrv, and C in rhostrw. D is stored in fw.
!
! rhostru, rhostrv and rhostrw are re-calculated from rhostr after
! they have been used.
!
!-----------------------------------------------------------------------
!
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
! tem4(i,j,k)=0.5*(tem10(i,j,k) +tem10(i,j,k-1))
tem4(i,j,k)=0.5*(rstj3i(i,j,k)+ rstj3i(i,j,k-1))
END DO
END DO
END DO
tema = (dtsml1*tacoef)**2 * wpprt*g*dzinv /cpdcv
temb = (dtsml1*tacoef*dzinv)**2
DO k=3,nz-2
DO j=1,ny-1
DO i=1,nx-1
! pk = tem5(i,j,k)*tema*tem9(i,j,k)
pk = tem5(i,j,k)*tema*pbzi(i,j,k)
nk = tem5(i,j,k)*temb*rstwi(i,j,k)
! nk = tem5(i,j,k)*temb*tem8(i,j,k)
tem1(i,j,k)= (-nk+pk)*csndsq(i,j,k-1)*j3inv(i,j,k-1)
tem3(i,j,k)=-( nk+pk)*csndsq(i,j,k )*j3inv(i,j,k )
tem2(i,j,k)= 1-tem4(i,j,k)*(tem1(i,j,k)+tem3(i,j,k))
tem1(i,j,k)=tem1(i,j,k)*tem4(i,j,k-1)
tem3(i,j,k)=tem3(i,j,k)*tem4(i,j,k+1)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Reset fw on the boundaries using the top and bottom boundary
! conditions of w.
!
!-----------------------------------------------------------------------
!
DO j=1,ny-1
DO i=1,nx-1
w(i,j,nz-1)=0.0
w(i,j,2) =-((u(i ,j,2)+u(i ,j,1))*j1(i,j,2) &
+(u(i+1,j,2)+u(i+1,j,1))*j1(i+1,j,2) &
+(v(i,j ,2)+v(i,j ,1))*j2(i,j,2) &
+(v(i,j+1,2)+v(i,j+1,1))*j2(i,j+1,2)) &
*mapfct(i,j,8)
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
fw(i,j,3) =fw(i,j, 3)-w(i,j, 2)*tem1(i,j, 3)
fw(i,j,nz-2)=fw(i,j,nz-2)-w(i,j,nz-1)*tem3(i,j,nz-2)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Call the tridiagonal solver for either a rigid or open upper
! boundary condition for w.
!
! tbc = 4, periodic fft for linearized hydrostatic radiation condition.
! tbc = 6, even fft for linearized hydrostatic radiation condition.
! The references are Klemp and Durran Jas (1983) and Chen (MWR) 1991.
!
! w1 is the coef for w(nz-1) and w2 is the coef. for w(nz-2 in the
! pressure equation. The pressure equation is solved at p(nz-2
! for w(nz-1). tem4(i,j,nz-2) is the summation of the known terms
! in the pressure equation. The only unknowns are p(nz-2), w(nz-1)
! and w(nz-1) at the new time level. In the tridiagonal solver
! the elimination step is performed and w(nz-2) is given in termsx
! of w(nz-1) and known terms. THe next step is to use the rad.
! condition (in fourier space):
!
!
! p = N * rhobar * w(nz-1) / abs(kx+ky)
!
! The end result is a relation for w(nz-1) given known quantities.
! Trigs1 and trigs2 are the predetermined trig functions used in the
! fft program in tridiag.
!
!-----------------------------------------------------------------------
IF(tbc == 4)THEN ! apply linear radiation condition to w at nz-1.
tema = rhostr(nx/2,ny/2,nz-2)*csndsq(nx/2,ny/2,nz-2)* &
j3inv(nx/2,ny/2,nz-2)
temb = dtsml1*j3inv(nx/2,ny/2,nz-2)
w2=temb*tacoef*(tema*dzinv+g05*prgw*rhostr(nx/2,ny/2,nz-2))
w1=temb*tacoef*(-tema*dzinv+g05*prgw*rhostr(nx/2,ny/2,nz-2))
DO j=1,ny-1
DO i=1,nx-1
tem4(i,j,2)=pprt(i,j,nz-2)+fp(i,j,nz-2)
END DO
END DO
nrho=SQRT(g/ptbar(nx/2,ny/2,nz-2)*(ptbar(nx/2,ny/2,nz-1) &
-ptbar(nx/2,ny/2,nz-3))*dzinv*0.5) &
*rhostr(nx/2,ny/2,nz-2)*j3inv(nx/2,ny/2,nz-2)
END IF
!
!-----------------------------------------------------------------------
!
! NOTE: tem4(1,1,4) is passed into subroutine tridiag and becomes
! a 2-D array of size (nx+1,ny+1).
!
! rhostrw is used as a work array in tridiag.
!
!-----------------------------------------------------------------------
!
CALL tridiag(nx,ny,nz,tem1,tem2,tem3,fw,rhostrw, &
tem4(1,1,1),tem4(1,1,2),w1,w2,nrho,tem4(1,1,3), &
trigs1,trigs2,ifax1,ifax2,wsave1,wsave2,vwork1,vwork2, &
tem4(1,1,4))
CALL avgsw(rhostr,nx,ny,nz, 1,nx-1, 1,ny-1, rhostrw)
! call TRIDIAG_old(nx,ny,nz,tem1,tem2,tem3,fw, tem4)
!
!-----------------------------------------------------------------------
!
! On exit of tridiag, the interior solution of w is stored in fw.
!
!-----------------------------------------------------------------------
IF(tbc == 4)THEN ! set itema = nz-1
itema = nz-1
ELSE ! set itema = nz-2
itema = nz-2
END IF
IF( nestgrd /= 1 .OR. mgrid == 1 ) THEN
! For non-nesting case or the
! base-grid of the nested case
DO k=3,itema
DO j=1,ny-1
DO i=1,nx-1
w(i,j,k) = fw(i,j,k)
END DO
END DO
END DO
!call test_dump (pprt,"XXX4solve_pprt",nx,ny,nz,0,0)
IF ( lbcopt == 1 ) THEN ! Internal boundary condition
ebc1=ebc
wbc1=wbc
sbc1=sbc
nbc1=nbc
IF( ebc == 4 ) ebc1=0
IF( wbc == 4 ) wbc1=0
IF( nbc == 4 ) nbc1=0
IF( sbc == 4 ) sbc1=0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mprecv2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mpsend2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
CALL mprecv2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
END IF
CALL acct_interrupt(mp_acct)
CALL lbcw(nx,ny,nz,dtsml1,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, &
ebc1,wbc1,nbc1,sbc1)
CALL acct_stop_inter
ELSE ! External boundary condition
ebc1 = 3
wbc1 = 3
sbc1 = 3
nbc1 = 3
IF( ebc == 0 ) ebc1 = 0
IF( wbc == 0 ) wbc1 = 0
IF( nbc == 0 ) nbc1 = 0
IF( sbc == 0 ) sbc1 = 0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mprecv2dew(w,nx,ny,nz,ebc1,wbc1,3,mptag,tem3)
CALL mpsend2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
CALL mprecv2dns(w,nx,ny,nz,nbc1,sbc1,3,mptag,tem3)
END IF
CALL acct_interrupt(mp_acct)
CALL lbcw(nx,ny,nz,dtsml1,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, &
ebc1,wbc1,nbc1,sbc1)
CALL acct_stop_inter
END IF
ELSE ! For nested interior grid
DO k=3,itema
DO j=2,ny-2
DO i=2,nx-2
w(i,j,k) = fw(i,j,k)
END DO
END DO
END DO
DO k=3,nz-2
DO j=1,ny-1
w( 1,j,k)=w( 1,j,k)+dtsml1 * wdtwb(j,k)
w(nx-1,j,k)=w(nx-1,j,k)+dtsml1 * wdteb(j,k)
END DO
END DO
DO k=3,nz-2
DO i=2,nx-2
w(i, 1,k)=w(i, 1,k)+dtsml1 * wdtsb(i,k)
w(i,ny-1,k)=w(i,ny-1,k)+dtsml1 * wdtnb(i,k)
END DO
END DO
END IF
!call test_dump (pprt,"XXX4Asolve_pprt",nx,ny,nz,0,1)
!
!-----------------------------------------------------------------------
!
! Calculate wcont at time tfuture, including the boundary
! points. Wcont at the lateral boundaries is calculated
! from boundary u, v and w. Wcont at the top and bottom
! depends on the boundary condition option.
!
!-----------------------------------------------------------------------
!
!call test_dump (w,"XXX3sBwc_w1",nx,ny,nz,1,0)
CALL wcontra(nx,ny,nz,u,v,w,mapfct,j1,j2,j3,aj3z, &
rhostr,rhostru,rhostrv,rhostrw,wcont,tem1,tem2)
!
!-----------------------------------------------------------------------
!
! Set the top and bottom boundary conditions for w based on u, v and
! wcont at the top and bottom boundaries.
!
!-----------------------------------------------------------------------
!
CALL acct_interrupt(bc_acct)
CALL vbcw(nx,ny,nz,w,wcont,tbc,bbc,u,v, &
rhostr,rhostru,rhostrv,rhostrw, &
j1,j2,j3)
CALL acct_stop_inter
!
!-----------------------------------------------------------------------
!
! Calculate the new pressure
!
! pprt(new)=pprt(old)+fp-dtsml/j3*tacoef*csndsq*difz(rhobar*w))
!
!-----------------------------------------------------------------------
!
tema = tacoef*dtsml1*dzinv
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
! tem1(i,j,k)= 0.5*(tem10(i,j,k )+tem10(i,j,k-1))
tem1(i,j,k)= 0.5*(rstj3i(i,j,k )+rstj3i(i,j,k-1))
END DO
END DO
END DO
IF( nestgrd == 1 .AND. mgrid /= 1 ) THEN
! For nested interior grid
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-2
pprt(i,j,k) = pprt(i,j,k) + fp(i,j,k) &
- tema*j3inv(i,j,k)*csndsq(i,j,k) &
* (tem1(i,j,k+1)*w(i,j,k+1)-tem1(i,j,k)*w(i,j,k))
END DO
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
pprt( 1,j,k)=pprt( 1,j,k)+dtsml1 * pdtwb(j,k)
pprt(nx-1,j,k)=pprt(nx-1,j,k)+dtsml1 * pdteb(j,k)
END DO
END DO
DO k=2,nz-2
DO i=2,nx-2
pprt(i, 1,k)=pprt(i, 1,k)+dtsml1 * pdtsb(i,k)
pprt(i,ny-1,k)=pprt(i,ny-1,k)+dtsml1 * pdtnb(i,k)
END DO
END DO
!
! Call the pprt bottom boundary condition subroutine to
! compute the hydrostatic pprt at k=1.
!
CALL acct_interrupt(bc_acct)
CALL pprtbbc(nx,ny,nz,g05,buoy2swt,rhostr,pprt,ptprt, &
pbari,ptbari,phydro, &
tem1,tem2) ! tem1 = new pprt at k=1.
CALL acct_stop_inter
! bc's for mp not implemented for nested grids
CALL acct_interrupt(mp_acct)
CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, &
tem1(1,1,1), 0,0,0,0 ,tbc,bbc)
CALL acct_stop_inter
ELSE ! For non-nesting case or the
! base-grid of the nested case
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
pprt(i,j,k) = pprt(i,j,k) + fp(i,j,k) &
- tema*j3inv(i,j,k)*csndsq(i,j,k) &
* (tem1(i,j,k+1)*w(i,j,k+1)-tem1(i,j,k)*w(i,j,k))
END DO
END DO
END DO
!
! Call the pprt bottom boundary condition subroutine to
! compute the hydrostatic pprt at k=1.
!
CALL acct_interrupt(bc_acct)
CALL pprtbbc(nx,ny,nz,g05,buoy2swt,rhostr,pprt,ptprt, &
pbari,ptbari,phydro, &
tem1,tem2) ! tem1 = new pprt at k=1.
CALL acct_stop_inter
!
!-----------------------------------------------------------------------
!
! Apply the boundary conditions for the pressure equation.
!
! For the open boundary case, the boundary value of pprt is
! predicted by the pressure equation.
!
!-----------------------------------------------------------------------
!
!call test_dump (pprt,"XXX5solve_pprt",nx,ny,nz,0,0)
IF( lbcopt == 1 ) THEN ! Internal boundary conditions
ebc1=ebc
wbc1=wbc
sbc1=sbc
nbc1=nbc
IF( ebc == 4 ) ebc1=0
IF( wbc == 4 ) wbc1=0
IF( nbc == 4 ) nbc1=0
IF( sbc == 4 ) sbc1=0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(pprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mprecv2dew(pprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mpsend2dns(pprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
CALL mprecv2dns(pprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
END IF
CALL acct_interrupt(mp_acct)
CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, &
tem1(1,1,1), ebc1,wbc1,nbc1,sbc1,tbc,bbc)
CALL acct_stop_inter
ELSE ! External boundary condition
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(pprt,nx,ny,nz,0,0,0,mptag,tem2)
CALL mprecv2dew(pprt,nx,ny,nz,0,0,0,mptag,tem2)
CALL mpsend2dns(pprt,nx,ny,nz,0,0,0,mptag,tem2)
CALL mprecv2dns(pprt,nx,ny,nz,0,0,0,mptag,tem2)
END IF
CALL acct_interrupt(mp_acct)
CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, &
tem1(1,1,1), 0,0,0,0,tbc,bbc)
CALL exbcp(nx,ny,nz, curtsml, pprt, &
exbcbuf(npr0exb),exbcbuf(nprdtexb))
CALL acct_stop_inter
END IF
END IF
!call test_dump (pprt,"XXX5Asolve_pprt",nx,ny,nz,0,1)
RETURN
END SUBROUTINE solvwpim1
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE TRIDIAG ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE tridiag(nx,ny,nz,a,b,c,d, tem2, tem1, temxy2,w1,w2, & 2,2
nrho,work,trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
temxy1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Solve a tridiagonal linear equation.
!
! The tridiagonal equation set to be solved is
!
! a(k)*w(k-1)+b(k)*w(k)+c(k)*w(k+1)=d(k) for k=3,nz-2
!
! given w(2) and w(nz-1) as the boundary conditions.
!
! The solution w(k) is stored directly into d(k).
!
! Reference: Numerical Recipes: FORTRAN Version, 1989. Page 40.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 9/9/94
!
! MODIFICATION HISTORY:
!
! 10/31/95 (D. Weber)
! Added linear hydrostatic upper radiation condition for w-p.
! References are Klemp and Durran (JAS, 1983) and Chen (1991).
! Includes the use of trigs1,trigs2,ifax1,ifax2.
!
! 07/22/97 (D. Weber)
! Added linear hydrostatic upper radiation condition for w-p using
! an even fft (fftopt=2).
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction
! ny Number of grid points in the y-direction
! nz Number of grid points in the vertical
!
! a lhs coefficient of tridigonal eqations
! b lhs coefficient of tridigonal eqations
! c lhs coefficient of tridigonal eqations
! d rhs coefficient of tridigonal eqations
! temxy2 Pforce array at nz-2.
! trigs1 Array containing pre-computed trig function for
! fftopt=1.
! trigs2 Array containing pre-computed trig function for
! fftopt=1.
! ifax1 Array containing the factors of nx for fftopt=1.
! ifax2 Array containing the factors of ny for fftopt=1.
!
! vwork1 2-D work array for fftopt=2 option.
! vwork2 2-D work array for fftopt=2 option.
! wsave1 Work array for fftopt=2.
! wsave2 Work array for fftopt=2.
!
! OUTPUT:
!
! d The solution w.
!
! WORK ARRAYS:
!
! temxy1 Work array at nz-2.
! tem1 Work array
! tem2 Work array
! work Work array for fft.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
!
! Include files:
!
INCLUDE 'bndry.inc'
INCLUDE 'globcst.inc'
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: a (nx,ny,nz) ! lhs coefficient of tridigonal eqations
REAL :: b (nx,ny,nz) ! lhs coefficient of tridigonal eqations
REAL :: c (nx,ny,nz) ! lhs coefficient of tridigonal eqations
REAL :: d (nx,ny,nz) ! rhs coefficient of tridigonal eqations
! The contains solution w on exit.
REAL :: temxy2(nx,ny) ! Pforce array at nz-2.
REAL :: trigs1(3*(nx-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
REAL :: trigs2(3*(ny-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
INTEGER :: ifax1(13) ! Array containing the factors of nx for
! fftopt=1.
INTEGER :: ifax2(13) ! Array containing the factors of ny for
! fftopt=1.
REAL :: vwork1 (nx+1,ny+1) ! 2-D work array for fftopt=2.
REAL :: vwork2 (ny,nx+1) ! 2-D work array for fftopt=2.
REAL :: wsave1 (3*(ny-1)+15) ! Work array for fftopt=2.
REAL :: wsave2 (3*(nx-1)+15) ! Work array for fftopt=2.
REAL :: temxy1(nx+1,ny+1) ! Work array at nz-2.
REAL :: tem1 (nx,ny) ! 2-D work array.
REAL :: tem2 (nx,ny,nz) ! 3-D work array.
REAL :: work (ny,nx) ! 2-D work array for fft code.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
REAL :: nrho
REAL :: w1,w2,wc
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO j=1,ny-1
DO i=1,nx-1
d(i,j,3)=d(i,j,3)/b(i,j,3)
tem1(i,j)=b(i,j,3)
END DO
END DO
DO k=4,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k) = c(i,j,k-1)/tem1(i,j)
tem1(i,j)=b(i,j,k)-a(i,j,k)*tem2(i,j,k)
d(i,j,k)=(d(i,j,k)-a(i,j,k)*d(i,j,k-1))/tem1(i,j)
END DO
END DO
END DO
IF(tbc == 4)THEN ! apply upper radiation boundary condition.
!
! Computing the new c(nz-2) (which is stored in nz-1)....
!
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,nz-1) = c(i,j,nz-2)/tem1(i,j)
END DO
END DO
wc = tem2(nx/2,ny/2,nz-1) ! coef. for w(nz-1) in tri. elim. phase
!
! Combining the pforce array and d(nz-2) from the elimination phase.
! Note: these arrays are a function of x and/or y. All others variables
! are determined from base state which are not a function
! of x,y. (note: zflat must be at or below the height of scalar nz-2)
!
DO j=1,ny-1
DO i=1,nx-1
temxy1(i,j) = temxy2(i,j) + w2*d(i,j,nz-2)
END DO
END DO
!-----------------------------------------------------------------------
!
! Call the upper radiation subroutine to update w at the top (nz-2)
!
!-----------------------------------------------------------------------
IF(fftopt == 2)THEN ! call the even vfftpack fft...
CALL uprad3(nx,ny,w1,w2,wc,nrho,wsave1,wsave2,vwork1,vwork2, &
temxy1,work)
ELSE IF(fftopt == 1)THEN ! call the periodic fft....
CALL uprad1(nx,ny,w1,w2,wc,nrho,ifax1,ifax2,trigs1,trigs2, &
temxy1,work)
END IF ! end of fftopt....
DO j=1,ny-1
DO i=1,nx-1
d(i,j,nz-1) = temxy1(i,j)
END DO
END DO
!-----------------------------------------------------------------------
!
! Perform the back substitution phase of the tridiagonal solver.
!
!-----------------------------------------------------------------------
DO k=nz-2,3,-1
DO j=1,ny-1
DO i=1,nx-1
d(i,j,k)=d(i,j,k)-tem2(i,j,k+1)*d(i,j,k+1)
END DO
END DO
END DO
ELSE IF(tbc /= 4)THEN ! perform backsubtitution for other bc's.
DO k=nz-3,3,-1
DO j=1,ny-1
DO i=1,nx-1
d(i,j,k)=d(i,j,k)-tem2(i,j,k+1)*d(i,j,k+1)
END DO
END DO
END DO
END IF ! end of tbc loop...
RETURN
END SUBROUTINE tridiag
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SOLVPT_LRG ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE solvpt_lrg(nx,ny,nz, exbcbufsz, dtbig1, & 1,20
ptprt, ptdteb,ptdtwb,ptdtnb,ptdtsb, &
rhostr,rhostri,ptforce, exbcbuf, tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinate the time integration of the potential temperature
! equation in large time steps.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/91
!
! MODIFICATION HISTORY:
!
! 9/17/94 (M. Xue)
! Rewritten for small time step integration of ptprt equation.
!
! 11/6/1995 (M. Xue)
! Added option for fourth order vertical advection for ptbar.
!
! 4/17/96 (M. Xue)
! Removed the block for 4th order advection of ptbar.
!
! 4/24/1997 (M. Xue)
! Rewrote as SOLVPT_LRG based on SOLVPT.
!
! 10/5/1998 (Dan Weber)
! Added rhostri for efficiency.
!
!-----------------------------------------------------------------------
!
! 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 for this call (s)
!
! ptprt Perturbation potential temperature at times tpast and
! tpresent (K)
!
! ptdteb Time tendency of the ptprt field at the east boundary
! ptdtwb Time tendency of the ptprt field at the west boundary
! ptdtnb Time tendency of the ptprt field at the north boundary
! ptdtsb Time tendency of the ptprt field at the south boundary
!
! rhostr Base state density rhobar times j3 (kg/m**3)
! rhostri Inverse base state density rhobar times j3 (kg/m**3)
!
! ptforce Gravity wave inactive forcing terms in pt-eq.
! (K*kg/(m**3*s))
!
! OUTPUT:
!
! ptprt Perturbation potential temperature at time tfuture (K)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INCLUDE 'timelvls.inc'
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: dtbig1 ! The big time step size for this call
! (s)
REAL :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: ptdteb(ny,nz) ! Time tendency of ptprt at east
! boundary
REAL :: ptdtwb(ny,nz) ! Time tendency of ptprt at west
! boundary
REAL :: ptdtnb(nx,nz) ! Time tendency of ptprt at north
! boundary
REAL :: ptdtsb(nx,nz) ! Time tendency of ptprt at south
! boundary
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: rhostri(nx,ny,nz) ! Inverse base state density rhobar times j3.
REAL :: ptforce(nx,ny,nz) ! Gravity wave inactive forcing terms
! in potential temperature eq.
! (K*kg/(m**3*s))
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: tem1(nx,ny,nz) ! Work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k, tstrtlvl
INTEGER :: ebc1,wbc1,nbc1,sbc1
REAL :: deltat
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Integrate forward by one timestep the potential temperature
! equation. When PT-eq. is integrated inside small time steps,
! it is stepped forward by a small time step, otherwise, it is
! stepped forward by a large time step using leapfrog scheme
! (i.e. 2*dtbig1 from ptprt at tpast).
!
!-----------------------------------------------------------------------
!
IF(sadvopt == 4) THEN ! Forward-based FCT scheme
deltat = dtbig1
tstrtlvl = tpresent
ELSE
deltat = dtbig1*2
tstrtlvl = tpast
END IF
IF( nestgrd == 1 .AND. mgrid /= 1 ) THEN
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-2
ptprt(i,j,k,tfuture)=ptprt(i,j,k,tstrtlvl) &
+deltat*ptforce(i,j,k)*rhostri(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Integrate PT equation and the boundary conditions for a nested grid.
!
!-----------------------------------------------------------------------
!
DO k=2,nz-2
DO i=1,nx-1
ptprt(i, 1,k,tfuture)= &
2*ptprt(i, 1,k,tpresent)-ptprt(i, 1,k,tpast)
ptprt(i,ny-1,k,tfuture)= &
2*ptprt(i,ny-1,k,tpresent)-ptprt(i,ny-1,k,tpast)
END DO
END DO
DO k=2,nz-2
DO j=2,ny-2
ptprt( 1,j,k,tfuture)= &
2*ptprt( 1,j,k,tpresent)-ptprt( 1,j,k,tpast)
ptprt(nx-1,j,k,tfuture)= &
2*ptprt(nx-1,j,k,tpresent)-ptprt(nx-1,j,k,tpast)
END DO
END DO
! bc's for mp not implemented for nested grids
CALL acct_interrupt(mp_acct)
CALL bcsclr(nx,ny,nz, deltat*0.5 , &
ptprt(1,1,1,tstrtlvl),ptprt(1,1,1,tpresent), &
ptprt(1,1,1,tfuture),ptdteb,ptdtwb,ptdtnb,ptdtsb, &
0,0,0,0 ,tbc,bbc)
CALL acct_stop_inter
ELSE
!
!-----------------------------------------------------------------------
!
! Integrate PT equation and the boundary conditions for a base grid.
!
!-----------------------------------------------------------------------
!
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
ptprt(i,j,k,tfuture)=ptprt(i,j,k,tstrtlvl) &
+deltat*ptforce(i,j,k)*rhostri(i,j,k)
END DO
END DO
END DO
!call test_dump (ptprt,"XXXsolve_ptprt",nx,ny,nz,0,0)
IF ( lbcopt == 1 ) THEN ! Internal boundary conditions
ebc1=ebc
wbc1=wbc
sbc1=sbc
nbc1=nbc
IF( ebc == 4 ) ebc1=0
IF( wbc == 4 ) wbc1=0
IF( nbc == 4 ) nbc1=0
IF( sbc == 4 ) sbc1=0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(ptprt(1,1,1,tfuture),nx,ny,nz,ebc1,wbc1,0,mptag,tem1)
CALL mprecv2dew(ptprt(1,1,1,tfuture),nx,ny,nz,ebc1,wbc1,0,mptag,tem1)
CALL mpsend2dns(ptprt(1,1,1,tfuture),nx,ny,nz,nbc1,sbc1,0,mptag,tem1)
CALL mprecv2dns(ptprt(1,1,1,tfuture),nx,ny,nz,nbc1,sbc1,0,mptag,tem1)
END IF
CALL acct_interrupt(mp_acct)
CALL bcsclr(nx,ny,nz, deltat*0.5 , &
ptprt(1,1,1,tstrtlvl),ptprt(1,1,1,tpresent), &
ptprt(1,1,1,tfuture),ptdteb,ptdtwb,ptdtnb,ptdtsb, &
ebc1,wbc1,nbc1,sbc1,tbc,bbc)
CALL acct_stop_inter
ELSE ! External boundary condition
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(ptprt(1,1,1,tfuture),nx,ny,nz,0,0,0,mptag,tem1)
CALL mprecv2dew(ptprt(1,1,1,tfuture),nx,ny,nz,0,0,0,mptag,tem1)
CALL mpsend2dns(ptprt(1,1,1,tfuture),nx,ny,nz,0,0,0,mptag,tem1)
CALL mprecv2dns(ptprt(1,1,1,tfuture),nx,ny,nz,0,0,0,mptag,tem1)
END IF
CALL acct_interrupt(mp_acct)
CALL bcsclr(nx,ny,nz, deltat*0.5 , &
ptprt(1,1,1,tstrtlvl),ptprt(1,1,1,tpresent), &
ptprt(1,1,1,tfuture),ptdteb,ptdtwb,ptdtnb,ptdtsb, &
0,0,0,0,tbc,bbc)
CALL exbcpt(nx,ny,nz, curtim+dtbig, ptprt(1,1,1,tfuture), &
exbcbuf(npt0exb),exbcbuf(nptdtexb))
CALL acct_stop_inter
END IF
END IF
!call test_dump (ptprt,"XXXAsolve_ptprt",nx,ny,nz,0,1)
RETURN
END SUBROUTINE solvpt_lrg
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SOLVPT_SML ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE solvpt_sml(nx,ny,nz, exbcbufsz, dtbig1,dtsml1,curtsml, & 1,21
ptprt,w, ptdteb,ptdtwb,ptdtnb,ptdtsb, &
ptbar,rhostr,rhostri,rhostrw,j3,aj3z,ptforce, &
exbcbuf, &
ptadv,tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinate the time integration of the potential temperature
! equation in small time steps.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/91
!
! MODIFICATION HISTORY:
!
! 9/17/94 (M. Xue)
! Rewritten for small time step integration of ptprt equation.
!
! 11/6/1995 (M. Xue)
! Added option for fourth order vertical advection for ptbar.
!
! 4/17/96 (M. Xue)
! Removed the block for 4th order advection of ptbar.
!
! 4/24/1997 (M. Xue)
! Rewrote as SOLVPT_SML based on SOLVPT.
!
! 9/18/1998 (D. Weber)
! Added array aj3z and rhostri,w for efficiency.
!
!-----------------------------------------------------------------------
!
! 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 for this call (s)
! dtsml1 The big time step size for this call (s)
! curtsml Current time during small time step integration.
!
! ptprt Perturbation potential temperature at times tpast and
! tpresent (K)
! w Vertical component of Cartesian velocity
!
! ptdteb Time tendency of the ptprt field at the east boundary
! ptdtwb Time tendency of the ptprt field at the west boundary
! ptdtnb Time tendency of the ptprt field at the north boundary
! ptdtsb Time tendency of the ptprt field at the south boundary
!
! ptbar Base state potential temperature (K)
! rhostr Base state density rhobar times j3 (kg/m**3)
! rhostri Inverse base state density rhobar times j3 (kg/m**3)
! rhostrw rhostr averaged to w-point.
!
! j3 Coordinate transformation Jacobian d(zp)/dz
! aj3z Avgz of the coordinate transformation Jacobian d(zp)/dz
!
! ptforce Gravity wave inactive forcing terms in pt-eq.
! (K*kg/(m**3*s))
!
! OUTPUT:
!
! ptprt Perturbation potential temperature at time tfuture (K)
!
! WORK ARRAYS:
!
! ptadv Advection of base state potential temperature
! tem1 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 ! The big time step size for this call
! (s)
REAL :: dtsml1 ! The big time step size for this call
! (s)
REAL :: curtsml ! Current time during small time step
! integration.
REAL :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: w (nx,ny,nz,nt) ! Total w-velocity (m/s)
REAL :: ptdteb(ny,nz) ! Time tendency of ptprt at east
! boundary
REAL :: ptdtwb(ny,nz) ! Time tendency of ptprt at west
! boundary
REAL :: ptdtnb(nx,nz) ! Time tendency of ptprt at north
! boundary
REAL :: ptdtsb(nx,nz) ! Time tendency of ptprt at south
! boundary
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) ! Inverse base state density rhobar times j3.
REAL :: rhostrw(nx,ny,nz) ! rhostr averaged to w-point.
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian
! defined as d( zp )/d( z ).
REAL :: aj3z (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
REAL :: ptforce(nx,ny,nz) ! Gravity wave inactive forcing terms
! in potential temperature eq.
! (K*kg/(m**3*s))
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: ptadv (nx,ny,nz) ! Temporary array to store base state
! potential temperature advection.
REAL :: tem1 (nx,ny,nz) ! Temporary work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k, onvf
INTEGER :: ebc1,wbc1,nbc1,sbc1
REAL :: deltat, dz05 , time
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Integrate forward by one timestep the potential temperature
! equation. When PT-eq. is integrated inside small time steps,
! it is stepped forward by a small time step, otherwise, it is
! stepped forward by a large time step using leapfrog scheme
! (i.e. 2*dtbig1 from ptprt at tpast).
!
!-----------------------------------------------------------------------
!
deltat = dtsml1
time = curtsml
!
!-----------------------------------------------------------------------
!
! Base state potential temperature advection. This term is added to
! the array ptadv to yield the total potential temperature advection.
!
! ptbar is assumed to be independent of physical x and y, therefore
! d(ptbar)/dx and d(ptbar)/dy for constant z are zero, the base state
! advection is -w*d(ptbar)/dzp = -w/j3*d(ptbar)/dz.
!
! This term is responsible for internal gravity waves, when potential
! temperature equation is solved inside the small time steps, this
! term is evaluated inside the small time steps.
!
!-----------------------------------------------------------------------
!
! dz05 = dz*0.5
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=rhostrw(i,j,k)*w(i,j,k,tfuture) &
*(ptbar(i,j,k)-ptbar(i,j,k-1)) &
*dzinv/aj3z(i,j,k)
END DO
END DO
END DO
onvf = 0
CALL avgz(tem1 , onvf, &
nx,ny,nz, 1,nx-1, 1,ny-1, 2,nz-2, ptadv)
IF( nestgrd == 1 .AND. mgrid /= 1 ) THEN
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-2
ptprt(i,j,k,tfuture)=ptprt(i,j,k,tfuture) &
+deltat*(ptforce(i,j,k)-ptadv(i,j,k))*rhostri(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Integrate PT equation and the boundary conditions for a nested grid.
!
!-----------------------------------------------------------------------
!
DO k=2,nz-2
DO j=1,ny-1
ptprt( 1,j,k,tfuture)=ptprt( 1,j,k,tfuture) &
+deltat* ptdtwb(j,k)
ptprt(nx-1,j,k,tfuture)=ptprt(nx-1,j,k,tfuture) &
+deltat* ptdteb(j,k)
END DO
END DO
DO k=2,nz-2
DO i=2,nx-2
ptprt(i, 1,k,tfuture)=ptprt(i, 1,k,tfuture) &
+deltat* ptdtsb(i,k)
ptprt(i,ny-1,k,tfuture)=ptprt(i,ny-1,k,tfuture) &
+deltat* ptdtnb(i,k)
END DO
END DO
! bc's for mp not implemented for nested grids
CALL acct_interrupt(mp_acct)
CALL bcsclr(nx,ny,nz, deltat*0.5 , &
ptprt(1,1,1,tfuture),ptprt(1,1,1,tpresent), &
ptprt(1,1,1,tfuture),ptdteb,ptdtwb,ptdtnb,ptdtsb, &
0,0,0,0 ,tbc,bbc)
CALL acct_stop_inter
ELSE
!
!-----------------------------------------------------------------------
!
! Integrate PT equation and the boundary conditions for a base grid.
!
!-----------------------------------------------------------------------
!
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
ptprt(i,j,k,tfuture)=ptprt(i,j,k,tfuture) &
+deltat*(ptforce(i,j,k)-ptadv(i,j,k))*rhostri(i,j,k)
END DO
END DO
END DO
!call test_dump (ptprt,"XXX1solve_ptprt",nx,ny,nz,0,0)
IF ( lbcopt == 1 ) THEN ! Internal boundary conditions
ebc1=ebc
wbc1=wbc
sbc1=sbc
nbc1=nbc
IF( ebc == 4 ) ebc1=0
IF( wbc == 4 ) wbc1=0
IF( nbc == 4 ) nbc1=0
IF( sbc == 4 ) sbc1=0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(ptprt(1,1,1,tfuture),nx,ny,nz,ebc1,wbc1,0,mptag,tem1)
CALL mprecv2dew(ptprt(1,1,1,tfuture),nx,ny,nz,ebc1,wbc1,0,mptag,tem1)
CALL mpsend2dns(ptprt(1,1,1,tfuture),nx,ny,nz,nbc1,sbc1,0,mptag,tem1)
CALL mprecv2dns(ptprt(1,1,1,tfuture),nx,ny,nz,nbc1,sbc1,0,mptag,tem1)
END IF
CALL acct_interrupt(mp_acct)
CALL bcsclr(nx,ny,nz, deltat*0.5 , &
ptprt(1,1,1,tfuture),ptprt(1,1,1,tpresent), &
ptprt(1,1,1,tfuture),ptdteb,ptdtwb,ptdtnb,ptdtsb, &
ebc1,wbc1,nbc1,sbc1,tbc,bbc)
CALL acct_stop_inter
ELSE ! External boundary condition
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(ptprt(1,1,1,tfuture),nx,ny,nz,0,0,0,mptag,tem1)
CALL mprecv2dew(ptprt(1,1,1,tfuture),nx,ny,nz,0,0,0,mptag,tem1)
CALL mpsend2dns(ptprt(1,1,1,tfuture),nx,ny,nz,0,0,0,mptag,tem1)
CALL mprecv2dns(ptprt(1,1,1,tfuture),nx,ny,nz,0,0,0,mptag,tem1)
END IF
CALL acct_interrupt(mp_acct)
CALL bcsclr(nx,ny,nz, deltat*0.5 , &
ptprt(1,1,1,tfuture),ptprt(1,1,1,tpresent), &
ptprt(1,1,1,tfuture),ptdteb,ptdtwb,ptdtnb,ptdtsb, &
0,0,0,0,tbc,bbc)
CALL exbcpt(nx,ny,nz, time , ptprt(1,1,1,tfuture), &
exbcbuf(npt0exb),exbcbuf(nptdtexb))
CALL acct_stop_inter
END IF
END IF
!call test_dump (ptprt,"XXX1Asolve_ptprt",nx,ny,nz,0,1)
RETURN
END SUBROUTINE solvpt_sml
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SOLVQV ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE solvqv(nx,ny,nz, exbcbufsz, dtbig1, & 1,32
qv,u,v,wcont, ustr,vstr,wstr, &
qvdteb,qvdtwb,qvdtnb,qvdtsb, &
rhostr,rhostri,qvbar,kmh,kmv,rprntl,qvsflx,pbldpth, &
x,y,z,zp,mapfct,j1,j2,j3,aj3x,aj3y,j3inv,qvcumsrc, &
usflx,vsflx,ptsflx,ptsfc,qvsfc,ptbar,ptprt, &
exbcbuf,bcrlx, &
qadv,qmix, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem1_0,tem2_0,tem3_0)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinate the time integration of the equation for water vapor
! specific humidity qv, and the setting of boundary conditions.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/91
!
! MODIFICATION HISTORY:
!
! 5/20/92 (M. Xue)
! Added full documentation.
!
! 2/10/93 (K. Droegemeier)
! Cleaned up documentation.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 8/22/95 (M. Xue)
! Added ptcumsrc term to the right hand side of qv equation.
! It was omitted before.
!
! 08/30/95 (M. Xue and Y. Liu)
! Changed EXBC for water variables to zero-gradient when the
! variables are missing in the boundary file.
!
! 08/30/95 (Yuhe Liu)
! Fixed a bug which double computes the qv advection and mixing.
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 4/1/96 (Donghai Wang, X. Song and M. Xue)
! Added the implicit treatment for the vertical mixing.
!
! 3/24/1997 (M. Xue)
! Code to handle the case of forward time integration added.
!
! 7/10/1997 (Fanyou Kong -- CMRP)
! Added the positive definite advection scheme option (sadvopt = 5)
!
! 7/17/1998 (M. Xue)
! Changed call to ADVQFCT.
!
! 9/18/1998 (D. Weber)
! Added arrays aj3x,y, u,v,wstr,rhostri do not alter!
!
!-----------------------------------------------------------------------
!
! 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)
!
! qv Water vapor specific humidity at tpast and tpresent
! (kg/kg)
!
! u x component of velocity at all time levels (m/s)
! v y component of velocity at all time levels (m/s)
! wcont Contravariant vertical velocity (m/s)
!
! ustr Appropriate ustr for advection.
! vstr Appropriate vstr for advection.
! wstr Appropriate wstr for advection.
!
! qvdteb Time tendency of the qv field at the east boundary
! qvdtwb Time tendency of the qv field at the west boundary
! qvdtnb Time tendency of the qv field at the north boundary
! qvdtsb Time tendency of the qv field at the south boundary
!
! rhostr Base state density rhobar times j3 (kg/m**3)
! rhostri Inverse base state density rhobar times j3 (kg/m**3)
! qvbar Base state water vapor specific humidity (kg/kg)
!
! 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
!
! qvsflx Surface flux of moisture (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
! qvcumsrc Source term in qv-equation due to cumulus parameterization
!
! OUTPUT:
!
! qv Water vapor specific humidity at time tfuture (kg/kg)
!
! WORK ARRAYS:
!
! qadv Advection term of water vapor eq. (kg/(m**3*s)).
! A local array.
! qmix Total mixing in water vapor eq. (kg/(m**3*s)).
! A local array.
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
! tem4 Temporary work array.
! tem5 Temporary work array.
! tem6 Temporary work array.
! tem6 Temporary work array.
! tem7 Temporary work array.
! tem8 Temporary work array.
!
! tem1_0 Temporary work array.
! tem2_0 Temporary work array.
! tem3_0 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
REAL :: qv (nx,ny,nz,nt) ! Water vapor specific humidity (kg/kg)
REAL :: u (nx,ny,nz,nt) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: ustr (nx,ny,nz) ! Appropriate ustr for advection.
REAL :: vstr (nx,ny,nz) ! Appropriate vstr for advection.
REAL :: wstr (nx,ny,nz) ! Appropriate wstr for advection.
REAL :: qvdteb(ny,nz) ! Time tendency of qv at east boundary
REAL :: qvdtwb(ny,nz) ! Time tendency of qv at west boundary
REAL :: qvdtnb(nx,nz) ! Time tendency of qv at north boundary
REAL :: qvdtsb(nx,nz) ! Time tendency of qv at south boundary
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: rhostri(nx,ny,nz) ! Inverse base state density rhobar times j3.
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific
! humidity (kg/kg)
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
REAL :: vsflx(nx,ny) ! Surface flux of v-momentum
REAL :: ptsflx(nx,ny) ! Surface heat flux (K*kg/(m*s**2))
REAL :: qvsflx(nx,ny) ! Surface flux of moisture (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 :: qvcumsrc(nx,ny,nz) ! Source in qv-equation due to cumulus
! parameterization
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature
REAL :: ptsfc(nx,ny) ! Potential temperature at the ground level (K)
REAL :: qvsfc(nx,ny) ! Effective qv at the surface (kg/kg)
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: bcrlx (nx,ny) ! EXBC relaxation coefficients
REAL :: qadv(nx,ny,nz) ! Advection of water vapor (kg/(m**3*s))
! A local array.
REAL :: qmix(nx,ny,nz) ! Total mixing of water vapor
! (kg/(m**3*s))
! A local array.
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 :: 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.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
INTEGER :: tstrtlvl ! Starting level of time integration
REAL :: deltat
INTEGER :: ebc1,wbc1,nbc1,sbc1, qflag
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Compute the advection term of the water vapor specific humidity
! equation and store the result in array qadv.
!
!-----------------------------------------------------------------------
!
CALL set_acct(advs_acct)
IF(sadvopt == 4) THEN ! Forward-based FCT scheme
deltat = dtbig1
tstrtlvl = tpresent
ELSE
deltat = dtbig1*2
tstrtlvl = tpast
END IF
IF (sadvopt == 1 .OR. sadvopt == 2 .OR. sadvopt == 3 ) THEN
! 2nd or 4th order advection
CALL advq(nx,ny,nz,qv,u,v,wcont,ustr,vstr,wstr, &
rhostr,qvbar, mapfct, &
qadv, &
tem1,tem2,tem3,tem4)
ELSE IF( sadvopt == 4.OR.sadvopt == 5) THEN ! FCT advection
CALL advqfct(nx,ny,nz,dtbig1,qv,u,v,wcont, ustr,vstr,wstr, &
rhostr,rhostri,mapfct,j3,qadv, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem1_0,tem2_0,tem3_0)
END IF
CALL set_acct(cmix_acct)
!
!-----------------------------------------------------------------------
!
! Compute the mixing terms for the water vapor specific humidity
! equation. This includes both physical and computational
! mixing. Store the result in array qmix.
!
!-----------------------------------------------------------------------
!
CALL mixqv(nx,ny,nz, &
qv(1,1,1,tstrtlvl),rhostr,qvbar,kmh,kmv,rprntl,qvsflx, &
pbldpth(1,1,tpresent), &
x,y,z,zp, mapfct, j1,j2,j3,aj3x,aj3y,j3inv, &
usflx,vsflx,ptsflx,ptsfc,qvsfc,ptbar,ptprt, &
qmix, &
tem1,tem2,tem3,tem4,tem5,tem6)
!
!-----------------------------------------------------------------------
!
! Call BRLXQ to calculate the boundary relaxation and computation
! mixing for qv
!
! qmix = qmix + qvexbc_term
!
!-----------------------------------------------------------------------
!
IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN
CALL acct_interrupt(bc_acct)
qflag = 1
CALL brlxq(nx,ny,nz, deltat*0.5, qflag, qv,rhostr, qmix, &
exbcbuf(nqv0exb), exbcbuf(nqc0exb), &
exbcbuf(nqr0exb), exbcbuf(nqi0exb), &
exbcbuf(nqs0exb), exbcbuf(nqh0exb), &
exbcbuf(nqvdtexb),exbcbuf(nqcdtexb), &
exbcbuf(nqrdtexb),exbcbuf(nqidtexb), &
exbcbuf(nqsdtexb),exbcbuf(nqhdtexb),bcrlx, &
tem1,tem2,tem3,tem4)
CALL acct_stop_inter
END IF
!
!-----------------------------------------------------------------------
!
! Calculate qvforce, store the result in tem1.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=-qadv(i,j,k)+rhostr(i,j,k)*qvcumsrc(i,j,k) &
+qmix(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Treat the vertically implicit mixing term.
!
!-----------------------------------------------------------------------
!
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
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,qv(1,1,1,tstrtlvl),rhostr, &
tem2,tem1,tem3,tem4,tem5,tem6)
END IF
!
!-----------------------------------------------------------------------
!
! Integrate forward by one timestep the water vapor specific humidity
! equation, yielding qv at time = tfuture.
!
!-----------------------------------------------------------------------
!
IF( nestgrd == 1 .AND. mgrid /= 1 ) THEN
CALL set_acct(tinteg_acct)
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-2
qv(i,j,k,tfuture)=qv(i,j,k,tstrtlvl) &
+deltat*tem1(i,j,k)*rhostri(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set the boundary conditions on qv for an adaptive (nested)
! grid run. If using only one grid, this portion of the code is
! skipped....proceed to next comment block.
!
!-----------------------------------------------------------------------
!
DO k=2,nz-2
DO i=1,nx-1
qv(i, 1,k,tfuture)= &
2*qv(i, 1,k,tpresent)-qv(i, 1,k,tpast)
qv(i,ny-1,k,tfuture)= &
2*qv(i,ny-1,k,tpresent)-qv(i,ny-1,k,tpast)
END DO
END DO
DO k=2,nz-2
DO j=2,ny-2
qv( 1,j,k,tfuture)= &
2*qv( 1,j,k,tpresent)-qv( 1,j,k,tpast)
qv(nx-1,j,k,tfuture)= &
2*qv(nx-1,j,k,tpresent)-qv(nx-1,j,k,tpast)
END DO
END DO
! bc's for mp not implemented for nested grids
CALL acct_interrupt(bc_acct)
CALL bcqv(nx,ny,nz, deltat*0.5, &
qv,qvbar, qvdteb,qvdtwb,qvdtnb,qvdtsb, &
0,0,0,0,tbc,bbc)
CALL acct_stop_inter
ELSE
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
qv(i,j,k,tfuture)=qv(i,j,k,tstrtlvl) &
+deltat*tem1(i,j,k)*rhostri(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set the boundary conditions on qv for a NON-adaptive (uniform)
! grid run.
!
!-----------------------------------------------------------------------
!
CALL set_acct(bc_acct)
!call test_dump (qv,"XXXsolve_qv",nx,ny,nz,0,0)
IF ( lbcopt == 1 ) THEN
ebc1=ebc
wbc1=wbc
sbc1=sbc
nbc1=nbc
IF( ebc == 4 ) ebc1=0
IF( wbc == 4 ) wbc1=0
IF( nbc == 4 ) nbc1=0
IF( sbc == 4 ) sbc1=0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(qv(1,1,1,tfuture),nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mprecv2dew(qv(1,1,1,tfuture),nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mpsend2dns(qv(1,1,1,tfuture),nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
CALL mprecv2dns(qv(1,1,1,tfuture),nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
END IF
CALL acct_interrupt(mp_acct)
CALL bcqv(nx,ny,nz, dtbig1, &
qv,qvbar, qvdteb,qvdtwb,qvdtnb,qvdtsb, &
ebc1,wbc1,nbc1,sbc1,tbc,bbc)
CALL acct_stop_inter
ELSE
ebc1 = 3
wbc1 = 3
sbc1 = 3
nbc1 = 3
IF( ebc == 0 ) ebc1 = 0
IF( wbc == 0 ) wbc1 = 0
IF( nbc == 0 ) nbc1 = 0
IF( sbc == 0 ) sbc1 = 0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(qv(1,1,1,tfuture),nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mprecv2dew(qv(1,1,1,tfuture),nx,ny,nz,ebc1,wbc1,0,mptag,tem2)
CALL mpsend2dns(qv(1,1,1,tfuture),nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
CALL mprecv2dns(qv(1,1,1,tfuture),nx,ny,nz,nbc1,sbc1,0,mptag,tem2)
END IF
CALL acct_interrupt(mp_acct)
CALL bcqv(nx,ny,nz, dtbig1, &
qv,qvbar, qvdteb,qvdtwb,qvdtnb,qvdtsb, &
ebc1,wbc1,nbc1,sbc1,tbc,bbc)
qflag = 1
CALL exbcq(nx,ny,nz, qflag, curtim+dtbig,qv(1,1,1,tfuture), &
exbcbuf(nqv0exb), exbcbuf(nqvdtexb))
CALL acct_stop_inter
END IF
END IF
!call test_dump (qv,"XXXAsolve_qv",nx,ny,nz,0,1)
RETURN
END SUBROUTINE solvqv
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SOLVQ ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag, & 5,31
q,u,v,wcont, ustr,vstr,wstr, &
qdteb,qdtwb,qdtnb,qdtsb, &
rhostr,rhostri,kmh,kmv,rprntl, &
x,y,z,zp, mapfct, j1,j2,j3,aj3x,aj3y,j3inv, &
qccumsrc,qrcumsrc,qicumsrc,qscumsrc, &
exbcbuf,bcrlx, &
qadv,qmix, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem1_0,tem2_0,tem3_0)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinate the time integration of the equations for the water
! substance quantities qc,qr,qi,qs and qh, i.e. all water variables
! except the water vapor specific humidity qv.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Ming Xue
! 10/10/91
!
! MODIFICATION HISTORY:
!
! 5/20/92 (K. Droegemeier and M. Xue)
! Added full documentation.
!
! 2/10/93 (K. Droegemeier)
! Cleaned up documentation.
!
! 8/12/95 (M. Xue)
! Added flag qflag.
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 4/1/96 (Donghai Wang, X. Song and M. Xue)
! Added the implicit treatment for the vertical mixing.
!
! 3/24/1997 (M. Xue)
! Code to handle the case of forward time integration added.
!
! 7/10/1997 (Fanyou Kong -- CMRP)
! Added the positive definite advection scheme option (sadvopt = 5)
!
! 4/15/1998 (Donghai Wang)
! Added the source terms to the right hand terms of the qc,qr,qi,qs
! equations due to the K-F cumulus parameterization.
!
! 7/17/1998 (M. Xue)
! Changed call to ADVQFCT.
!
! 9/18/1998 (D. Weber)
! Added arrays aj3x,y, u,v,wstr,rhostri, do not alter!
!
!-----------------------------------------------------------------------
!
! 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)
! qflag Indicator for water/ice species
!
! q One of the liquid or ice variables at tpast and
! tpresent (kg/kg)
!
! u x component of velocity at all time levels (m/s)
! v y component of velocity at all time levels (m/s)
! wcont Contravariant vertical velocity (m/s)
! ustr Appropriate ustr for advection.
! vstr Appropriate vstr for advection.
! wstr Appropriate wstr for advection.
!
! qdteb Time tendency of liquid/ice variable q at east boundary
! qdtwb Time tendency of liquid/ice variable q at west boundary
! qdtnb Time tendency of liquid/ice variable q at north
! boundary
! qdtsb Time tendency of liquid/ice variable q at south
! boundary
!
! rhostr Base state density rhobar times j3 (kg/m**3)
! rhostri Inverse 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
!
! 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
! qccumsrc Source term in qc-equation due to cumulus parameterization
! qrcumsrc Source term in qr-equation due to cumulus parameterization
! qicumsrc Source term in qi-equation due to cumulus parameterization
! qscumsrc Source term in qs-equation due to cumulus parameterization
!
!
! OUTPUT:
!
! q Liquid/ice variable q at tfuture (kg/kg)
!
! WORK ARRAYS:
!
! qadv Advection term of water variable eq. (kg/(m**3*s)).
! A local array.
! qmix Total mixing in water variable eq. (kg/(m**3*s)).
! A local array.
! 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.
!
! tem1_0 Temporary work array.
! tem2_0 Temporary work array.
! tem3_0 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
INTEGER :: qflag ! Indicator for water/ice species
REAL :: q (nx,ny,nz,nt) ! One of the water/ice variables (kg/kg)
REAL :: u (nx,ny,nz,nt) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: ustr (nx,ny,nz) ! Appropriate ustr for advection (m/s)
REAL :: vstr (nx,ny,nz) ! Appropriate vstr for advection (m/s)
REAL :: wstr (nx,ny,nz) ! Appropriate wstr for advection (m/s)
REAL :: qdteb(ny,nz) ! Time tendency of liquid/ice at
! e-boundary
REAL :: qdtwb(ny,nz) ! Time tendency of liquid/ice at
! w-boundary
REAL :: qdtnb(nx,nz) ! Time tendency of liquid/ice at
! n-boundary
REAL :: qdtsb(nx,nz) ! Time tendency of liquid/ice at
! s-boundary
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: rhostri(nx,ny,nz) ! Inverse 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 :: 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 :: qccumsrc(nx,ny,nz) ! Source in qc-equation due to cumulus
! parameterization
REAL :: qrcumsrc(nx,ny,nz) ! Source in qr-equation due to cumulus
! parameterization
REAL :: qicumsrc(nx,ny,nz) ! Source in qi-equation due to cumulus
! parameterization
REAL :: qscumsrc(nx,ny,nz) ! Source in qs-equation due to cumulus
! parameterization
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: bcrlx (nx,ny) ! EXBC relaxation coefficients
REAL :: qadv(nx,ny,nz) ! Advection of water/ice substance
! (kg/(m**3*s)). A local array.
REAL :: qmix(nx,ny,nz) ! Total mixing of water/ice substance
! (kg/(m**3*s)). A local array.
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 :: tem1_0(0:nx,0:ny,0:nz) ! automatic work array
REAL :: tem2_0(0:nx,0:ny,0:nz) ! automatic work array
REAL :: tem3_0(0:nx,0:ny,0:nz) ! automatic work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
INTEGER :: tstrtlvl
REAL :: deltat
INTEGER :: ebc1,wbc1,nbc1,sbc1
INTEGER :: nq0exb,nqdtexb
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Compute the advection term for a general water substance
! q and store the result in qadv.
!
!-----------------------------------------------------------------------
!
CALL set_acct(advs_acct)
IF(sadvopt == 4) THEN ! Forward-based FCT scheme
deltat = dtbig1
tstrtlvl = tpresent
ELSE
deltat = dtbig1*2
tstrtlvl = tpast
END IF
IF (sadvopt == 1 .OR. sadvopt == 2 .OR. sadvopt == 3 ) THEN
! 2nd or 4th order centerd sc
DO i=1,nx
DO j=1,ny
DO k=1,nz
tem7(i,j,k) = 0.0
END DO
END DO
END DO
CALL advq(nx,ny,nz,q,u,v,wcont,ustr,vstr,wstr, &
rhostr, tem7, mapfct, &
qadv, &
tem1,tem2,tem3,tem4)
ELSE IF( sadvopt == 4 .OR. sadvopt == 5 ) THEN ! FCT advection
CALL advqfct(nx,ny,nz,dtbig1,q,u,v,wcont, ustr,vstr,wstr, &
rhostr,rhostri,mapfct,j3,qadv, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem1_0,tem2_0,tem3_0)
END IF
!
!-----------------------------------------------------------------------
!
! Compute the mixing terms for the general water substance (q)
! equation, including both physical and computational mixing.
! Store the result in array qmix.
!
!-----------------------------------------------------------------------
!
CALL set_acct(cmix_acct)
CALL mixq(nx,ny,nz, &
q(1,1,1,tstrtlvl),rhostr,kmh,kmv,rprntl, &
x,y,z,zp,mapfct, j1,j2,j3,aj3x,aj3y,j3inv, &
qmix, &
tem1,tem2,tem3,tem4,tem5,tem6)
!
!-----------------------------------------------------------------------
!
! Call BRLXQ to added to qmix the additional boundary relaxation and
! spatial mixing in the boundary zone
!
! qmix = qmix + qexbc_term
!
!-----------------------------------------------------------------------
!
IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN
CALL acct_interrupt(bc_acct)
CALL brlxq(nx,ny,nz, deltat*0.5, qflag, q,rhostr, qmix, &
exbcbuf(nqv0exb), exbcbuf(nqc0exb), &
exbcbuf(nqr0exb), exbcbuf(nqi0exb), &
exbcbuf(nqs0exb), exbcbuf(nqh0exb), &
exbcbuf(nqvdtexb),exbcbuf(nqcdtexb), &
exbcbuf(nqrdtexb),exbcbuf(nqidtexb), &
exbcbuf(nqsdtexb),exbcbuf(nqhdtexb),bcrlx, &
tem1,tem2,tem3,tem4)
CALL acct_stop_inter
END IF
!
!-----------------------------------------------------------------------
!
! Calculate qforce, store the result in tem1.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=-qadv(i,j,k)+qmix(i,j,k)
IF(qflag == 2) THEN
tem1(i,j,k)= tem1(i,j,k)+rhostr(i,j,k)*qccumsrc(i,j,k)
ELSE IF(qflag == 3) THEN
tem1(i,j,k)= tem1(i,j,k)+rhostr(i,j,k)*qrcumsrc(i,j,k)
ELSE IF(qflag == 4) THEN
tem1(i,j,k)= tem1(i,j,k)+rhostr(i,j,k)*qicumsrc(i,j,k)
ELSE IF(qflag == 5) THEN
tem1(i,j,k)= tem1(i,j,k)+rhostr(i,j,k)*qscumsrc(i,j,k)
ELSE
END IF
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Treat the vertically implicit mixing term.
!
!-----------------------------------------------------------------------
!
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
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,q(1,1,1,tstrtlvl),rhostr, &
tem2,tem1,tem3,tem4,tem5,tem6)
END IF
!
!-----------------------------------------------------------------------
!
! Integrate forward by one timestep the general water substance
! (q) equation, yielding q at time = tfuture.
!
!-----------------------------------------------------------------------
!
CALL set_acct(tinteg_acct)
IF( nestgrd == 1 .AND. mgrid /= 1 ) THEN
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-2
q(i,j,k,tfuture)=q(i,j,k,tstrtlvl) &
+deltat*tem1(i,j,k)*rhostri(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set the boundary conditions on q for an adaptive (nested)
! grid run. If using only one grid, this portion of the code is
! skipped....proceed to next comment block.
!
!-----------------------------------------------------------------------
!
DO k=2,nz-2
DO i=1,nx-1
q(i, 1,k,tfuture)=2*q(i, 1,k,tpresent)-q(i, 1,k,tpast)
q(i,ny-1,k,tfuture)=2*q(i,ny-1,k,tpresent)-q(i,ny-1,k,tpast)
END DO
END DO
DO k=2,nz-2
DO j=2,ny-2
q( 1,j,k,tfuture)=2*q( 1,j,k,tpresent)-q( 1,j,k,tpast)
q(nx-1,j,k,tfuture)=2*q(nx-1,j,k,tpresent)-q(nx-1,j,k,tpast)
END DO
END DO
! bc's for mp not implemented for nested grids
CALL acct_interrupt(bc_acct)
CALL bcq(nx,ny,nz, dtbig1, q, qdteb,qdtwb,qdtnb,qdtsb, &
0,0,0,0,tbc,bbc)
CALL acct_stop_inter
ELSE
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
q(i,j,k,tfuture)=q(i,j,k,tstrtlvl) &
+deltat*tem1(i,j,k)*rhostri(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set the boundary conditions on q for a NON-adaptive (uniform)
! grid run.
!
!-----------------------------------------------------------------------
!
!call test_dump (q,"XXXsolve_q",nx,ny,nz,0,0)
IF ( lbcopt == 1 ) THEN
ebc1=ebc
wbc1=wbc
sbc1=sbc
nbc1=nbc
IF( ebc == 4 ) ebc1=0
IF( wbc == 4 ) wbc1=0
IF( nbc == 4 ) nbc1=0
IF( sbc == 4 ) sbc1=0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew( CALL mprecv2dew( CALL mpsend2dns( CALL mprecv2dns( END IF
CALL acct_interrupt(bc_acct)
CALL bcq(nx,ny,nz, dtbig1, q, qdteb,qdtwb,qdtnb,qdtsb, &
ebc1,wbc1,nbc1,sbc1,tbc,bbc)
CALL acct_stop_inter
ELSE
ebc1 = 3
wbc1 = 3
sbc1 = 3
nbc1 = 3
IF( ebc == 0 ) ebc1 = 0
IF( wbc == 0 ) wbc1 = 0
IF( nbc == 0 ) nbc1 = 0
IF( sbc == 0 ) sbc1 = 0
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew( CALL mprecv2dew( CALL mpsend2dns( CALL mprecv2dns( END IF
CALL acct_interrupt(bc_acct)
CALL bcq(nx,ny,nz, dtbig1, q, qdteb,qdtwb,qdtnb,qdtsb, &
ebc1,wbc1,nbc1,sbc1,tbc,bbc)
! Zero-gradient condition will be
! reset if exbc data is available for q
IF ( qflag == 1 ) THEN
nq0exb = nqv0exb
nqdtexb = nqvdtexb
ELSE IF ( qflag == 2 ) THEN
nq0exb = nqc0exb
nqdtexb = nqcdtexb
ELSE IF ( qflag == 3 ) THEN
nq0exb = nqr0exb
nqdtexb = nqrdtexb
ELSE IF ( qflag == 4 ) THEN
nq0exb = nqi0exb
nqdtexb = nqidtexb
ELSE IF ( qflag == 5 ) THEN
nq0exb = nqs0exb
nqdtexb = nqsdtexb
ELSE IF ( qflag == 6 ) THEN
nq0exb = nqh0exb
nqdtexb = nqhdtexb
END IF
CALL exbcq(nx,ny,nz, qflag, curtim+dtbig,q(1,1,1,tfuture), &
exbcbuf(nq0exb),exbcbuf(nqdtexb))
CALL acct_stop_inter
END IF
END IF
!call test_dump (q,"XXXAsolve_q",nx,ny,nz,0,1)
RETURN
END SUBROUTINE solvq
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE UPRAD1 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE uprad1(nx,ny,w1,w2,wc,nrho,ifax1,ifax2,trigs1,trigs2, & 1,8
temxy1,work)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! To apply the upper radiation condition using a periodic 1 or 2-D
! fft. The variable temxy1 contains input information and is
! overwritten with the final w (at nz-1) on exit.
!
! In the case of the open boundary, a linear hydrostatic relation
! between vertical velocity and pressure is used in fourier space.
! The use of fast fourier transform is required for this condition.
!
! The fft routine used is FFT991 obtained from NCAR. It is a
! version of the one used by ECMWF. Before fft991 is called, set99
! is called to setup the variables needed for the transform routine.
! (init3d.f)
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber
! 07/22/1997.
!
! Modification History:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction.
! ny Number of grid points in the y-direction.
!
! w1 Pre-fft w(nz-1) coefficient.
! w2 Pre-fft w(nz-2) coefficient.
! wc Pre-fft w(nz-1) coefficient from the reduction
! phase of the tridiaginal solver.
!
! trigs1 Array containing pre-computed trig function for
! fftopt=1.
! trigs2 Array containing pre-computed trig function for
! fftopt=1.
! ifax1 Array containing the factors of nx for fftopt=1.
! ifax2 Array containing the factors of ny for fftopt=1.
!
! OUTPUT:
!
!
! temxy1 On input known forcing term, on output final w at nz-1.
!
!
! WORK ARRAYS:
!
!
! work 2-D work array for fft code.
!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx ! Number of grid points in the x-direction.
INTEGER :: ny ! Number of grid points in the y-direction.
REAL :: w1 ! pre-fft w(nz-1) coefficient.
REAL :: w2 ! pre-fft w(nz-2) coefficient.
REAL :: wc ! pre-fft w(nz-1) coefficient from the
! reduction phase of the tridiagonal solver.
REAL :: trigs1(3*(nx-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
REAL :: trigs2(3*(ny-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
INTEGER :: ifax1(13) ! Array containing the factors of nx for
! fftopt=1.
INTEGER :: ifax2(13) ! Array containing the factors of ny for
! fftopt=1.
REAL :: temxy1 (nx+1,ny+1) ! On input known forcing term,
! on output final w at nz-1.
REAL :: work (ny,nx) ! 2-D work array for fft code.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,itema,itemb
REAL :: eps,kx,ky,pi,nrho,tema,temb,temc,temd,teme,temf,temg
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'bndry.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!-----------------------------------------------------------------------
!
! Perform the fft on temxy1
!
!-----------------------------------------------------------------------
IF(ny == 4)THEN ! 1-d transform in x, note nx must be ODD!
itema = 1
itemb = nx
CALL fft991(temxy1,work,trigs1,ifax1,1,nx+1,nx-1,1,-1)
ELSE IF(nx == 4)THEN ! 1-d transform in y, note ny must be ODD!
itema = ny
itemb = 1
CALL fft991(temxy1,work,trigs2,ifax2,nx+1,1,ny-1,1,-1)
ELSE ! Do 2-d transform, note nx,ny must be ODD!
itema = ny
itemb = nx
CALL fft991(temxy1,work,trigs1,ifax1,1,nx+1,nx-1,ny-1,-1)
CALL fft991(temxy1,work,trigs2,ifax2,nx+1,1,ny-1,nx,-1)
END IF ! end of run type if block.
!-----------------------------------------------------------------------
!
! Compute the wave space w at nz-1.
!
!-----------------------------------------------------------------------
eps= 1.0E-8
pi = 4.0*ATAN(1.0)
tema = pi/(nx-2)
temb = pi/(ny-2)
temc = 2.0*dxinv
temd = 2.0*dyinv
temf = w1-w2*wc
temg = eps - nrho
DO j=1,itema ! note: kx and ky are in finite difference form.
ky = temd*SIN(INT((j-1)*0.5+0.001)*temb)
ky = ky*ky
DO i=1,itemb
kx = temc*SIN(INT((i-1)*0.5+0.001)*tema)
teme =SQRT(kx*kx+ky)
temxy1(i,j)=-teme*temxy1(i,j)/(teme*temf+temg)
END DO
END DO
!-----------------------------------------------------------------------
!
! Perform the reverse fft on temxy1 to obtain w(nz-1) in real space.
!
!-----------------------------------------------------------------------
IF(ny == 4)THEN ! 1-d transform in x, note nx must be ODD!
CALL fft991(temxy1,work,trigs1,ifax1,1,nx+1,nx-1,1,1)
DO j=2,ny-1
DO i=1,nx-1
temxy1(i,j) = temxy1(i,1)
END DO
END DO
ELSE IF(nx == 4)THEN ! 1-d transform in y, note ny must be ODD!
CALL fft991(temxy1,work,trigs2,ifax2,nx+1,1,ny-1,1,1)
DO j=1,ny-1
DO i=2,nx-1
temxy1(i,j) = temxy1(1,j)
END DO
END DO
ELSE ! Do 2-d transform, note nx,ny must be ODD!
CALL fft991(temxy1,work,trigs2,ifax2,nx+1,1,ny-1,nx, 1)
CALL fft991(temxy1,work,trigs1,ifax1,1,nx+1,nx-1,ny-1, 1)
END IF ! end of run type if block.
RETURN
END SUBROUTINE uprad1
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE UPRAD3 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE uprad3(nx,ny,w1,w2,wc,nrho,wsave1,wsave2,vwork1,vwork2, & 1,8
temxy1,work)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! To apply the upper radiation condition using an even (vfftpack)
! sequence fft. The variable temxy1 contains input information and
! is overwritten with the updated w (at nz-1) on exit.
!
! In the case of the open boundary, a linear, hydrostatic relation
! between vertical velocity to pressure is used in fourier space.
! The use of fast fourier transform is NOT required for this
! option but is recommended. If nx and ny are not of special
! character, then a basic fourier transform is used.
!
! The fft routine used is vcost from the vfftpack software suite
! available at PSC.
!
! Before vcost is called, the arrays wsave1,wsave2,vwork1,vwork2
! must be initialized (init3d.f).
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber
! 07/22/1997.
!
!
! Modification History:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction.
! ny Number of grid points in the y-direction.
!
! w1 Pre-fft w(nz-1) coefficient.
! w2 Pre-fft w(nz-2) coefficient.
! wc Pre-fft w(nz-1) coefficient from the reduction
! phase of the tridiaginal solver.
!
! vwork1 2-D work array for fftopt=2 option.
! vwork2 2-D work array for fftopt=2 option.
! wsave1 Work array for fftopt=2 option.
! wsave2 Work array for fftopt=2 option.
!
! OUTPUT:
!
!
! temxy1 On input known forcing term, on output final w at nz-1.
!
!
! WORK ARRAYS:
!
!
! work 2-D work array for fft code.
!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx ! Number of grid points in the x
! direction.
INTEGER :: ny ! Number of grid points in the y
! direction.
REAL :: w1 ! Pre-fft w(nz-1) coefficient.
REAL :: w2 ! Pre-fft w(nz-2) coefficient.
REAL :: wc ! Pre-fft w(nz-1) coefficient from the
! reduction phase of the tridiagonal
! solver.
REAL :: vwork1 (nx+1,ny+1) ! 2-D work array for fftopt=2.
REAL :: vwork2 (ny,nx+1) ! 2-D work array for fftopt=2.
REAL :: wsave1 (3*(ny-1)+15) ! Work array for fftopt=2.
REAL :: wsave2 (3*(nx-1)+15) ! Work array for fftopt=2.
REAL :: temxy1 (nx+1,ny+1) ! On input known forcing term,
! on output final w at nz-1.
REAL :: work (ny,nx) ! 2-D work array (global)
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'bndry.inc'
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j
REAL :: kx,ky,pi,tema,temb,temc,temd,teme,temf
REAL :: eps,nrho
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!-----------------------------------------------------------------------
!
! Perform the fft on temxy1.
!
!-----------------------------------------------------------------------
IF(ny == 4)THEN ! 1-d transform in x.
DO i=1,nx-1 ! swap temxy1(i,j) array into work(j,i)
work(1,i) = temxy1(i,1)
END DO
CALL vcost(1,nx-1,work,vwork2,ny,wsave2)
DO i=1,nx-1 ! swap back for computations....
temxy1(i,1) = work(1,i)
END DO
ELSE IF(nx == 4)THEN ! 1-d transform in y.
CALL vcost(1,ny-1,temxy1,vwork1,nx+1,wsave1)
ELSE ! Do 2-d transform.
CALL vcost(nx-1,ny-1,temxy1,vwork1,nx+1,wsave1)
DO j=1,ny-1 ! swapping arrays.
DO i=1,nx-1
work(j,i) = temxy1(i,j)
END DO
END DO
CALL vcost(ny-1,nx-1,work,vwork2,ny,wsave2)
DO j=1,ny-1 ! swapping arrays for computations.
DO i=1,nx-1
temxy1(i,j) = work(j,i)
END DO
END DO
END IF ! end of run type if block.
!-----------------------------------------------------------------------
!
! Compute the wave space w at nz-1.
!
!-----------------------------------------------------------------------
eps = 1.0E-13
pi = 4.0*ATAN(1.0)
tema = 0.5*pi/(nx-2) ! finite difference form.
temb = 0.5*pi/(ny-2) ! finite difference form.
temc = 2.0*dxinv
temd = 2.0*dyinv
teme = w1-w2*wc
temf = eps - nrho
temxy1(1,1) = -temxy1(1,1)/(teme+temf) ! first part of 2-d case.
IF(ny == 4)THEN ! perform computations in x direction.
DO i=2,nx-1 ! compute from i=2,nx-1 and j=1...
kx = temc*SIN(INT(i-1)*tema)
kx = SQRT(kx*kx)
temxy1(i,1)=-kx*temxy1(i,1)/(kx*teme+temf)
END DO
ELSE IF(nx == 4)THEN ! perform computations in y direction.
DO j=2,ny-1 ! compute from j=2,ny-1 and i=1.
ky = temd*SIN(INT(j-1)*temb)
ky = SQRT(ky*ky)
temxy1(1,j)=-ky*temxy1(1,j)/(ky*teme+temf)
END DO
ELSE ! compute the 2-D version.
DO i=2,nx-1 ! compute from i=2,nx-1 and j=1.
kx = temc*SIN(INT(i-1)*tema)
kx = SQRT(kx*kx)
temxy1(i,1)=-kx*temxy1(i,1)/(kx*teme+temf)
END DO
DO j=2,ny-1 ! compute from i=1,j=2,ny-1.
ky = temd*SIN(INT(j-1)*temb)
ky = SQRT(ky*ky)
temxy1(1,j)=-ky*temxy1(1,j)/(ky*teme+temf)
END DO
DO j=2,ny-1 ! compute from i=2,nx-1,j=2,ny-1.
ky = temd*SIN(INT(j-1)*temb)
ky = ky*ky
DO i=2,nx-1
kx = temc*SIN(INT(i-1)*tema)
kx =SQRT(kx*kx+ky)
temxy1(i,j)=-kx*temxy1(i,j)/(kx*teme+temf)
END DO
END DO
END IF ! end of the run type if block.
!-----------------------------------------------------------------------
!
! Perform the reverse fft on temxy1 to obtain w in real space.
!
!-----------------------------------------------------------------------
IF(ny == 4)THEN ! 1-d transform in x, note nx must be EVEN!
DO i=1,nx-1
work(1,i) = temxy1(i,1)
END DO
CALL vcost(1,nx-1,work,vwork2,ny,wsave2)
DO j=1,ny-1
DO i=1,nx-1
temxy1(i,j) = work(1,i)
END DO
END DO
ELSE IF(nx == 4)THEN ! 1-d transform in y, note ny must be EVEN!
CALL vcost(1,ny-1,temxy1,vwork1,nx+1,wsave1)
DO j=1,ny-1
DO i=1,nx-1
temxy1(i,j) = temxy1(1,j)
END DO
END DO
ELSE ! Do 2-d transform, note nx,ny must be EVEN!
DO j=1,ny-1 ! swapping arrays.
DO i=1,nx-1
work(j,i) = temxy1(i,j)
END DO
END DO
CALL vcost(ny-1,nx-1,work,vwork2,ny,wsave2)
DO j=1,ny-1 ! swapping arrays.
DO i=1,nx-1
temxy1(i,j) = work(j,i)
END DO
END DO
CALL vcost(nx-1,ny-1,temxy1,vwork1,nx+1,wsave1)
END IF ! end of runopt if block, w (nz-1) is updated.
RETURN
END SUBROUTINE uprad3
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE PPRTBBC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE pprtbbc(nx,ny,nz,g05,buoy2swt,rhostr,pprt,ptprt, & 6,8
pbari,ptbari,phydro, &
pprt1,tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! To compute the perturbation pressure below the ground (scalar k=1)
! The technique uses the perturbation hydrostatic relation obtained
! from the model vertical momentum equation. The technique is
! summarized here.
!
! The calculation of the hydrostatic pressure bottom boundary
! condition involves two steps. The first step is used to
! approximate the second order pprt term in the buoyancy term
! and the second step obtains the final pprt at k=1.
!
! Step A: pprt(1)=D=(pprt(2) - F - dz*phydro + 2.0*A*B)/(1-C)
!
! where, A = dz*g*avgz(rhobar*J3)/(2*cpdcv)
! B = pprt(i,j,2)/pbar(i,j,2)
! C = A/pbar(i,j,1)
! F = alpha*div*(1)-alpha*div*(2) not included....
! note:
! D = the intermediate result pprt(1) for use in the
! second order pprt term
!
! Step B: pprt(1) = D + A*(1/cpdcv -1)*(B*B + E*E)/(1.0-C)
!
! where, E = D/pbar(i,j,1)
!
! The result is stored in phydro(i,j) and passed into bcp to set
! pprt(i,j,1).
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber
! 11/06/1997.
!
!
! Modification History:
!
! 3/18/99 (D. Weber)
! Bug fix to second order terms.
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! 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
!
! ptprt Perturbation potential temperature at time tpresent (K)
! pprt Perturbation pressure at time tpresent (Pascal)
!
! phydro Big time step forcing term for use in computing the
! hydrostatic pressure at k=1.
!
! rhostr Base state density rhobar times j3 (kg/m**3)
! ptbari Inverse base state potential temperature (K)
! pbari Inverse base state pressure (Pascal)
!
! OUTPUT:
!
! pprt1 Holds pprt(i,j,1) via hydrostatic balance.
!
! WORK ARRAYS:
!
! tem2 Work array.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure at tfuture
! (Pascal)
REAL :: phydro(nx,ny) ! Big time step forcing for computing
! hydrostatic pprt at k=1.
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: ptbari(nx,ny,nz) ! Inverse base state pot. temperature (K)
REAL :: pbari (nx,ny,nz) ! Inverse base state pressure (Pascal).
REAL :: pprt1 (nx,ny) ! pprt1
!
!-----------------------------------------------------------------------
!
! Temporary WORK ARRAYS:
!
!-----------------------------------------------------------------------
!
REAL :: tem1 (nx,ny,nz) ! Temporary work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j
REAL :: g05
REAL :: tema,temb,a,b,c,d,e
INTEGER :: buoy2swt !Switch for 1st-order or 2nd-order in buoyancy
REAL :: ptemk,ptemk1,pttemk,pttemk1
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'bndry.inc'
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'phycst.inc' ! Physical constants
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF( ptsmlstp == 1 ) THEN ! add in the 1st and 2nd order pt terms.
tema = 0.5/cpdcv
DO j=1,ny-1
DO i=1,nx-1
ptemk = pprt(i,j,2)*pbari(i,j,2) ! new code.....
ptemk1= pprt(i,j,1)*pbari(i,j,1)
pttemk = ptprt(i,j,2)*ptbari(i,j,2)
pttemk1= ptprt(i,j,1)*ptbari(i,j,1)
temb = 0.5*(rhostr(i,j,1)+rhostr(i,j,2))
pprt1(i,j)=phydro(i,j)+temb*g05*(pttemk+pttemk1 &
-buoy2swt*(tema*(ptemk*pttemk+ptemk1*pttemk1) &
+pttemk*pttemk+pttemk1*pttemk1))
END DO
END DO
ELSE
DO j=1,ny-1
DO i=1,nx-1
pprt1(i,j)= phydro(i,j)
END DO
END DO
END IF
tema =buoy2swt*0.5*(1.0/cpdcv -1.0)
temb = dz*g*0.25/cpdcv
DO j=1,ny-1
DO i=1,nx-1
! compute the intermediate result...
a = temb*(rhostr(i,j,2)+rhostr(i,j,1))
b = pprt(i,j,2)*pbari(i,j,2)
c = a*pbari(i,j,1)
d = (pprt(i,j,2) - dz*pprt1(i,j) + a*b) / (1.0-c)
! compute the final pprt(i,j,1)....
e = d*pbari(i,j,1)
pprt1(i,j) = d + a*tema*(b*b + e*e) / (1.0-c)
END DO
END DO
!call test_dump (pprt1,"XXXsolve_pprt1",nx,ny,1,0,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend1dew(pprt1,nx,ny,ebc,wbc,0,mptag,tem1)
CALL mprecv1dew(pprt1,nx,ny,ebc,wbc,0,mptag,tem1)
CALL mpsend1dns(pprt1,nx,ny,nbc,sbc,0,mptag,tem1)
CALL mprecv1dns(pprt1,nx,ny,nbc,sbc,0,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
CALL bcs2d(nx,ny,pprt1, ebc,wbc,nbc,sbc)
CALL acct_stop_inter
!call test_dump (pprt1,"XXXAsolve_pprt1",nx,ny,1,0,1)
RETURN
END SUBROUTINE pprtbbc