!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ADVUVW ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE advuvw(nx,ny,nz,u,v,w,wcont,rhostr,ubar,vbar, mapfct, & 1,4
uforce,vforce,wforce,uadv,vadv,wadv, &
ustr,vstr,wstr,tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinates the calculation of the advection terms uadv, vadv and
! wadv of the u, v and w equations. These terms are written in
! equivalent advection form.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 11/20/91.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (M. Xue)
! Added full documentation.
!
! 5/29/92 (K. Brewster)
! Further facelift.
!
! 4/9/93 (M. Xue & K. Brewster)
! Some index bounds for operators and loop bounds were corrected.
! Redundant calculations were done before at the boundaries that
! caused out-of-bound calculations, which should not have
! affected the results in normal situations though.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 2/9/94 (D. Jahn)
! Add tem3 work array in conjunction w/ changes to advu,advv,
! advw, and advcts.
!
! 9/1/94 (D. Weber & Y. Lu)
! Cleaned up documentation
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 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/11/98 (D. Weber)
! Removed most single operation loops and merged the advection
! contributions into the forcing arrays prior to the return.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! mapfct Map factors at scalar, u and v points
!
! u x component of velocity at all time levels (m/s)
! v y component of velocity at all time levels (m/s)
! w Vertical component of velocity in Cartesian
! coordinates at all time levels (m/s).
! wcont Contravariant vertical velocity (m/s)
! rhostr Base state density rhobar times j3 (kg/m**3)
! ubar Base state x component of velocity (m/s)
! vbar Base state y component of velocity (m/s)
!
! OUTPUT:
!
! uadv u-equation advection terms (kg/m**3)*(m/s)/s
! vadv v-equation advection terms (kg/m**3)*(m/s)/s
! wadv w-equation advection terms (kg/m**3)*(m/s)/s
!
! WORK ARRAYS:
!
! ustr u * rhostr
! vstr v * rhostr
! wstr w * rhostr
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'timelvls.inc'
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: mapfct(nx,ny,8) ! Map factors at scalar, u and v points
REAL :: u (nx,ny,nz,nt) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity (m/s)
REAL :: w (nx,ny,nz,nt) ! Total w-velocity (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s)
REAL :: uadv (nx,ny,nz) ! u-eqn advection terms
! (kg/m**3)*(m/s)/s
REAL :: vadv (nx,ny,nz) ! v-eqn advection terms
! (kg/m**3)*(m/s)/s
REAL :: wadv (nx,ny,nz) ! w-eqn advection terms
! (kg/m**3)*(m/s)/s
REAL :: uforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in u-momentum equation (kg/(m*s)**2)
! uforce= -uadv + umix + ucorio
REAL :: vforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in v-momentum equation (kg/(m*s)**2)
! vforce= -vadv + vmix + vcorio
REAL :: wforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in w-momentum equation (kg/(m*s)**2)
! wforce= -wadv + wmix + wbuoy
REAL :: ustr (nx,ny,nz) ! u * rhostr
REAL :: vstr (nx,ny,nz) ! v * rhostr
REAL :: wstr (nx,ny,nz) ! w * rhostr
REAL :: tem1 (nx,ny,nz) ! Temporary work array
REAL :: tem2 (nx,ny,nz) ! Temporary work array
REAL :: tem3 (nx,ny,nz) ! Temporary work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
REAL :: tema
INTEGER :: tlevel
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
tlevel=tpresent
!
!-----------------------------------------------------------------------
!
! Calculate ustr=rhostr*u, vstr=rhostr*v, wstr=rhostr*w
!
!-----------------------------------------------------------------------
!
CALL uvwrho
(nx,ny,nz,u(1,1,1,tlevel),v(1,1,1,tlevel), &
wcont,rhostr, &
ustr,vstr,wstr)
!
!-----------------------------------------------------------------------
!
! Calculate uadv:
!
!-----------------------------------------------------------------------
!
CALL advu(nx,ny,nz, u,v,w, &
ustr,vstr,wstr, ubar, mapfct(1,1,2), uadv, uforce, &
tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! Calculate vadv:
!
!-----------------------------------------------------------------------
!
CALL advv(nx,ny,nz, u,v,w, &
ustr,vstr,wstr, vbar, mapfct(1,1,3), vadv, vforce, &
tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! Calculate wadv:
!
!-----------------------------------------------------------------------
!
CALL advw(nx,ny,nz, u,v,w, &
ustr,vstr,wstr, mapfct(1,1,1), wadv, wforce, &
tem1,tem2,tem3)
RETURN
END SUBROUTINE advuvw
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ADVP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE advp(nx,ny,nz,pprt,u,v,w,wcont,rhostr, mapfct, & 1,1
j3,aj3x,aj3y,aj3z, &
padv, &
uj3,vj3,wj3,tem1,tem2,tem3,tem4,mp_tem)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the advection of perturbation pressure.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 11/20/91.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (M. Xue)
! Added full documentation.
!
! 5/29/92 (K. Brewster)
! Further facelift.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 5/25/93 (M. Xue)
! Fixed the bug in the base state vertical pressure term.
! A j3 was missing before.
!
! 5/26/93 (M. Xue)
! Corrected pbar advection term.
!
! 9/1/94 (D. Weber & Y. Lu)
! Cleaned up documentation
!
! 9/7/94 (M. Xue)
! This subroutine now calculates the perturbation pressure advection
! only, by calling subroutine ADVCTS.
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 9/18/98 (D. Weber)
! Added aj3x,y,z arrays.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! pprt Perturbation pressure at all time levels (Pascal)
! u x component of velocity at all time levels (m/s)
! v y component of velocity at all time levels (m/s)
! w Vertical component of velocity in Cartesian
! coordinates at al time levels (m/s).
! wcont Contravariant vertical velocity (m/s)
!
! rhostr Base state air density (kg/m**3)
!
! mapfct Map factors at scalar points
!
! aj3x Avgx of the coordinate transformation Jacobian d(zp)/dz
! aj3y Avgy of the coordinate transformation Jacobian d(zp)/dz
! aj3z Avgz of the coordinate transformation Jacobian d(zp)/dz
!
! OUTPUT:
!
! padv Pressure equation advection terms (kg/(m*s**3))
!
! WORK ARRAYS:
!
! uj3 Temporary work array.
! vj3 Temporary work array.
! wj3 Temporary work array.
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
! tem4 Temporary work array.
! mp_tem Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'timelvls.inc'
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: u (nx,ny,nz,nt) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity (m/s)
REAL :: w (nx,ny,nz,nt) ! Total w-velocity (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: pprt (nx,ny,nz,nt) ! Perturbation pressure from that
! of base state atmosphere (Pascal)
REAL :: rhostr(nx,ny,nz) ! Base state air density (kg/m**3)
REAL :: mapfct(nx,ny) ! Map factors at scalar points
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 :: padv (nx,ny,nz) ! Advection term of the pressure
! equation (kg/(m*s**3)
REAL :: uj3 (nx,ny,nz) ! Temporary work array to carry u*j3
REAL :: vj3 (nx,ny,nz) ! Temporary work array to carry v*j3
REAL :: wj3 (nx,ny,nz) ! Temporary work array to carry wcont*j3
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 :: mp_tem(MAX(nx,ny)*nz) ! Temporary message passing array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
INTEGER :: tlevel
REAL :: tema
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
tlevel=tpresent
!
!-----------------------------------------------------------------------
!
! Advection of the perturbation pressure
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Calculate the advection of perturbation pressure.
! (j3*u*dp/dx+j3*v*dp/dy+j3*wcont*dp/dz) =
! avgx( avgx(j3) * u * difx(p) )
! + avgy( avgy(j3) * u * dify(p) )
! + avgz( avgz(j3) * wcont * difz(p) )
!
! Subroutine ADVCTS is called to evaluate this term.
!
!-----------------------------------------------------------------------
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx
uj3(i,j,k)=u(i,j,k,tlevel)*aj3x(i,j,k)
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny
DO i=1,nx-1
vj3(i,j,k)=v(i,j,k,tlevel)*aj3y(i,j,k)
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
wj3(i,j,k)=wcont(i,j,k)*aj3z(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Perturbation potential temperature advection
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem4(i,j,k) = 0.0
END DO
END DO
END DO
CALL advcts(nx,ny,nz, pprt, &
u(1,1,1,tlevel),v(1,1,1,tlevel), &
uj3,vj3,wj3,tem4, mapfct, &
padv, tem1,tem2,tem3,mp_tem)
!
!-----------------------------------------------------------------------
!
! Note that the base state pressure advection term is evaluated
! inside the small time step, therefore padv does not include
! this part on exit of this subroutine.
!
!-----------------------------------------------------------------------
!
RETURN
END SUBROUTINE advp
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ADVQ ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE advq(nx,ny,nz,q,u,v,wcont,ustr,vstr,wstr, & 3,1
rhostr,qbar, mapfct, &
qadv, &
tem1,tem2,tem3,mp_tem)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate water/ice advection.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 05/20/92.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (M. Xue)
! Added full documentation.
!
! 5/29/92 (K. Brewster)
! Further facelift and order of argument list changed.
!
! 2/9/94 (D. Jahn)
! Add tem3 to the call of advcts.
!
! 9/1/94 (D. Weber & Y. Lu)
! Cleaned up documentation
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 9/28/98 (Dan Weber)
! Moved u,v,wstr computations outside of this routine.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! mapfct Map factors at scalar points
!
! q Water or ice mixing ratio (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)
! rhostr Base state air density (kg/m**3)
! qbar Base state of water oy ice mixing ratio (kg/kg)
!
! OUTPUT:
!
! qadv Water or ice mixing ratio advection term in the
! conservation equation (kg/m**3)*(kg/kg)/s
!
! WORK ARRAYS:
!
! ustr u * rhostr
! vstr v * rhostr
! wstr w * rhostr
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
! mp_tem Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'timelvls.inc'
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: mapfct(nx,ny) ! Map factors at scalar points
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 :: rhostr(nx,ny,nz) ! Base state air density (kg/m**3)
REAL :: qbar (nx,ny,nz) ! Base state of the water/ice variable
REAL :: qadv (nx,ny,nz) ! Advection of a water/ice variable
REAL :: ustr (nx,ny,nz) ! u * rhostr
REAL :: vstr (nx,ny,nz) ! v * rhostr
REAL :: wstr (nx,ny,nz) ! w * rhostr
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 :: mp_tem(MAX(nx,ny)*nz) ! Temporary message passing array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: tlevel
INTEGER :: i, j, k
REAL :: tema
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global model control parameters
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
tlevel=tpresent
!
!-----------------------------------------------------------------------
!
! Calculate ustr=rhostr*u, vstr=rhostr*v, wstr=rhostr*w
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! Advection of the water/ice mixing ratio
!
!-----------------------------------------------------------------------
CALL advcts(nx,ny,nz, q, &
u(1,1,1,tlevel),v(1,1,1,tlevel), &
ustr,vstr,wstr, qbar, mapfct, &
qadv, tem1,tem2,tem3,mp_tem)
RETURN
END SUBROUTINE advq
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ADVU ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE advu(nx,ny,nz,u,v,w,ustr,vstr,wstr, ubar, mapfct, & 1,10
uadv, uforce, &
tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the advection terms of the u-equation. These terms are
! written in equivalent advection form.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 11/20/91.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (M. Xue)
! Added full documentation.
!
! 5/29/92 (K. Brewster)
! Further facelift.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 2/9/94 (D. Jahn)
! Add 4th-order momentum advection.
!
! 9/1/94 (D. Weber & Y. Lu)
! Cleaned up documentation
!
! 11/5/1995 (M. Xue)
! Added vertical fourth order advection for u.
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 9/1/98 (D. Weber and T. Chung)
! Removed operators and merged loops 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
!
! mapfct Map factors at u points
!
! u x component of velocity at a given time level (m/s)
! v y component of velocity at all time levels (m/s)
! w Vertical component of velocity in Cartesian
! coordinates at all time levels (m/s).
!
! ustr u * rhostr
! vstr v * rhostr
! wstr w * rhostr
!
! ubar Base state x component of velocity (m/s)
!
! OUTPUT:
!
! uforce Acoustically inactive forcing terms in u-momentum
! equation (kg/(m*s)**2).
!
! WORK ARRAYS:
!
! uadv u-equation advection terms (kg/m**3)*(m/s)/s
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'timelvls.inc'
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: mapfct(nx,ny) ! Map factors at u points
REAL :: u (nx,ny,nz,nt) ! Total u-velocity at a given time level
! (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity at a given time level
! (m/s)
REAL :: w (nx,ny,nz,nt) ! Total w-velocity at a given time level
! (m/s)
REAL :: ustr (nx,ny,nz) ! u*rhostr
REAL :: vstr (nx,ny,nz) ! v*rhostr
REAL :: wstr (nx,ny,nz) ! w*rhostr
REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s)
REAL :: uadv (nx,ny,nz) ! Advection term of the u-eq.
! (kg/m**3)*(m/s)/s
REAL :: uforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in u-momentum equation (kg/(m*s)**2)
REAL :: tem1 (nx,ny,nz) ! Temporary work array
REAL :: tem2 (nx,ny,nz) ! Temporary work array
REAL :: tem3 (nx,ny,nz) ! Temporary work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,onvf
INTEGER :: ibgn,iend,jbgn,jend
REAL :: vsb, vnb
REAL :: fourthirds,foursixth,onesix
REAL :: tema
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global model control parameters
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'bndry.inc' ! Parameters for boundary conditions.
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
fourthirds = 4.0/3.0
foursixth = 4.0/6.0
onesix = 1.0/6.0
!-----------------------------------------------------------------------
!
! Calculate rho*u* du/dx = avgx( avgx( ustr )*difx( u ) ):
!
!-----------------------------------------------------------------------
tema = 0.25*dxinv
DO k=2,nz-2
DO j=1,ny-1
DO i=2,nx-1
uadv(i,j,k)=tema*((ustr(i+1,j,k)+ustr(i,j,k))* &
(u(i+1,j,k,2)-u(i,j,k,2)) &
+(ustr(i,j,k)+ustr(i-1,j,k))* &
(u(i,j,k,2)-u(i-1,j,k,2)))
END DO
END DO
END DO
!call test_dump (uadv,"XXX2advu_uadv",nx,ny,nz,1,0)
!call test_dump (ustr,"XXXadvu_ustr",nx,ny,nz,1,0)
!-----------------------------------------------------------------------
!
! tem2 contains avgx((avgx(u*rho))*du/dx) (second order term)
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! If momentum advection 4th-order, calculate adjustment term:
! avg2x(avg2x(ustr)*dif2x(u))
!
!-----------------------------------------------------------------------
IF (madvopt == 2 .OR. madvopt == 3) THEN
tema = 0.25*dxinv
DO k=2,nz-2
DO j=1,ny-1
DO i=2,nx-1
tem2(i,j,k)=tema*(ustr(i+1,j,k)+ustr(i-1,j,k))* &
(u(i+1,j,k,2)-u(i-1,j,k,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! At this point tem2 contains avg2x(ustr)*dif2x(u)
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! In the case of periodic boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
! call test_dump (tem2,"XXXadvct_tem2",nx,ny,nz,1,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(tem2,nx,ny,nz,ebc,wbc,1,mptag,tem1)
CALL mprecv2dew(tem2,nx,ny,nz,ebc,wbc,1,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
ibgn = 3
iend = nx-2
jbgn = 1
jend = ny-1
IF (wbc == 0) ibgn=2
IF (wbc == 2) THEN
ibgn=2
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO j=jbgn,jend
tem2(1,j,k) = tem2(nx-2,j,k)
END DO
END DO
END IF
END IF
IF (ebc == 0) iend = nx-1
IF (ebc == 2) THEN
iend = nx-1
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO j=jbgn,jend
tem2(nx,j,k) = tem2(3,j,k)
END DO
END DO
END IF
END IF
!-----------------------------------------------------------------------
!
! In the case of rigid boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
IF (wbc == 1) THEN
ibgn=2
DO k=2,nz-2
DO j=jbgn,jend
tem2(1,j,k) = -tem2(3,j,k)
END DO
END DO
END IF
IF (ebc == 1) THEN
iend = nx-1
DO k=2,nz-2
DO j=jbgn,jend
tem2(nx,j,k) = -tem2(nx-2,j,k)
END DO
END DO
END IF
! call test_dump (tem2,"XXXAadvct_tem2",nx,ny,nz,1,0) ! k=1,nz-1 not set
CALL acct_stop_inter
!-----------------------------------------------------------------------
!
! Concatenate 2nd-order and adjusting term for 4th-order
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=jbgn,jend
DO i=ibgn,iend
uadv(i,j,k) = fourthirds*uadv(i,j,k) - &
onesix*(tem2(i+1,j,k)+tem2(i-1,j,k))
END DO
END DO
END DO
!call test_dump (uadv,"XXX1advu_uadv",nx,ny,nz,1,0)
END IF
!-----------------------------------------------------------------------
!
! Calculate rho*v* du/dy = avgy( avgx( vstr )* dify( u ) )
!
!-----------------------------------------------------------------------
tema = 0.25*dyinv
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-1
tem2(i,j,k)=tema*((vstr(i,j,k)+vstr(i-1,j,k))* &
(u(i,j,k,2)-u(i,j-1,k,2)) &
+(vstr(i,j+1,k)+vstr(i-1,j+1,k))* &
(u(i,j+1,k,2)-u(i,j,k,2)))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! tem1 contains avgy( avgx( vstr )* dify( u ) ) second order term
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! For 4th-order momentum advection, calculate the adjustment term.
! avg2y( avgx( avgy (vstr) )* dif2y( u ))
!
!-----------------------------------------------------------------------
IF (madvopt == 2 .OR. madvopt == 3) THEN
DO k=2,nz-2 ! 2 bugs found here..multiply sb plus...
DO j=2,ny-2
DO i=2,nx-1
tem1(i,j,k)=0.25*((vstr(i,j+1,k)+vstr(i,j,k)) &
+(vstr(i-1,j+1,k)+vstr(i-1,j,k)))
END DO
END DO
END DO
tema = 0.5*dyinv
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-1
tem3(i,j,k)=tema*(u(i,j+1,k,2)-u(i,j-1,k,2))*tem1(i,j,k)
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! At this point tem3 contains avgx( avgy (vstr) )* dif2y( u ))
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! In the case of periodic boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
! call test_dump (tem3,"XXXadvct_tem3",nx,ny,nz,0,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dns(tem3,nx,ny,nz,nbc,sbc,0,mptag,tem1)
CALL mprecv2dns(tem3,nx,ny,nz,nbc,sbc,0,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
ibgn = 2
iend = nx-1
jbgn = 3
jend = ny-3
IF (nbc == 0) jend = ny-2
IF (nbc == 2) THEN
jend = ny-2
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,ny-1,k) = tem3(i,2,k)
END DO
END DO
END IF
END IF
IF (sbc == 0) jbgn = 2
IF (sbc == 2) THEN
jbgn = 2
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,1,k) = tem3(i,ny-2,k)
END DO
END DO
END IF
END IF
!-----------------------------------------------------------------------
!
! In the case of rigid boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
IF (nbc == 1) THEN
jend = ny-2
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,ny-1,k) = tem3(i,ny-2,k)
END DO
END DO
END IF
IF (sbc == 1) THEN
jbgn = 2
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,1,k) = tem3(i,2,k)
END DO
END DO
END IF
! call test_dump (tem3,"XXXAadvct_tem3",nx,ny,nz,0,0) ! k=1,nz-1 not set
CALL acct_stop_inter
!-----------------------------------------------------------------------
!
! concatenate 2nd-order and adjusting term for 4th-order
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=jbgn,jend
DO i=ibgn,iend
tem2(i,j,k) = fourthirds*tem2(i,j,k) - &
onesix*(tem3(i,j+1,k)+tem3(i,j-1,k))
END DO
END DO
END DO
END IF
!-----------------------------------------------------------------------
!
! Calculate rho*v* du/dy on the north and south boundaries using
! one sided advection.
!
!-----------------------------------------------------------------------
!call test_dump (vstr,"XXXadvu_vstr",nx,ny,nz,1,0)
!call test_dump (u,"XXXadvu_u",nx,ny,3*nz,1,0)
!call test_dump (ubar,"XXXadvu_ubar",nx,ny,nz,1,0)
DO k=2,nz-2
DO i=2,nx-1
vsb=(vstr(i-1,1,k)+vstr(i-1,2,k) &
+vstr(i,1,k)+vstr(i,2,k))*0.25
vnb=(vstr(i-1,ny-1,k)+vstr(i-1,ny,k) &
+vstr(i,ny-1,k)+vstr(i,ny,k))*0.25
IF (vsb < 0.0) THEN
tem2(i,1,k)=vsb*(u(i,2,k,tpresent)-u(i,1,k,tpresent))*dyinv
ELSE
tem2(i,1,k)= rlxlbc* vsb*(u(i,1,k,tpast)-ubar(i,1,k))*dyinv
END IF
IF (vnb > 0.0) THEN
tem2(i,ny-1,k)=vnb* &
(u(i,ny-1,k,tpresent)-u(i,ny-2,k,tpresent))*dyinv
ELSE
tem2(i,ny-1,k)=-rlxlbc*vnb* &
(u(i,ny-1,k,tpast)-ubar(i,ny-1,k))*dyinv
END IF
END DO
END DO
!-----------------------------------------------------------------------
!
! Add the u-equation horizontal advection terms, rho*u*du/dx (uadv) and
! rho*v*du/dy (tem1).
!
!-----------------------------------------------------------------------
!call test_dump (uforce,"XXXBadvu_uforce",nx,ny,nz,1,0)
DO k=2,nz-2
DO j=1,ny-1
DO i=2,nx-1
uforce(i,j,k)=uforce(i,j,k) - &
(uadv(i,j,k)+tem2(i,j,k)) * mapfct(i,j)
END DO
END DO
END DO
!call test_dump (uforce,"XXXadvu_uforce",nx,ny,nz,1,0)
!call test_dump (uadv,"XXXadvu_uadv",nx,ny,nz,1,0)
!call test_dump (tem2,"XXXadvu_tem2",nx,ny,nz,1,0)
!call test_dump (mapfct,"XXXadvu_mapfct",nx,ny,1,1,0)
!-----------------------------------------------------------------------
!
! Calculate rho*w* du/dz = avgz( avgx( wstr ) * difz ( u ) )
!
!-----------------------------------------------------------------------
tema=0.5*dzinv
DO k=2,nz-1
DO j=1,ny-1
DO i=2,nx-1
tem2(i,j,k)=tema*(wstr(i,j,k)+wstr(i-1,j,k))* &
(u(i,j,k,2)-u(i,j,k-1,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! tem2 contains avgz( avgx( wstr ) * difz ( u ) )
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! For 4th-order momentum advection, calculate the fourth order
! contribution avg2z( avgz( avgx (wstr) )* dif2z( u ))
!
!-----------------------------------------------------------------------
IF (madvopt == 3) THEN
tema = 0.125*dzinv
DO k=2,nz-2
DO j=1,ny-1
DO i=2,nx-1
tem3(i,j,k)=tema*((wstr(i,j,k+1)+wstr(i-1,j,k+1))+ &
(wstr(i,j,k)+wstr(i-1,j,k))) * &
(u(i,j,k+1,2)-u(i,j,k-1,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Add the fourth order adjustment term to the 2nd order term
! and add the u-equation vertical advection term to the
! horizontal advection terms.
!
!-----------------------------------------------------------------------
DO k=3,nz-3
DO j=1,ny-1
DO i=2,nx-1
uforce(i,j,k) = uforce(i,j,k) - &
foursixth*(tem2(i,j,k)+tem2(i,j,k+1)) + &
onesix*(tem3(i,j,k+1)+tem3(i,j,k-1))
END DO
END DO
END DO
!call test_dump (uforce,"XXX1advu_uforce",nx,ny,nz,1,0)
!call test_dump (tem2,"XXX1advu_tem2",nx,ny,nz,1,0)
!call test_dump (tem3,"XXX1advu_tem3",nx,ny,nz,1,0)
DO k=2,nz-2,nz-4
DO j=1,ny-1
DO i=2,nx-1
uforce(i,j,k)=uforce(i,j,k)-(tem2(i,j,k+1)+tem2(i,j,k))*0.5
END DO
END DO
END DO
!call test_dump (uforce,"XXX2advu_uforce",nx,ny,nz,1,0)
ELSE ! perform second order advection....
DO k=2,nz-2
DO j=1,ny-1
DO i=2,nx-1
uforce(i,j,k)=uforce(i,j,k)-(tem2(i,j,k+1)+tem2(i,j,k))*0.5
END DO
END DO
END DO
!call test_dump (uforce,"XXX3advu_uforce",nx,ny,nz,1,0)
END IF
RETURN
END SUBROUTINE advu
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ADVV ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE advv(nx,ny,nz,u,v,w,ustr,vstr,wstr, vbar, mapfct, & 1,10
vadv, vforce, &
tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the advection terms of the v-equation. These terms are
! written in equivalent advection form.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 11/20/91.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (M. Xue)
! Added full documentation.
!
! 5/29/92 (K. Brewster)
! Further facelift.
!
! 4/9/93 (M. Xue & K. Brewster)
! Some index bounds for operators and loop bounds were corrected.
! Redundant calculations were done before at the boundaries that
! caused out-of-bound calculations, which should not have
! affected the results in normal situations though.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 2/9/94 (D. Jahn)
! Add 4th-order momentum advection.
!
! 9/1/94 (D. Weber & Y. Lu)
! Cleaned up documentation
!
! 11/5/1995 (M. Xue)
! Added vertical fourth order advection for v.
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 9/1/98 (D. Weber and T. Chung)
! Removed operators and merged loops 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
!
! mapfct Map factors at v points
!
! v y component of velocity at a given time level (m/s)
!
! ustr u * rhostr
! vstr v * rhostr
! wstr w * rhostr
!
! vbar Base state y component of velocity (m/s)
!
! OUTPUT:
!
! vforce Acoustically inactive forcing terms in v-momentum
! equation (kg/(m*s)**2).
!
! WORK ARRAYS:
!
! vadv v equation advection terms (kg/m**3)*(m/s)/s
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INCLUDE 'timelvls.inc'
REAL :: mapfct(nx,ny) ! Map factors at v points
REAL :: u (nx,ny,nz,nt) ! Total u-velocity at a given time level
! (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity at a given time level
! (m/s)
REAL :: w (nx,ny,nz,nt) ! Total w-velocity at a given time level
! (m/s)
REAL :: ustr (nx,ny,nz) ! u*rhostr
REAL :: vstr (nx,ny,nz) ! v*rhostr
REAL :: wstr (nx,ny,nz) ! w*rhostr
REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s)
REAL :: vadv (nx,ny,nz) ! Advection term of the v-eq.
! (kg/m**3)*(m/s)/s
REAL :: vforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in v-momentum equation (kg/(m*s)**2)
REAL :: tem1 (nx,ny,nz) ! Temporary work array
REAL :: tem2 (nx,ny,nz) ! Temporary work array
REAL :: tem3 (nx,ny,nz) ! Temporary work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,onvf
INTEGER :: ibgn,iend,jbgn,jend
REAL :: ueb, uwb
REAL :: fourthirds,foursixth,onesix
REAL :: tema
REAL :: temb
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global model control parameters
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'bndry.inc' ! Parameters for boundary conditions.
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
fourthirds = (4.0/3.0)
foursixth = (4.0/6.0)
onesix = (1.0/6.0)
!-----------------------------------------------------------------------
!
! Calculate rho*u* dv/dx = avgx( avgy( ustr ) * difx( v ) )
!
!-----------------------------------------------------------------------
tema = 0.25*dxinv
DO k=2,nz-1
DO j=2,ny-1
DO i=2,nx-2
vadv(i,j,k)=tema*((ustr(i+1,j,k)+ustr(i+1,j-1,k))* &
(v(i+1,j,k,2)-v(i,j,k,2)) &
+(ustr(i,j,k)+ustr(i,j-1,k))* &
(v(i,j,k,2)-v(i-1,j,k,2)))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! For 4th-order momentum advection, calculate the adjustment term
! avg2x( avgx(avgy(ustr))*dif2x(v) )
!
!-----------------------------------------------------------------------
IF (madvopt == 2 .OR. madvopt == 3) THEN
tema = 0.125 * dxinv
DO k=2,nz-2
DO j=2,ny-1
DO i=2,nx-2
tem2(i,j,k)=tema*((ustr(i,j,k)+ustr(i,j-1,k)) + &
(ustr(i+1,j,k)+ustr(i+1,j-1,k))) * &
(v(i+1,j,k,2)-v(i-1,j,k,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! In the case of periodic boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
! call test_dump (tem2,"XXX1advct_tem2",nx,ny,nz,0,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(tem2,nx,ny,nz,ebc,wbc,0,mptag,tem1)
CALL mprecv2dew(tem2,nx,ny,nz,ebc,wbc,0,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
ibgn = 3
iend = nx-3
jbgn = 2
jend = ny-1
IF (wbc == 0) ibgn = 2
IF (wbc == 2) THEN
ibgn = 2
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO j=jbgn,jend
tem2(1,j,k) = tem2(nx-2,j,k)
END DO
END DO
END IF
END IF
IF (ebc == 0) iend = nx-2
IF (ebc == 2) THEN
iend = nx-2
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO j=jbgn,jend
tem2(nx-1,j,k) = tem2(2,j,k)
END DO
END DO
END IF
END IF
!-----------------------------------------------------------------------
!
! In the case of rigid boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
IF (wbc == 1) THEN
ibgn = 2
DO k=2,nz-2
DO j=jbgn,jend
tem2(1,j,k) = tem2(2,j,k)
END DO
END DO
END IF
IF (ebc == 1) THEN
iend = nx-2
DO k=2,nz-2
DO j=jbgn,jend
tem2(nx-1,j,k) = tem2(nx-2,j,k)
END DO
END DO
END IF
! call test_dump (tem2,"XXX1Aadvct_tem2",nx,ny,nz,0,0) ! k=1,nz-1 not set
CALL acct_stop_inter
!-----------------------------------------------------------------------
!
! Concatenate the 2nd and 4th-order terms for urho*dv/dx
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=jbgn,jend
DO i=ibgn,iend
vadv(i,j,k) = fourthirds*vadv(i,j,k) - &
onesix*(tem2(i+1,j,k)+tem2(i-1,j,k))
END DO
END DO
END DO
END IF
!-----------------------------------------------------------------------
!
! Calculate rho*u* dv/dx at the east and west boundaries
! using one-sided advection.
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=2,ny-1
uwb=(ustr(1,j,k)+ustr(2,j,k)+ustr(1,j-1,k)+ustr(2,j-1,k))*0.25
ueb=(ustr(nx-1,j,k)+ustr(nx,j,k)+ustr(nx-1,j-1,k) &
+ustr(nx,j-1,k))*0.25
IF (uwb < 0.0) THEN
vadv(1,j,k)=uwb*(v(2,j,k,tpresent)-v(1,j,k,tpresent))*dxinv
ELSE
vadv(1,j,k)=rlxlbc*uwb*(v(1,j,k,tpast)-vbar(1,j,k))*dxinv
END IF
IF (ueb > 0.0) THEN
vadv(nx-1,j,k)=ueb* &
(v(nx-1,j,k,tpresent)-v(nx-2,j,k,tpresent))*dxinv
ELSE
vadv(nx-1,j,k)=-rlxlbc*ueb* &
(v(nx-1,j,k,tpast)-vbar(nx-1,j,k))*dxinv
END IF
END DO
END DO
!-----------------------------------------------------------------------
!
! Calculate rho*v* dv/dy = avgy( avgy( vstr ) * dify( v ) )
!
!-----------------------------------------------------------------------
tema = 0.25*dyinv
DO k=2,nz-2
DO j=2,ny-1
DO i=1,nx-1
tem2(i,j,k)=tema*((vstr(i,j+1,k)+vstr(i,j,k))* &
(v(i,j+1,k,2)-v(i,j,k,2)) &
+(vstr(i,j,k)+vstr(i,j-1,k))* &
(v(i,j,k,2)-v(i,j-1,k,2)))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! tem2 contains avgy( avgy( vstr ) * dify( v ) ) second order term
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! For 4th-order momentum advection, calculate the adjustment term:
! avg2y (avg2y(vstr) * dif2y(v) )
!
!-----------------------------------------------------------------------
IF (madvopt == 2 .OR. madvopt == 3) THEN
tema = 0.25*dyinv
DO k=2,nz-2
DO j=2,ny-1
DO i=1,nx-1
tem3(i,j,k)=tema*(vstr(i,j+1,k)+vstr(i,j-1,k))* &
(v(i,j+1,k,2)-v(i,j-1,k,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! At this point tem3 contains avg2y(vstr) * dif2y(v)
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! In the case of periodic boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
! call test_dump (tem3,"XXX1advct_tem3",nx,ny,nz,2,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dns(tem3,nx,ny,nz,nbc,sbc,2,mptag,tem1)
CALL mprecv2dns(tem3,nx,ny,nz,nbc,sbc,2,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
ibgn = 1
iend = nx-1
jbgn = 3
jend = ny-2
IF (nbc == 0) jend = ny-1
IF (nbc == 2) THEN
jend = ny-1
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,ny,k) = tem3(i,3,k)
END DO
END DO
END IF
END IF
IF (sbc == 0) jbgn = 2
IF (sbc == 2) THEN
jbgn = 2
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,1,k) = tem3(i,ny-2,k)
END DO
END DO
END IF
END IF
!-----------------------------------------------------------------------
!
! In the case of rigid boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
IF (nbc == 1) THEN
jend = ny-1
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,ny,k) = -tem3(i,ny-2,k)
END DO
END DO
END IF
IF (sbc == 1) THEN
jbgn = 2
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,1,k) = -tem3(i,3,k)
END DO
END DO
END IF
! call test_dump (tem3,"XXX1Aadvct_tem3",nx,ny,nz,2,0) ! k=1,nz-1 not set
CALL acct_stop_inter
!-----------------------------------------------------------------------
!
! Concatenate the 2nd- and 4th-order terms: vstr*dv/dy
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=jbgn,jend
DO i=ibgn,iend
tem2(i,j,k) = fourthirds*tem2(i,j,k) - &
onesix*(tem3(i,j+1,k)+tem3(i,j-1,k))
END DO
END DO
END DO
END IF
!-----------------------------------------------------------------------
!
! Sum the x and y contributions to v advection
! At this point tem1 contains the y contribution to the v advection
!
!-----------------------------------------------------------------------
!
! Add the v-equation horizontal advection terms, rho*u*dv/dx (vadv) and
! rho*v*dv/dy (tem1).
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=2,ny-1
DO i=1,nx-1
vforce(i,j,k)=vforce(i,j,k) - &
(vadv(i,j,k)+tem2(i,j,k)) * mapfct(i,j)
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Calculate rho*w* dv/dz = avgz( avgy( wstr ) * difz( v ) )
!
!-----------------------------------------------------------------------
tema = 0.5*dzinv
DO k=2,nz-1
DO j=2,ny-1
DO i=1,nx-1
tem2(i,j,k)=tema*(wstr(i,j,k)+wstr(i,j-1,k))* &
(v(i,j,k,2)-v(i,j,k-1,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! tem2 contains avgz( avgy( wstr ) * difz( v ) )
!
!-----------------------------------------------------------------------
IF (madvopt == 3) THEN
!-----------------------------------------------------------------------
!
! For 4th-order momentum advection, calculate the adjustment term.
! avg2z( avgz( avgy (wstr) )* dif2z( v ))
!
!-----------------------------------------------------------------------
tema = 0.125 * dzinv
DO k=2,nz-2
DO j=2,ny-1
DO i=1,nx-1
tem3(i,j,k) = tema*((wstr(i,j,k+1)+wstr(i,j-1,k+1))+ &
(wstr(i,j,k)+wstr(i,j-1,k))) * &
(v(i,j,k+1,2)-v(i,j,k-1,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Add the fourth order contribution to the 2nd order term and
! add the v-equation vertical advection term rho*w*dv/dz (tem1) to
! the horizontal advection terms.
!
!-----------------------------------------------------------------------
DO k=3,nz-3
DO j=2,ny-1
DO i=1,nx-1
vforce(i,j,k) = vforce(i,j,k) - &
foursixth*(tem2(i,j,k+1)+tem2(i,j,k)) + &
onesix*(tem3(i,j,k+1)+tem3(i,j,k-1))
END DO
END DO
END DO
DO k=2,nz-2,nz-4
DO j=2,ny-1
DO i=1,nx-1
vforce(i,j,k)=vforce(i,j,k) - (tem2(i,j,k+1)+tem2(i,j,k))*0.5
END DO
END DO
END DO
ELSE ! perform second order advection...
DO k=2,nz-2
DO j=2,ny-1
DO i=1,nx-1
vforce(i,j,k)=vforce(i,j,k) - (tem2(i,j,k+1)+tem2(i,j,k))*0.5
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE advv
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ADVW ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE advw(nx,ny,nz,u,v,w,ustr,vstr,wstr, mapfct, & 1,10
wadv, wforce, &
tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
!
! Calculate the advection terms of the w-equation. These terms are
! written in equivalent advection form.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 11/20/91.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (M. Xue)
! Added full documentation.
!
! 5/29/92 (K. Brewster)
! Further facelift.
!
! 4/9/93 (M. Xue & K. Brewster)
! Some index bounds for operators and loop bounds were corrected.
! Redundant calculations were done before at the boundaries that
! caused out-of-bound calculations, which should not have
! affected the results in normal situations though.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 2/9/94 (D. Jahn)
! Add 4th-order momentum advection.
!
! 9/1/94 (D. Weber & Y. Lu)
! Cleaned up documentation
!
! 11/5/1995 (M. Xue)
! Added vertical fourth order advection for w.
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 9/1/98 (D. Weber and T. Chung)
! Removed operators and merged loops 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
!
! mapfct Map factors at scalar points
!
! w z component of velocity at a given time level (m/s)
!
! ustr u * rhostr
! vstr v * rhostr
! wstr w * rhostr
!
! OUTPUT:
!
! wforce Acoustically inactive forcing terms in w-momentum
! equation (kg/(m*s)**2).
!
! WORK ARRAYS:
!
! wadv Advection term in w equation (kg/m**3)*(m/s)/s
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INCLUDE 'timelvls.inc'
REAL :: mapfct(nx,ny) ! Map factors at scalar points
REAL :: u (nx,ny,nz,nt) ! Total u-velocity at a given time level
! (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity at a given time level
! (m/s)
REAL :: w (nx,ny,nz,nt) ! Total w-velocity at a given time level
! (m/s)
REAL :: ustr (nx,ny,nz) ! u*rhostr
REAL :: vstr (nx,ny,nz) ! v*rhostr
REAL :: wstr (nx,ny,nz) ! w*rhostr
REAL :: wadv (nx,ny,nz) ! w-eqn advection term (kg/m**3)*(m/s)/s
REAL :: wforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in w-momentum equation (kg/(m*s)**2)
REAL :: tem1 (nx,ny,nz) ! Temporary work array
REAL :: tem2 (nx,ny,nz) ! Temporary work array
REAL :: tem3 (nx,ny,nz) ! Temporary work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,onvf
INTEGER :: ibgn,iend,jbgn,jend
REAL :: ueb, uwb, vsb, vnb
REAL :: fourthirds,foursixth,onesix
REAL :: tema
REAL :: temb
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global model control parameters
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'bndry.inc' ! Parameters for boundary conditions.
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
fourthirds = (4.0/3.0)
foursixth = (4.0/6.0)
onesix = (1.0/6.0)
!
!-----------------------------------------------------------------------
!
! Calculate rho*u* dw/dx = avgx( avgz( ustr ) * difx( w ) )
!
!-----------------------------------------------------------------------
tema = 0.25 * dxinv
DO k=2,nz-1
DO j=1,ny-1
DO i=2,nx-2
wadv(i,j,k)=tema*((ustr(i+1,j,k)+ustr(i+1,j,k-1))* &
(w(i+1,j,k,2)-w(i,j,k,2)) &
+(ustr(i,j,k)+ustr(i,j,k-1))* &
(w(i,j,k,2)-w(i-1,j,k,2)))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! At this point wadv contains avgx(avgz(u*rho)*dw/dx)
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! For 4th-order momentum advection, calculate the 4th-order
! adjustment term for ustr*dw/dx:
! avg2x( avgx(avgz(ustr)) * dif2x(w) )
!
!-----------------------------------------------------------------------
IF (madvopt == 2 .OR. madvopt == 3) THEN
tema = 0.125 * dxinv
DO k=2,nz-1
DO j=1,ny-1
DO i=2,nx-2
tem2(i,j,k)=tema*((ustr(i,j,k-1)+ustr(i,j,k)) + &
(ustr(i+1,j,k-1)+ustr(i+1,j,k))) * &
(w(i+1,j,k,2)-w(i-1,j,k,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! At this point tem2 contains avgxz(ustr)*dw/dx
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! In the case of periodic boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
! call test_dump (tem2,"XXX2advct_tem2",nx,ny,nz,0,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(tem2,nx,ny,nz,ebc,wbc,0,mptag,tem1)
CALL mprecv2dew(tem2,nx,ny,nz,ebc,wbc,0,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
ibgn = 3
iend = nx-3
jbgn = 1
jend = ny-1
IF (ebc == 0) iend = nx-2
IF (ebc == 2) THEN
iend = nx-2
IF (mp_opt == 0) THEN
DO k=2,nz-1
DO j=jbgn,jend
tem2(nx-1,j,k) = tem2(2,j,k)
END DO
END DO
END IF
END IF
IF (wbc == 0) ibgn = 2
IF (wbc == 2) THEN
ibgn = 2
IF (mp_opt == 0) THEN
DO k=2,nz-1
DO j=jbgn,jend
tem2(1,j,k) = tem2(nx-2,j,k)
END DO
END DO
END IF
END IF
!-----------------------------------------------------------------------
!
! In the case of rigid boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
IF (ebc == 1) THEN
iend = nx-2
DO k=2,nz-1
DO j=jbgn,jend
tem2(nx-1,j,k) = tem2(nx-2,j,k)
END DO
END DO
END IF
IF (wbc == 1) THEN
ibgn = 2
DO k=2,nz-1
DO j=jbgn,jend
tem2(1,j,k) = tem2(2,j,k)
END DO
END DO
END IF
! call test_dump (tem2,"XXX2Aadvct_tem2",nx,ny,nz,0,0) ! k=1,nz-1 not set
CALL acct_stop_inter
!-----------------------------------------------------------------------
!
! Concatenate the 2nd- and 4th-order terms.
!
!-----------------------------------------------------------------------
DO k=2,nz-1
DO j=jbgn,jend
DO i=ibgn,iend
wadv(i,j,k) = fourthirds*wadv(i,j,k)- &
onesix*(tem2(i+1,j,k)+tem2(i-1,j,k))
END DO
END DO
END DO
END IF
!-----------------------------------------------------------------------
!
! Calculate rho*u* dw/dx on the east and west boundaries
! using one-sided advection.
!
!-----------------------------------------------------------------------
DO k=2,nz-1
DO j=1,ny-1
uwb=(ustr(1,j,k)+ustr(2,j,k)+ustr(1,j,k-1)+ustr(2,j,k-1))*0.25
ueb=(ustr(nx-1,j,k)+ustr(nx,j,k)+ustr(nx-1,j,k-1) &
+ustr(nx,j,k-1))*0.25
IF (uwb < 0.0) THEN
wadv(1,j,k)=uwb*(w(2,j,k,tpresent)-w(1,j,k,tpresent))*dxinv
ELSE
wadv(1,j,k)=rlxlbc*uwb*w(1,j,k,tpast)*dxinv
END IF
IF (ueb > 0.0) THEN
wadv(nx-1,j,k)=ueb* &
(w(nx-1,j,k,tpresent)-w(nx-2,j,k,tpresent))*dxinv
ELSE
wadv(nx-1,j,k)=-rlxlbc*ueb*w(nx-1,j,k,tpast)*dxinv
END IF
END DO
END DO
!-----------------------------------------------------------------------
!
! Calculate rho*v* dw/dy = avgy( avgz( vstr ) * dify ( w ) )
!
!-----------------------------------------------------------------------
tema=0.25*dyinv
DO k=2,nz-1
DO j=2,ny-2
DO i=1,nx-1
tem2(i,j,k)=tema*((vstr(i,j+1,k)+vstr(i,j+1,k-1)) * &
(w(i,j+1,k,2)-w(i,j,k,2)) &
+(vstr(i,j,k)+vstr(i,j,k-1)) * &
(w(i,j,k,2)-w(i,j-1,k,2)))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! At this point tem2 contains avgz(v*rho)*dw/dy second order term
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! For 4th-order momentum advection, calculate the adjustment
! term for vstr*dw/dy:
! avg2y( avgz(avgy(vstr)) * dif2y(w) )
!
!-----------------------------------------------------------------------
IF (madvopt == 2 .OR. madvopt == 3) THEN
tema = 0.125*dyinv
DO k=2,nz-1
DO j=2,ny-2
DO i=1,nx-1
tem3(i,j,k)=tema*((vstr(i,j+1,k)+vstr(i,j,k)) + &
(vstr(i,j+1,k-1)+vstr(i,j,k-1))) * &
(w(i,j+1,k,2)-w(i,j-1,k,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! At this point tem3 contains avgyz(vstr)*dw/dy
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! In the case of periodic boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
! call test_dump (tem3,"XXX2advct_tem3",nx,ny,nz,0,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dns(tem3,nx,ny,nz,nbc,sbc,0,mptag,tem1)
CALL mprecv2dns(tem3,nx,ny,nz,nbc,sbc,0,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
ibgn = 1
iend = nx-1
jbgn = 3
jend = ny-3
IF (nbc == 0) jend = ny-2
IF (nbc == 2) THEN
jend = ny-2
IF (mp_opt == 0) THEN
DO k=2,nz-1
DO i=ibgn,iend
tem3(i,ny-1,k) = tem3(i,2,k)
END DO
END DO
END IF
END IF
IF (sbc == 0) jbgn = 2
IF (sbc == 2) THEN
jbgn = 2
IF (mp_opt == 0) THEN
DO k=2,nz-1
DO i=ibgn,iend
tem3(i,1,k)=tem3(i,ny-2,k)
END DO
END DO
END IF
END IF
!-----------------------------------------------------------------------
!
! In the case of rigid boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
IF (nbc == 1) THEN
jend = ny-2
DO k=2,nz-1
DO i=ibgn,iend
tem3(i,ny-1,k) = tem3(i,ny-2,k)
END DO
END DO
END IF
IF (sbc == 1) THEN
jbgn = 2
DO k=2,nz-1
DO i=ibgn,iend
tem3(i,1,k)=tem3(i,2,k)
END DO
END DO
END IF
! call test_dump (tem3,"XXX2Aadvct_tem3",nx,ny,nz,0,0) ! k=1,nz-1 not set
CALL acct_stop_inter
!-----------------------------------------------------------------------
!
! Concatenate the 2nd- and 4th- order terms.
!
!-----------------------------------------------------------------------
DO k=2,nz-1
DO j=jbgn,jend
DO i=ibgn,iend
tem2(i,j,k) = fourthirds*tem2(i,j,k) - &
onesix*(tem3(i,j+1,k)+tem3(i,j-1,k))
END DO
END DO
END DO
END IF
!-----------------------------------------------------------------------
!
! Calculate rho*v* dw/dy on the north and south boundaries
! using one-sided advection.
!
!-----------------------------------------------------------------------
DO k=2,nz-1
DO i=1,nx-1
vsb=(vstr(i,1,k)+vstr(i,2,k)+vstr(i,1,k-1)+vstr(i,2,k-1))*0.25
vnb=(vstr(i,ny-1,k)+vstr(i,ny,k)+vstr(i,ny-1,k-1) &
+vstr(i,ny,k-1))*0.25
IF (vsb < 0.0) THEN
tem2(i,1,k)=vsb*(w(i,2,k,tpresent)-w(i,1,k,tpresent))*dyinv
ELSE
tem2(i,1,k)=rlxlbc*vsb*w(i,1,k,tpast)*dyinv
END IF
IF (vnb > 0.0) THEN
tem2(i,ny-1,k)=vnb* &
(w(i,ny-1,k,tpresent)-w(i,ny-2,k,tpresent))*dyinv
ELSE
tem2(i,ny-1,k)=-rlxlbc*vnb*w(i,ny-1,k,tpast)*dyinv
END IF
END DO
END DO
!-----------------------------------------------------------------------
!
! Add the w-equation horizontal advection terms, rho*u*dw/dx (wadv) and
! rho*v*dw/dy (tem1).
!
!-----------------------------------------------------------------------
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
wforce(i,j,k)=wforce(i,j,k) - &
(wadv(i,j,k)+tem2(i,j,k)) * mapfct(i,j)
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Calculate rho*w* dw/dz = avgz( avgz( wstr ) * difz( w ) )
!
!-----------------------------------------------------------------------
tema = 0.5*dzinv
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=tema*(wstr(i,j,k+1)+wstr(i,j,k))* &
(w(i,j,k+1,2)-w(i,j,k,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! tem2 contains avgz(w*rho)*dw/dz
!
!-----------------------------------------------------------------------
IF (madvopt == 3) THEN
!-----------------------------------------------------------------------
!
! If momentum advection 4th-order, calculate adjustment term:
! avg2z(avg2z(wstr)*dif2z(w))
!
!-----------------------------------------------------------------------
tema = 0.25*dzinv
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)=tema*(wstr(i,j,k+1)+wstr(i,j,k-1))* &
(w(i,j,k+1,2)-w(i,j,k-1,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Add the fourth order contribution to the 2nd order term and
! add the w-equation vertical advection term rho*w*dw/dz (tem1)
! to the horizontal advection terms.
!
!-----------------------------------------------------------------------
DO k=3,nz-2
DO j=1,ny-1
DO i=1,nx-1
wforce(i,j,k) = wforce(i,j,k) - &
foursixth*(tem1(i,j,k-1)+tem1(i,j,k)) + &
onesix*(tem2(i,j,k+1)+tem2(i,j,k-1))
END DO
END DO
END DO
DO k=2,nz-1,nz-3
DO j=1,ny-1
DO i=1,nx-1
wforce(i,j,k)=wforce(i,j,k) - (tem1(i,j,k)+tem1(i,j,k-1))*0.5
END DO
END DO
END DO
ELSE ! perform second order advection....
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
wforce(i,j,k)=wforce(i,j,k) - (tem1(i,j,k)+tem1(i,j,k-1))*0.5
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE advw
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ADVPT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE advpt(nx,ny,nz,ptprt,u,v,w,wcont, & 1,2
rhostr,ptbar, mapfct,j3,j3inv, &
ptadv, &
ustr,vstr,wstr,tem1,tem2,tem3,tem4,mp_tem)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the advection of total potential temperature
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 11/20/91.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (M. Xue)
! Added full documentation.
!
! 5/29/92 (K. Brewster)
! Further facelift and order of argument list changed.
!
! 4/10/93 (M. Xue & Hao Jin)
! Added the terrain.
!
! 5/26/93 (M. Xue)
! Corrected ptbar advection term and added w and j3 into the
! argument list.
!
! 9/4/94 (D. Jahn)
! Add tem4 to the call of advcts.
!
! 9/1/94 (D. Weber & Y. Lu)
! Cleaned up documentation
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 3/8/96 (Ming Xue)
! Fixed a bug in DO 510 i=1,nx-1, which was written as DO 510 i=i,nx-1.
! Effectively, the vertical advection remained 2nd order with the error.
!
! 9/1/98 (D. Weber and T. Chung)
! Removed operators and merged loops 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
!
! ptprt Perturbation potential temperature at all time levels
! (K)
! u x component of velocity at all time levels (m/s)
! v y component of velocity at all time levels (m/s)
! w Vertical component of velocity in Cartesian
! coordinates at all time levels (m/s).
! wcont Contravariant vertical velocity (m/s)
! rhostr Base state air density (kg/m**3)
! ptbar Base state potential temperature (K)
!
! mapfct Map factors at scalar points
!
! j3 Coordinate transformation Jacobian d(zp)/dz
!
! OUTPUT:
!
! ptadv Advection term of potential temperature eqn
! (kg/m**3)*K/s
!
! WORK ARRAYS:
!
! ustr u * rhostr
! vstr v * rhostr
! wstr w * rhostr
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
! tem4 Temporary work array.
! mp_tem Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INCLUDE 'timelvls.inc'
REAL :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: u (nx,ny,nz,nt) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nt) ! Total v-velocity (m/s)
REAL :: w (nx,ny,nz,nt) ! Total w-velocity (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: rhostr(nx,ny,nz) ! Base state air density (kg/m**3)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: mapfct(nx,ny) ! Map factors at scalar points
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian
! d(zp)/d(z)
REAL :: j3inv (nx,ny,nz) ! Coordinate transformation Jacobian
! d(zp)/d(z)
REAL :: ptadv(nx,ny,nz) ! Potential temperature advection
REAL :: ustr (nx,ny,nz) ! u * rhostr
REAL :: vstr (nx,ny,nz) ! v * rhostr
REAL :: wstr (nx,ny,nz) ! w * rhostr
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 :: mp_tem(MAX(nx,ny)*nz) ! Temporary message passing array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,onvf,tlevel
REAL :: tema
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global model control parameters
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'bndry.inc' ! Parameters for boundary conditions.
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
tlevel=tpresent
!
!-----------------------------------------------------------------------
!
! Calculate ustr=rhostr*u, vstr=rhostr*v, wstr=rhostr*w
!
!-----------------------------------------------------------------------
CALL rhouvw(nx,ny,nz,rhostr,tem1,tem2,tem3)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx
ustr(i,j,k)=u(i,j,k,tlevel)*tem1(i,j,k)
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny
DO i=1,nx-1
vstr(i,j,k)=v(i,j,k,tlevel)*tem2(i,j,k)
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
wstr(i,j,k)=wcont(i,j,k)*tem3(i,j,k)
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Perturbation potential temperature advection
!
!-----------------------------------------------------------------------
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem3(i,j,k) = 0.0
END DO
END DO
END DO
CALL advcts(nx,ny,nz, ptprt, &
u(1,1,1,tlevel),v(1,1,1,tlevel), &
ustr,vstr,wstr,tem3, mapfct, &
ptadv, tem1,tem2,tem4,mp_tem)
IF( ptsmlstp == 0 ) THEN
!-----------------------------------------------------------------------
!
! 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.
!
!-----------------------------------------------------------------------
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k) = rhostr(i,j,k)*j3inv(i,j,k)
END DO
END DO
END DO
tema=0.5*dzinv
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=tema*(tem2(i,j,k)+tem2(i,j,k-1))*w(i,j,k,2)* &
(ptbar(i,j,k)-ptbar(i,j,k-1))
END DO
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
ptadv(i,j,k)=ptadv(i,j,k)+0.5*(tem1(i,j,k+1)+tem1(i,j,k))
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE advpt
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ADVCTS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE advcts(nx,ny,nz,s,u,v,ustr,vstr,wstr,sbar, mapfct, & 3,10
sadv, &
tem1,tem2,tem3,mp_tem)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the scalar equation advection terms. These terms are
! written in equivalent advection form.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 11/20/91.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (M. Xue)
! Added full documentation.
!
! 4/9/93 (M. Xue & K. Brewster)
! Some index bounds for operators and loop bounds were corrected.
! Redundant calculations were done before at the boundaries that
! caused out-of-bound calculations, which should not have
! affected the results in normal situations though.
!
! 5/29/92 (K. Brewster)
! Further facelift.
!
! 2/9/94 (D. Jahn)
! Add 4th-order momentum advection.
!
! 9/1/94 (D. Weber & Y. Lu)
! Cleaned up documentation
!
! 11/5/1995 (M. Xue)
! Added vertical fourth order advection for scalars.
!
! 1/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 3/8/96 (Ming Xue)
! Fixed a bug in DO 510 i=1,nx-1, which was written as DO 510 i=i,nx-1.
! Effectively, the vertical advection remained 2nd order with the error.
!
! 9/1/98 (D. Weber and T. Chung)
! Removed operators and merged loops 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
!
! mapfct Map factors at scalar points
!
! s scalar field at a given time level (scalar units)
!
! ustr u * rhostr
! vstr v * rhostr
! wstr w * rhostr
!
! OUTPUT:
!
! sadv Advection term in scalar equation (kg/m**3)*scalar
! unit/s
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
!
! mp_tem Temporary work array.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INCLUDE 'timelvls.inc'
REAL :: mapfct(nx,ny) ! Map factors at scalar points
REAL :: s (nx,ny,nz,nt) ! a scalar field to be advected.
REAL :: u (nx,ny,nz) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz) ! Total v-velocity (m/s)
REAL :: ustr (nx,ny,nz) ! u * rhostr
REAL :: vstr (nx,ny,nz) ! v * rhostr
REAL :: wstr (nx,ny,nz) ! w * rhostr
REAL :: sbar (nx,ny,nz) ! A state of 's' towards which 's' is
! relaxed on the inflow boundaries.
REAL :: sadv (nx,ny,nz) ! advection term in scalar equation
! (kg/m**3)*(scalar units)/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 :: mp_tem(MAX(nx,ny)*nz) ! Temporary message passing array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,onvf
INTEGER :: ibgn,iend,jbgn,jend
REAL :: ueb, uwb, vsb, vnb
REAL :: fourthirds,foursixth,onesix
REAL :: tema
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Physical constants
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'bndry.inc' ! Parameters for boundary conditions.
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
fourthirds = (4.0/3.0)
foursixth = (4.0/6.0)
onesix = (1.0/6.0)
!-----------------------------------------------------------------------
!
! Calculate rho*u* ds/dx = avgx( ustr * difx ( s ) )
!
!-----------------------------------------------------------------------
tema = dxinv*0.5
DO k=2,nz-2
DO j=1,ny-1
DO i=2,nx-2
sadv(i,j,k)=tema*((s(i,j,k,2)-s(i-1,j,k,2))*ustr(i,j,k) &
+(s(i+1,j,k,2)-s(i,j,k,2))*ustr(i+1,j,k))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! For 4th-order scalar advection, calculate the 4th-order
! adjustment term for ustr*ds/dx:
! avg2x( avgx(ustr) * dif2x(s) )
!
!-----------------------------------------------------------------------
IF (sadvopt == 2.OR.sadvopt == 3) THEN
tema = 0.25*dxinv
DO k=2,nz-2
DO j=1,ny-1
DO i=2,nx-2
tem2(i,j,k)=tema*(ustr(i+1,j,k)+ustr(i,j,k))* &
(s(i+1,j,k,2)-s(i-1,j,k,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! In the case of periodic boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
! call test_dump (tem2,"XXX3advct_tem2",nx,ny,nz,0,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(tem2,nx,ny,nz,ebc,wbc,0,mptag,tem1)
CALL mprecv2dew(tem2,nx,ny,nz,ebc,wbc,0,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
ibgn = 3
iend = nx-3
jbgn = 1
jend = ny-1
IF (ebc == 0) iend = nx-2
IF (ebc == 2) THEN
iend = nx-2
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO j=jbgn,jend
tem2(nx-1,j,k) = tem2(2,j,k)
END DO
END DO
END IF
END IF
IF (wbc == 0) ibgn = 2
IF (wbc == 2) THEN
ibgn = 2
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO j=jbgn,jend
tem2(1,j,k) = tem2(nx-2,j,k)
END DO
END DO
END IF
END IF
!-----------------------------------------------------------------------
!
! In the case of rigid boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
IF (ebc == 1) THEN
iend = nx-2
DO k=2,nz-2
DO j=jbgn,jend
tem2(nx-1,j,k) = tem2(nx-2,j,k)
END DO
END DO
END IF
IF (wbc == 1) THEN
ibgn = 2
DO k=2,nz-2
DO j=jbgn,jend
tem2(1,j,k) = tem2(2,j,k)
END DO
END DO
END IF
! call test_dump (tem2,"XXX3Aadvct_tem2",nx,ny,nz,0,0) ! k=1,nz-1 not set
CALL acct_stop_inter
!-----------------------------------------------------------------------
!
! Concatenate the 2nd- and 4th- order terms for ustr*ds/dx
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=jbgn,jend
DO i=ibgn,iend
sadv(i,j,k) = fourthirds*sadv(i,j,k) - &
onesix*(tem2(i+1,j,k)+tem2(i-1,j,k))
END DO
END DO
END DO
END IF
!-----------------------------------------------------------------------
!
! Calculate rho*u* ds/dx on the east and west boundaries
! using one-sided advection.
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=1,ny-1
uwb=(ustr(1,j,k)+ustr(2,j,k))*0.5
ueb=(ustr(nx-1,j,k)+ustr(nx,j,k))*0.5
IF (uwb < 0.0) THEN
sadv(1,j,k)=uwb*(s(2,j,k,tpresent)-s(1,j,k,tpresent))*dxinv
ELSE
sadv(1,j,k)=rlxlbc*uwb*(s(1,j,k,tpast)-sbar(1,j,k))*dxinv
END IF
IF (ueb > 0.0) THEN
sadv(nx-1,j,k)=ueb* &
(s(nx-1,j,k,tpresent)-s(nx-2,j,k,tpresent))*dxinv
ELSE
sadv(nx-1,j,k)=-rlxlbc*ueb* &
(s(nx-1,j,k,tpast)-sbar(nx-1,j,k))*dxinv
END IF
END DO
END DO
!-----------------------------------------------------------------------
!
! Calculate rho*v* ds/dy = avgy( vstr * dify( s ) )
!
!-----------------------------------------------------------------------
tema = 0.5*dyinv
DO k=2,nz-2
DO j=2,ny-2
DO i=1,nx-1
tem1(i,j,k)=tema*((s(i,j,k,2)-s(i,j-1,k,2))*vstr(i,j,k) &
+(s(i,j+1,k,2)-s(i,j,k,2))*vstr(i,j+1,k))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! For 4th-order scalar advection, calculate the 4th-order
! adjustment term for vstr*ds/dy:
! avg2y( avgy(vstr) * dif2y(s) )
!
!-----------------------------------------------------------------------
IF (sadvopt == 2.OR.sadvopt == 3) THEN
tema = 0.25*dyinv
DO k=2,nz-2
DO j=2,ny-2
DO i=1,nx-1
tem3(i,j,k)=tema*(vstr(i,j+1,k)+vstr(i,j,k))* &
(s(i,j+1,k,2)-s(i,j-1,k,2))
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! In the case of periodic boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
! call test_dump (tem3,"XXX3advct_tem3",nx,ny,nz,0,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dns(tem3,nx,ny,nz,nbc,sbc,0,mptag,mp_tem)
CALL mprecv2dns(tem3,nx,ny,nz,nbc,sbc,0,mptag,mp_tem)
END IF
CALL acct_interrupt(bc_acct)
ibgn = 1
iend = nx-1
jbgn = 3
jend = ny-3
IF (nbc == 0) jend = ny-2
IF (nbc == 2) THEN
jend = ny-2
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,ny-1,k) = tem3(i,2,k)
END DO
END DO
END IF
END IF
IF (sbc == 0) jbgn = 2
IF (sbc == 2) THEN
jbgn = 2
IF (mp_opt == 0) THEN
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,1,k) = tem3(i,ny-2,k)
END DO
END DO
END IF
END IF
!-----------------------------------------------------------------------
!
! In the case of rigid boundary condition, need to calculate
! the 4th-order term for one additional point near boundary.
!
!-----------------------------------------------------------------------
IF (nbc == 1) THEN
jend = ny-2
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,ny-1,k) = tem3(i,ny-2,k)
END DO
END DO
END IF
IF (sbc == 1) THEN
jbgn = 2
DO k=2,nz-2
DO i=ibgn,iend
tem3(i,1,k) = tem3(i,2,k)
END DO
END DO
END IF
! call test_dump (tem3,"XXX3Aadvct_tem3",nx,ny,nz,0,0) ! k=1,nz-1 not set
CALL acct_stop_inter
!-----------------------------------------------------------------------
!
! Concatenate the 2nd- and 4th-order terms for vstr*ds/dy.
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=jbgn,jend
DO i=ibgn,iend
tem1(i,j,k) = fourthirds*tem1(i,j,k) - &
onesix*(tem3(i,j+1,k)+tem3(i,j-1,k))
END DO
END DO
END DO
END IF
!-----------------------------------------------------------------------
!
! Calculate rho*v* ds/dy on the north and south boundaries
! using one-sided advection.
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO i=1,nx-1
vsb=(vstr(i,1,k)+vstr(i,2,k))*0.5
vnb=(vstr(i,ny-1,k)+vstr(i,ny,k))*0.5
IF (vsb < 0.0) THEN
tem1(i,1,k)=vsb*(s(i,2,k,tpresent)-s(i,1,k,tpresent))*dyinv
ELSE
tem1(i,1,k)=rlxlbc*vsb*(s(i,1,k,tpast)-sbar(i,1,k))*dyinv
END IF
IF (vnb > 0.0) THEN
tem1(i,ny-1,k)=vnb* &
(s(i,ny-1,k,tpresent)-s(i,ny-2,k,tpresent))*dyinv
ELSE
tem1(i,ny-1,k)=-rlxlbc*vnb* &
(s(i,ny-1,k,tpast)-sbar(i,ny-1,k))*dyinv
END IF
END DO
END DO
!-----------------------------------------------------------------------
!
! Add the scalar equation horizontal advection terms, rho*u*ds/dx (sadv)
! and rho*v*ds/dy (tem2).
!
!-----------------------------------------------------------------------
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
sadv(i,j,k)=(sadv(i,j,k)+tem1(i,j,k)) * mapfct(i,j)
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Calculate rho*w* ds/dz = avgz( wstr * difz( s ) )
!
!-----------------------------------------------------------------------
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=dzinv*(s(i,j,k,2)-s(i,j,k-1,2))*wstr(i,j,k)
END DO
END DO
END DO
IF (sadvopt == 3) THEN
!-----------------------------------------------------------------------
!
! For 4th-order scalar advection, calculate the 4th-order
! adjustment term for wstr*ds/dz: avg2z( avgz(wstr) * dif2z(s) )
! and combine it with the second order term.
!
! Note that we apply the 4th order scheme only upto the second point
! from the top and bottom boundary. We do not give special treatment
! to the periodic or wall boundary cases as we do for the side
! boundaries.
!
! Add the scalar equation vertical advection term rho*w*ds/dz (tem2)
! to the horizontal advection terms.
!
!-----------------------------------------------------------------------
tema = 0.25*dzinv
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem3(i,j,k)=tema*(wstr(i,j,k+1)+wstr(i,j,k))* &
(s(i,j,k+1,2)-s(i,j,k-1,2))
END DO
END DO
END DO
DO k=3,nz-3
DO j=1,ny-1
DO i=1,nx-1
sadv(i,j,k) = sadv(i,j,k)+ &
foursixth*(tem1(i,j,k+1)+tem1(i,j,k)) - &
onesix*(tem3(i,j,k+1)+tem3(i,j,k-1))
END DO
END DO
END DO
DO k=2,nz-2,nz-4
DO j=1,ny-1
DO i=1,nx-1
sadv(i,j,k)= sadv(i,j,k) + (tem1(i,j,k+1)+tem1(i,j,k))*0.5
END DO
END DO
END DO
ELSE ! perform second order advection...
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
sadv(i,j,k)= sadv(i,j,k) + (tem1(i,j,k+1)+tem1(i,j,k))*0.5
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE advcts