!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GRDTRAN ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE grdtran(nx,ny,nz,ubar,vbar,u,v,w,ptprt,pprt, & 2,3
qv,qc,qr,qi,qs,qh,qvbar,rhostr,x,y,zp,j3,j3inv, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
tem8,tem9,tem10,tem11)
!
!--------------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinate the calculation and adjustment of domain translation
! speed using either cell-tracking or optimal pattern translation
! algorithm.
!
!---------------------------------------------------------------------------
!
! AUTHORS: Ming Xue.
! 3/29/1995.
!
! MODIFICATION HISTORY:
!
! 04/18/1995 (Y. Liu)
! Set the default values for variable uchange and vchange.
!
! 06/24/1995 (Alan Shapiro)
! Documentation clean-up.
!
! 06/27/95 (M. Xue)
! Added qvbar to argument list of GRDTRAN and ADJUVWV.
!
!-----------------------------------------------------------------------
!
! 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 z-direction (vertical)
!
! ubar Model's base state x velocity component (m/s)
! vbar Model's base state y velocity component (m/s)
!
! u Total u-velocity (m/s)
! v Total v-velocity (m/s)
! w Total w-velocity (m/s)
! ptprt Perturbation potential temperature (K)
! pprt Perturbation pressure (Pascal)
! qv Water vapor specific humidity (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rain water mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
!
! qvbar Base-state water vapor mixing ratio (kg/kg)
! rhostr Base state density rhobar times j3 (kg/m**3)
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! zp Vertical coordinate of grid points in physical space (m)
! j3 Coordinate transformation Jacobian d(zp)/dz
!
! OUTPUT:
!
! ubar Redefined base state x velocity component (m/s)
! vbar Redefined base state y velocity component (m/s)
!
! u Redefined total u-velocity (m/s)
! v Redefined total v-velocity (m/s)
! w Redefined total w-velocity (m/s)
! ptprt Redefined perturbation potential temperature at tpast (K)
! pprt Redefined perturbation pressure at tpast (Pascal)
! qv Redefined water vapor specific humidity at tpast (kg/kg)
! qc Redefined cloud water mixing ratio at tpast (kg/kg)
! qr Redefined rain water mixing ratio at tpast (kg/kg)
! qi Redefined cloud ice mixing ratio at tpast (kg/kg)
! qs Redefined snow mixing ratio at tpast (kg/kg)
! qh Redefined hail mixing ratio at tpast (kg/kg)
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
! tem4 Temporary work array.
! tem5 Temporary work array.
! tem6 Temporary work array.
! tem7 Temporary work array.
! tem8 Temporary work array.
! tem9 Temporary work array.
! tem10 Temporary work array.
! tem11 Temporary work array.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'timelvls.inc'
INTEGER :: nx,ny,nz ! Number of grid points in x, y, z directions
REAL :: ubar (nx,ny,nz) ! Base state x velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state y velocity (m/s)
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 :: ptprt(nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz,nt) ! Perturbation pressure (Pascal)
REAL :: qv (nx,ny,nz,nt) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz,nt) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz,nt) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz,nt) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz,nt) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz,nt) ! Hail mixing ratio (kg/kg)
REAL :: qvbar(nx,ny,nz) ! Base-state water vapor mixing ratio (kg/kg)
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: x (nx) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL :: y (ny) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at
! w-point of staggered grid.
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ).
REAL :: j3inv (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ).
REAL :: tem1 (nx,ny,nz) ! Temporary work array.
REAL :: tem2 (nx,ny,nz) ! Temporary work array.
REAL :: tem3 (nx,ny,nz) ! Temporary work array.
REAL :: tem4 (nx,ny,nz) ! Temporary work array.
REAL :: tem5 (nx,ny,nz) ! Temporary work array.
REAL :: tem6 (nx,ny,nz) ! Temporary work array.
REAL :: tem7 (nx,ny,nz) ! Temporary work array.
REAL :: tem8 (nx,ny,nz) ! Temporary work array.
REAL :: tem9 (nx,ny,nz) ! Temporary work array.
REAL :: tem10(nx,ny,nz) ! Temporary work array.
REAL :: tem11(nx,ny,nz) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global constants that control model
! execution
INCLUDE 'grid.inc' ! Grid parameters
INTEGER :: i,j,k,ntrcell,ireturn
REAL :: uvmvmag, scal ,ucelcntr,vcelcntr,xcelcntr,ycelcntr
REAL :: uchange,vchange,xcctr_old,ycctr_old
INTEGER :: called
SAVE called, xcelcntr, ycelcntr
DATA called /0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
uchange = 0.0
vchange = 0.0
IF(MOD(nstep,nceltrk) == 0 .AND.(cltkopt == 1.OR.grdtrns == 2))THEN
!
!-----------------------------------------------------------------------
!
! Current model grid origin in absolute coord.
!
!-----------------------------------------------------------------------
!
xgrdorg = xgrdorg + tceltrk * umove
ygrdorg = ygrdorg + tceltrk * vmove
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem11(i,j,k) = rhostr(i,j,k)*j3inv(i,j,k)
END DO
END DO
END DO
IF( called /= 0 ) THEN
xcctr_old = xcelcntr
ycctr_old = ycelcntr
END IF
!
!-----------------------------------------------------------------------
!
! Find new cell center.
!
!-----------------------------------------------------------------------
!
CALL celtrk
(nx,ny,nz,2,nx-2,2,ny-2,2,nz-2, &
w(1,1,1,tpresent), qr(1,1,1,tpresent), &
tem11, x,y,zp, ntrcell, xcelcntr, ycelcntr, &
ireturn, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9,tem10)
!
!-----------------------------------------------------------------------
!
! Cell movement speed between the past two celtrk calls.
!
!-----------------------------------------------------------------------
!
IF( ntrcell > 0 ) THEN
IF( called /= 0 ) THEN
ucelcntr = (xcelcntr-xcctr_old)/tceltrk
vcelcntr = (ycelcntr-ycctr_old)/tceltrk
ELSE
ucelcntr = 0.0
vcelcntr = 0.0
END IF
called = 1
END IF
END IF
IF (grdtrns == 2 .AND. MOD(nstep,nceltrk) == 0) THEN
!
!-----------------------------------------------------------------------
!
! Calculate adjustment to umove and vmove (uchange,vchange) which
! tries to bring the center of mass of cells back
! to the center of model domain during time period tcrestr
! assuming the cell center moves at past speed.
!
!-----------------------------------------------------------------------
!
IF( ireturn == 0 .AND. called /= 0 ) THEN
uchange = ( xcelcntr+ucelcntr*tcrestr &
-((x(1)+x(nx))*0.5+xgrdorg+umove*tcrestr))/tcrestr
vchange = ( ycelcntr+vcelcntr*tcrestr &
-((y(1)+y(ny))*0.5+ygrdorg+vmove*tcrestr))/tcrestr
uvmvmag = SQRT( uchange*uchange + vchange*vchange )
IF( uvmvmag > 10.0 ) THEN
scal = 10.0 /uvmvmag
uchange = uchange * scal
vchange = vchange * scal
END IF
ELSE
uchange = 0.0
vchange = 0.0
END IF
ELSE IF (grdtrns == 3) THEN
!
!-----------------------------------------------------------------------
!
! Automatic domain translation.
!
! Call AUTOTRANS to compute running mean of the optimal pattern
! translation umove and vmove, and return the appropriate increments
! (uchange, vchange) in umove and vmove at the end of time window.
!
! The average domain translation speed in a given time window
! (twindow) is calculated from the movement speed of traced properties
! during each time step (autotrans is called every time step).
! Only at the end of this time window are non-zero values of umove
! and vmove sent out. At other times, they are zero, indicating
! no change should be made to umove or vmove.
!
!-----------------------------------------------------------------------
!
CALL autotrans(nx,ny,nz, &
ubar,vbar,u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
uchange, vchange, tem1)
END IF
!
!-----------------------------------------------------------------------
!
! Reset the model u and v velocity values using the new domain
! translation speed.
! (ADJUVMV should only be called when it's time to use it.)
!
!-----------------------------------------------------------------------
!
IF( uchange /= 0. .OR. vchange /= 0.) THEN
WRITE(6,'(3(/1x,a)/)') &
'ATTENTION: UMOVE or VMOVE has just been changed. ', &
'Subroutine ADJUVMV will now be called to adjust the time-', &
'dependent variables.'
CALL adjuvmv(nx,ny,nz, &
ubar,vbar,u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,qvbar, &
uchange, vchange, tem1, tem11)
END IF
RETURN
END SUBROUTINE grdtran
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ADJUVMV ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE adjuvmv(nx,ny,nz, & 2,99
ubar,vbar,u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,qvbar, &
uchange,vchange , tem1, tem2)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Adjust the model variables per change in the domain translation
! speed.
!
!-----------------------------------------------------------------------
!
! AUTHORS: Ming Xue
! Based on Alan Shapiro and Steve Lazarus' GRIDTRANS.
! 9/7/1993
!
! MODIFICATION HISTORY:
!
! 06/27/95 (M. Xue)
! Added qvbar to the argument list. bcqv called for qv.
!
! 11/20/97 (fanyou Kong - CMRP)
! Add a do loop to remove possible negative Qx values after
! calling galilei subroutine
!
! 7/16/98 (M. Xue and Y. Richardson)
! Fixed a minor bug in call to bcp. It was introduced
! when subroutine BCP was modified.
!
! 2000/04/21 (Gene Bassett)
! Adjusted BC calls so that tfuture is update for all variables (rather
! than u, v, w & pprt for tpast).
!
!-----------------------------------------------------------------------
!
! 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 z-direction (vertical)
!
! ubar Model's base state x velocity component (m/s)
! vbar Model's base state y velocity component (m/s)
!
! u Total u-velocity (m/s)
! v Total v-velocity (m/s)
! w Total w-velocity (m/s)
! ptprt Perturbation potential temperature (K)
! pprt Perturbation pressure (Pascal)
! qv Water vapor specific humidity (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rain water mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
! qvbar Base-state water vapor mixing ratio (kg/kg)
!
! uchange Change in the domain translation speed in x-dir.
! vchange Change in the domain translation speed in y-dir.
!
! OUTPUT:
!
! ubar Redefined base state x velocity component (m/s)
! vbar Redefined base state y velocity component (m/s)
!
! u Redefined total u-velocity (m/s)
! v Redefined total v-velocity (m/s)
! w Redefined total w-velocity (m/s)
! ptprt Redefined perturbation potential temperature at tpast (K)
! pprt Redefined perturbation pressure at tpast (Pascal)
! qv Redefined water vapor specific humidity at tpast (kg/kg)
! qc Redefined cloud water mixing ratio at tpast (kg/kg)
! qr Redefined rain water mixing ratio at tpast (kg/kg)
! qi Redefined cloud ice mixing ratio at tpast (kg/kg)
! qs Redefined snow mixing ratio at tpast (kg/kg)
! qh Redefined hail mixing ratio at tpast (kg/kg)
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
! tem2 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INCLUDE 'timelvls.inc'
INTEGER :: nx,ny,nz ! Number of grid points in x, y, z directions
REAL :: ubar (nx,ny,nz) ! Base state x velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state y velocity (m/s)
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 :: ptprt(nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz,nt) ! Perturbation pressure (Pascal)
REAL :: qv (nx,ny,nz,nt) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz,nt) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz,nt) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz,nt) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz,nt) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz,nt) ! Hail mixing ratio (kg/kg)
REAL :: qvbar(nx,ny,nz) ! Base-state water vapor mixing ratio (kg/kg)
REAL :: uchange ! Change in domain translation speed in x-dir.
REAL :: vchange ! Change in domain translation speed in y-dir.
REAL :: tem1 (nx,ny,nz) ! Temporary work array.
REAL :: tem2 (nx,ny,nz) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
INTEGER :: newebc ! East boundary condition parameter.
INTEGER :: newwbc ! West boundary condition parameter.
INTEGER :: newnbc ! North boundary condition parameter.
INTEGER :: newsbc ! South boundary condition parameter.
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global constants that control model
INCLUDE 'bndry.inc' ! Boundary condition parameters
INCLUDE 'mp.inc' ! Message passing parameters.
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Reset the domain translation speed if umove or vmove given in the
! input file is different from those in the restart data.
!
!-----------------------------------------------------------------------
!
WRITE(6,'((/1x,a,f10.3)/)') &
'ATTENTION: the domain translation speed was changed at t=', &
curtim
WRITE(6,'(2(/1x,2(a,f10.4))/)') &
'umove_new=',umove ,' vmove_new=',vmove , &
'umove_old=',umove-uchange,' vmove_old=',vmove-vchange
DO k=1,nz
DO j=1,ny
DO i=1,nx
ubar(i,j,k)=ubar(i,j,k)-uchange
vbar(i,j,k)=vbar(i,j,k)-vchange
u(i,j,k,tpast )=u(i,j,k,tpast )-uchange
u(i,j,k,tpresent)=u(i,j,k,tpresent)-uchange
u(i,j,k,tfuture )=u(i,j,k,tfuture )-uchange
v(i,j,k,tpast )=v(i,j,k,tpast )-vchange
v(i,j,k,tpresent)=v(i,j,k,tpresent)-vchange
v(i,j,k,tfuture )=v(i,j,k,tfuture )-vchange
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Adjust tpast values of all the time dependent variables
! such that their local time derivatives in the new moving reference
! frame satisfies the Galilean transformation:
!
! time derivative of a variable in the old system =
!
! time derivative of the variable in the new system
! - uchange *d(var)/dx - vchange*d(var)/dy
!
! where the time derivative is defined as
!
! ( var( tpresent ) - var( tpast ) )/dt.
!
!-----------------------------------------------------------------------
!
CALL galilei(nx,ny,nz, u, 2,nx-1,2,ny-2,1,nz-1, &
uchange,vchange,tem1)
CALL galilei(nx,ny,nz, v, 2,nx-2,2,ny-1,1,nz-1, &
uchange,vchange,tem1)
CALL galilei(nx,ny,nz, w, 2,nx-2,2,ny-2,1,nz , &
uchange,vchange,tem1)
CALL galilei(nx,ny,nz, ptprt,2,nx-2,2,ny-2,1,nz-1, &
uchange,vchange,tem1)
CALL galilei(nx,ny,nz, pprt, 2,nx-2,2,ny-2,1,nz-1, &
uchange,vchange,tem1)
CALL galilei(nx,ny,nz, qv, 2,nx-2,2,ny-2,1,nz-1, &
uchange,vchange,tem1)
CALL galilei(nx,ny,nz, qc, 2,nx-2,2,ny-2,1,nz-1, &
uchange,vchange,tem1)
CALL galilei(nx,ny,nz, qr, 2,nx-2,2,ny-2,1,nz-1, &
uchange,vchange,tem1)
CALL galilei(nx,ny,nz, qi, 2,nx-2,2,ny-2,1,nz-1, &
uchange,vchange,tem1)
CALL galilei(nx,ny,nz, qs, 2,nx-2,2,ny-2,1,nz-1, &
uchange,vchange,tem1)
CALL galilei(nx,ny,nz, qh, 2,nx-2,2,ny-2,1,nz-1, &
uchange,vchange,tem1)
!
!****************************************************************
! Remove possible negative values for Q-fields resulted from
! the calling of 'galilei' subroutine
!****************************************************************
!
DO i=2,nx-2
DO j=2,ny-2
DO k=1,nz-1
qv(i,j,k,1) = AMAX1(qv(i,j,k,1), 0.0)
qc(i,j,k,1) = AMAX1(qc(i,j,k,1), 0.0)
qr(i,j,k,1) = AMAX1(qr(i,j,k,1), 0.0)
qi(i,j,k,1) = AMAX1(qi(i,j,k,1), 0.0)
qs(i,j,k,1) = AMAX1(qs(i,j,k,1), 0.0)
qh(i,j,k,1) = AMAX1(qh(i,j,k,1), 0.0)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Apply boundary conditions to redefined dynamical, thermodynamical
! and microphysical variables at tpast. Note: there is no
! provision for radiation lateral boundary conditions; if radiation
! conditions are chosen, zero-gradient conditions are applied instead.
!
!-----------------------------------------------------------------------
!
IF ((ebc /= 1).AND.(ebc /= 2).AND.(ebc /= 3).AND.(ebc /= 0)) THEN
newebc = 3
ELSE
newebc = ebc
END IF
IF ((wbc /= 1).AND.(wbc /= 2).AND.(wbc /= 3).AND.(wbc /= 0)) THEN
newwbc = 3
ELSE
newwbc = wbc
END IF
IF ((nbc /= 1).AND.(nbc /= 2).AND.(nbc /= 3).AND.(nbc /= 0)) THEN
newnbc = 3
ELSE
newnbc = nbc
END IF
IF ((sbc /= 1).AND.(sbc /= 2).AND.(sbc /= 3).AND.(sbc /= 0)) THEN
newsbc = 3
ELSE
newsbc = sbc
END IF
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(u,nx,ny,nz,newebc,newwbc,1,mptag,tem2)
CALL mprecv2dew(u,nx,ny,nz,newebc,newwbc,1,mptag,tem2)
CALL mpsend2dns(u,nx,ny,nz,newnbc,newsbc,1,mptag,tem2)
CALL mprecv2dns(u,nx,ny,nz,newnbc,newsbc,1,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL bcu(nx,ny,nz,0., u(1,1,1,tfuture), &
0.,0.,0.,0.,newebc,newwbc,newnbc,newsbc,tbc,bbc)
CALL acct_stop_inter
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(v,nx,ny,nz,newebc,newwbc,2,mptag,tem2)
CALL mprecv2dew(v,nx,ny,nz,newebc,newwbc,2,mptag,tem2)
CALL mpsend2dns(v,nx,ny,nz,newnbc,newsbc,2,mptag,tem2)
CALL mprecv2dns(v,nx,ny,nz,newnbc,newsbc,2,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL bcv(nx,ny,nz,0., v(1,1,1,tfuture), &
0.,0.,0.,0.,newebc,newwbc,newnbc,newsbc,tbc,bbc)
CALL acct_stop_inter
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(tem1,nx,ny,nz,newebc,newwbc,3,mptag,tem2)
CALL mprecv2dew(tem1,nx,ny,nz,newebc,newwbc,3,mptag,tem2)
CALL mpsend2dns(tem1,nx,ny,nz,newnbc,newsbc,3,mptag,tem2)
CALL mprecv2dns(tem1,nx,ny,nz,newnbc,newsbc,3,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL lbcw(nx,ny,nz,0.0,w,tem1, &
0.,0.,0.,0.,newebc,newwbc,newnbc,newsbc)
CALL acct_stop_inter
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(ptprt,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mprecv2dew(ptprt,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mpsend2dns(ptprt,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
CALL mprecv2dns(ptprt,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL bcpt(nx,ny,nz,0.,ptprt, &
0.,0.,0.,0.,newebc,newwbc,newnbc,newsbc,tbc,bbc)
CALL acct_stop_inter
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,1)=pprt(i,j,2,1)
END DO
END DO
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(pprt,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mprecv2dew(pprt,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mpsend2dns(pprt,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
CALL mprecv2dns(pprt,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL bcp(nx,ny,nz,0., &
pprt(1,1,1,tfuture),0.,0.,0.,0.,tem1(1,1,1), &
newebc,newwbc,newnbc,newsbc,tbc,bbc)
CALL acct_stop_inter
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(qv,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mprecv2dew(qv,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mpsend2dns(qv,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
CALL mprecv2dns(qv,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL bcqv(nx,ny,nz,0., qv,qvbar, &
0.,0.,0.,0.,newebc,newwbc,newnbc,newsbc,tbc,bbc)
CALL acct_stop_inter
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(qc,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mprecv2dew(qc,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mpsend2dns(qc,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
CALL mprecv2dns(qc,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL bcq(nx,ny,nz,0., qc, &
0.,0.,0.,0.,newebc,newwbc,newnbc,newsbc,tbc,bbc)
CALL acct_stop_inter
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(qr,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mprecv2dew(qr,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mpsend2dns(qr,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
CALL mprecv2dns(qr,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL bcq(nx,ny,nz,0., qr, &
0.,0.,0.,0.,newebc,newwbc,newnbc,newsbc,tbc,bbc)
CALL acct_stop_inter
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(qi,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mprecv2dew(qi,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mpsend2dns(qi,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
CALL mprecv2dns(qi,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL bcq(nx,ny,nz,0., qi, &
0.,0.,0.,0.,newebc,newwbc,newnbc,newsbc,tbc,bbc)
CALL acct_stop_inter
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(qs,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mprecv2dew(qs,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mpsend2dns(qs,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
CALL mprecv2dns(qs,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL bcq(nx,ny,nz,0., qs, &
0.,0.,0.,0.,newebc,newwbc,newnbc,newsbc,tbc,bbc)
CALL acct_stop_inter
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(qh,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mprecv2dew(qh,nx,ny,nz,newebc,newwbc,0,mptag,tem2)
CALL mpsend2dns(qh,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
CALL mprecv2dns(qh,nx,ny,nz,newnbc,newsbc,0,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL bcq(nx,ny,nz,0., qh, &
0.,0.,0.,0.,newebc,newwbc,newnbc,newsbc,tbc,bbc)
CALL acct_stop_inter
RETURN
END SUBROUTINE adjuvmv
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GALILEI ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE galilei(nx,ny,nz, var, ibgn,iend,jbgn,jend,kbgn,kend, & 11
uchange,vchange, tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Adjust tpast values of all the time dependent variables
! such that their local time derivatives in the new moving reference
! frame satisfies the Galilean transformation:
!
! time derivative of a variable in the old system =
!
! time derivative of the variable in the new system
! - vchange*d(var)/dx - vchange*d(var)/dy
!
! where the time derivative is defined as
!
! ( var( tpresent ) - var( tpast ) )/dt.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
! AUTHORS: Alan Shapiro and Steve Lazarus
! 5/23/93.
!
! MODIFICATION HISTORY:
! 9/7/93 (MX)
! Code cleanup with modifications.
!
! 3/25/94 (MX)
! A bug fixed to remove the dependency of var(tpart,i) on
! var(tpast,n-i) in loop 100.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! 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 z-direction (vertical)
!
! var Any one of the four-dimensional variable that is advected
! using leapfrog-centered scheme.
! uchange New-old domain translation speed in x direction
! vchange New-old domain translation speed in y direction
! ibgn i-index where multiplication begins.
! iend i-index where multiplication ends.
! jbgn j-index where multiplication begins.
! jend j-index where multiplication ends.
! kbgn k-index where multiplication begins.
! kend k-index where multiplication ends.
!
! OUTPUT:
!
! var var at tpast.
!
! WORK ARRAY:
!
! tem1 Work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INCLUDE 'timelvls.inc'
INTEGER :: nx,ny,nz ! Number of grid points in x, y, z directions
REAL :: var(nx,ny,nz,nt) ! Any one of the advect variables.
INTEGER :: ibgn,iend,jbgn,jend,kbgn,kend
REAL :: uchange ! New-old domain translation speed in x-dir
REAL :: vchange ! New-old domain translation speed in y-dir
REAL :: tem1(nx,ny,nz) ! Local work array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
REAL :: dvardx, dvardy
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global constants that control model
INCLUDE 'grid.inc' ! Grid parameters
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO k=kbgn,kend
DO j=jbgn,jend
DO i=ibgn,iend
dvardx = 0.25*dxinv* &
(var(i+1,j,k,tpast) + var(i+1,j,k,tpresent) &
- var(i-1,j,k,tpast) - var(i-1,j,k,tpresent))
dvardy = 0.25*dyinv* &
(var(i,j+1,k,tpast) + var(i,j+1,k,tpresent) &
- var(i,j-1,k,tpast) - var(i,j-1,k,tpresent))
tem1(i,j,k) = var(i,j,k,tpast) &
- dtbig*(uchange*dvardx+vchange*dvardy)
END DO
END DO
END DO
DO k=kbgn,kend
DO j=jbgn,jend
DO i=ibgn,iend
var(i,j,k,tpast) = tem1(i,j,k)
END DO
END DO
END DO
RETURN
END SUBROUTINE galilei
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE AUTOTRANS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE autotrans(nx,ny,nz, & 1
ubar,vbar,u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
uchange, vchange, tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! This subroutine computes the optimum (least squares) pattern translation
! speed based on vertical velocity data. This technique is described in:
!
! T. Gal-Chen, 1982: "Errors in Fixed and Moving Frame of References:
! Applications for Conventional and Doppler Radar Analysis",
! J. Atmos. Sci., Vol. 39, 1982. 2279-2300.
!
!
! Special note: Many of the arrays in the AUTOTRANS subroutine call
! are not actually used in AUTOTRANS. However, these arrays will
! likely be used in future tests of this subroutine.
!
!-----------------------------------------------------------------------
!
! AUTHORS: Alan Shapiro and Steve Lazarus
! 5/02/93.
!
! MODIFICATION HISTORY:
!
! 4/18/94 Alan Shapiro and Robb Randall
! Changed weighting so that stronger storm contributes more to
! the calculation of bigu and bigv.
! Also, only regions of updraft (w>0) are used in the calculation.
!
!
!-----------------------------------------------------------------------
!
! 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 z-direction (vertical)
!
! ubar Model's base state x velocity component (m/s)
! vbar Model's base state y velocity component (m/s)
!
! u Total u-velocity (m/s)
! v Total v-velocity (m/s)
! w Total w-velocity (m/s)
! ptprt Perturbation potential temperature (K)
! pprt Perturbation pressure (Pascal)
! qv Water vapor specific humidity (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rain water mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
!
!
! OUTPUT:
!
! uchange Change in umove (=0 if it's not time to change).
! vchange Change in vmove (=0 if it's not time to change).
!
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INCLUDE 'timelvls.inc'
INTEGER :: nx,ny,nz ! Number of grid points in x, y, z directions
REAL :: ubar (nx,ny,nz) ! Base state x velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state y velocity (m/s)
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 :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz,nt) ! Perturbation pressure (Pascal)
REAL :: qv (nx,ny,nz,nt) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz,nt) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz,nt) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz,nt) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz,nt) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz,nt) ! Hail mixing ratio (kg/kg)
REAL :: tem1 (nx,ny,nz) ! Temporary work array.
!
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc' ! Global constants that control model
INCLUDE 'grid.inc' ! Grid parameters
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL :: w1111, w1110, w2111, &
w0111, w1211, w1011 ! Vertical velocity.
INTEGER :: i, j, k
INTEGER :: kmin ! Minimum k-level for calculation of optimal
! grid-translation.
INTEGER :: kmax ! Maximum k-level for calculation of optimal
! grid-translation.
REAL :: dwdt, dwdx, dwdy ! Spatial and temporal derivatives of w.
REAL :: bigu ! Optimum perturbation translation speed in
! x direction (computed each time step).
REAL :: bigv ! Optimum perturbation translation speed in
! y direction (computed each time step).
REAL :: speed ! speed = sqrt(bigu**2 + bigv**2)
REAL :: a, b, c, d, e ! Coefficients in solution for bigu and bigv.
REAL :: dtinv ! Reciprocal of dtbig
REAL :: denom ! Denominator of formula for bigu, bigv
INTEGER :: ntlev ! Number of big time steps in every time
! window (rounded down)
REAL :: velmax ! Maximum allowable perturbation translation
! speed.
REAL :: uchange ! Optimum perturbation translation speed in
! x direction (averaged over a time window).
REAL :: vchange ! Optimum perturbation translation speed in
! y direction (averaged over a time window).
REAL :: sumu ! Running sum of bigu
REAL :: sumv ! Running sum of bigv
SAVE sumu, sumv
INTEGER :: ichnge ! Parameter that signifies whether grid
! translation will be changed this time step.
SAVE ichnge
DATA sumu /0./
DATA sumv /0./
DATA ichnge /0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Initialize uchange and vchange to zero.
!
!-----------------------------------------------------------------------
!
uchange = 0.
vchange = 0.
!
!-----------------------------------------------------------------------
!
! Define maximum value of grid translation in terms of CFL criterion.
!
!-----------------------------------------------------------------------
!
dtinv = 1./dtbig
velmax = dtinv*SQRT(dx*dy)
!-----------------------------------------------------------------------
!
! Compute optimal grid translation components based on low-level
! vertical velocity data.
!
!-----------------------------------------------------------------------
!
a = 0.
b = 0.
c = 0.
d = 0.
e = 0.
kmin = 3
kmax = 2 + chkdpth/dz
DO k=kmin, kmax
DO j=2,ny-2
DO i=2,nx-2
w1111 = MAX(0.0, w(i,j,k,tpresent))
w1110 = MAX(0.0, w(i,j,k,tpast))
w2111 = MAX(0.0, w(i+1,j,k,tpresent))
w0111 = MAX(0.0, w(i-1,j,k,tpresent))
w1211 = MAX(0.0, w(i,j+1,k,tpresent))
w1011 = MAX(0.0, w(i,j-1,k,tpresent))
dwdt = dtinv*(w1111**4 - w1110**4)
dwdx = 0.5*dxinv*(w2111**4 - w0111**4)
dwdy = 0.5*dyinv*(w1211**4 - w1011**4)
a = a + dwdt*dwdx
b = b + dwdx*dwdx
c = c + dwdx*dwdy
d = d + dwdt*dwdy
e = e + dwdy*dwdy
END DO
END DO
END DO
denom = c*c - b*e
IF (denom == 0.) THEN
PRINT *, 'Warning: denom = 0. bigu, bigv set to 0'
bigu = 0.
bigv = 0.
ELSE
bigu = (a*e - c*d)/denom
bigv = (b*d - c*a)/denom
PRINT *, 'bigu, bigv = ', bigu, bigv
END IF
!
!-----------------------------------------------------------------------
!
! Check that bigu and bigv are not too big. If the speed of
! translation exceeds velmax (based on CFL criterion) then set
! bigu and bigv to zero.
!
!-----------------------------------------------------------------------
!
speed = SQRT(bigu*bigu + bigv*bigv)
IF (speed > velmax) THEN
bigu = 0.
bigv = 0.
PRINT *, 'Warning! speed= ', speed, ' exceeds velmax= ', velmax
PRINT *, 'Set bigu and bigv to zero.'
PRINT *, 'bigu, bigv = ', bigu, bigv
END IF
!
!-----------------------------------------------------------------------
!
! Sum bigu and bigv each time step.
!
!-----------------------------------------------------------------------
!
sumu = sumu + bigu
sumv = sumv + bigv
!
!-----------------------------------------------------------------------
!
! Update umove and vmove at the end of every time window.
!
!-----------------------------------------------------------------------
!
ntlev = IFIX(twindow/dtbig)
ichnge = MOD(nstep, ntlev)
IF (ichnge == 0) THEN
uchange = sumu/ntlev
vchange = sumv/ntlev
PRINT *, 'umove and vmove are being updated at time =', curtim
PRINT *, 'Old umove and vmove = ', umove, vmove
PRINT *, 'ntlev, uchange, vchange = ', ntlev, uchange, vchange
umove = umove + uchange
vmove = vmove + vchange
PRINT *, 'New umove and vmove will be = ', umove, vmove
!
!-----------------------------------------------------------------------
!
! Set sumu and sumv to zero for the beginning of next time window.
!
!-----------------------------------------------------------------------
!
sumu = 0.
sumv = 0.
END IF
RETURN
END SUBROUTINE autotrans