!######################################################################## !######################################################################## !###### ###### !###### SUBROUTINE SOLVQ ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !######################################################################## !######################################################################## SUBROUTINE solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag, & q,u,v,wcont, qdteb,qdtwb,qdtnb,qdtsb,rhostr,kmh,kmv,rprntl, & x,y,z,zp, mapfct, j1,j2,j3,j3inv, & qccumsrc,qrcumsrc,qicumsrc,qscumsrc,exbcbuf,bcrlx, & qadv,qmix,tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, & tem9,tem10,tem11, tem1_0,tem2_0,tem3_0) !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Coordinate the time integration of the equations for the water ! substance quantities qc,qr,qi,qs and qh, i.e. all water variables ! except the water vapor specific humidity qv. ! !----------------------------------------------------------------------- ! ! AUTHOR: Ming Xue ! 10/10/91 ! ! MODIFICATION HISTORY: ! ! 5/20/92 (K. Droegemeier and M. Xue) ! Added full documentation. ! ! 2/10/93 (K. Droegemeier) ! Cleaned up documentation. ! ! 8/12/95 (M. Xue) ! Added flag qflag. ! ! 1/25/96 (Donghai Wang & Yuhe Liu) ! Added the map projection factor to ARPS governing equations. ! ! 4/1/96 (Donghai Wang, X. Song and M. Xue) ! Added the implicit treatment for the vertical mixing. ! ! 3/24/1997 (M. Xue) ! Code to handle the case of forward time integration added. ! ! 7/10/1997 (Fanyou Kong -- CMRP) ! Added the positive definite advection scheme option (sadvopt = 5) ! ! 4/15/1998 (Donghai Wang) ! Added the source terms to the right hand terms of the qc,qr,qi,qs ! equations due to the K-F cumulus parameterization. ! ! 7/17/1998 (M. Xue) ! Changed call to ADVQFCT. ! !----------------------------------------------------------------------- ! Use modules !----------------------------------------------------------------------- USE bndry USE exbc USE globcst !----------------------------------------------------------------------- ! Variable Declarations: !----------------------------------------------------------------------- IMPLICIT NONE ! Force explicit declarations INTEGER, PARAMETER :: nt=3 ! The no. of t-levels of t-dependent arrays. INTEGER, PARAMETER :: tpast =1 ! Index of time level for the past time. INTEGER, PARAMETER :: tpresent=2 ! Index of time level for the present time. INTEGER, PARAMETER :: tfuture =3 ! Index of time level for the future time. INTEGER, INTENT(IN) :: nx, ny, nz ! Number of grid points in 3 directions REAL, INTENT(IN) :: dtbig1 ! Local big time step INTEGER, INTENT(IN) :: qflag ! Indicator for water/ice species REAL, INTENT(INOUT) :: q (nx,ny,nz,nt) ! One of the water/ice variables (kg/kg) REAL, INTENT(IN ) :: u (nx,ny,nz,nt) ! Total u-velocity (m/s) REAL, INTENT(IN ) :: v (nx,ny,nz,nt) ! Total v-velocity (m/s) REAL, INTENT(IN ) :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s) REAL, INTENT(IN ) :: qdteb(ny,nz) ! Time tendency of liquid/ice at e-boundary REAL, INTENT(IN ) :: qdtwb(ny,nz) ! Time tendency of liquid/ice at w-boundary REAL, INTENT(IN ) :: qdtnb(nx,nz) ! Time tendency of liquid/ice at n-boundary REAL, INTENT(IN ) :: qdtsb(nx,nz) ! Time tendency of liquid/ice at s-boundary REAL, INTENT(IN ) :: rhostr(nx,ny,nz) ! Base state density rhobar times j3. REAL, INTENT(IN ) :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for momentum. ( m**2/s ) REAL, INTENT(IN ) :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for momentum. ( m**2/s ) REAL, INTENT(IN ) :: rprntl(nx,ny,nz) ! Reciprocal of Prandtl number REAL, INTENT(IN ) :: x (nx) ! x-coord. of the physical and computational grid. ! Defined at u-point. REAL, INTENT(IN ) :: y (ny) ! y-coord. of the physical and computational grid. ! Defined at v-point. REAL, INTENT(IN ) :: z (nz) ! z-coord. of the computational grid. ! Defined at w-point on the staggered grid. REAL, INTENT(IN ) :: zp (nx,ny,nz) ! Physical height coordinate defined at w-point of ! the staggered grid. REAL, INTENT(IN ) :: mapfct(nx,ny) ! Map factors at scalar points REAL, INTENT(IN ) :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian defined as - d( zp )/d( x ). REAL, INTENT(IN ) :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian defined as - d( zp )/d( y ). REAL, INTENT(IN ) :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian defined as d( zp )/d( z ). REAL, INTENT(IN ) :: j3inv (nx,ny,nz) ! Coordinate transformation Jacobian defined as d( zp )/d( z ). REAL, INTENT(IN ) :: qccumsrc(nx,ny,nz) ! Source in qc-equation due to cumulus parameterization REAL, INTENT(IN ) :: qrcumsrc(nx,ny,nz) ! Source in qr-equation due to cumulus parameterization REAL, INTENT(IN ) :: qicumsrc(nx,ny,nz) ! Source in qi-equation due to cumulus parameterization REAL, INTENT(IN ) :: qscumsrc(nx,ny,nz) ! Source in qs-equation due to cumulus parameterization INTEGER, INTENT(IN) :: exbcbufsz ! EXBC buffer size REAL, INTENT(IN ) :: exbcbuf( exbcbufsz ) ! EXBC buffer array REAL, INTENT(IN ) :: bcrlx (nx,ny) ! EXBC relaxation coefficients !----------------------------------------------------------------------- ! Work arrays: !----------------------------------------------------------------------- REAL, INTENT(INOUT) :: qadv(nx,ny,nz) ! Advection of water/ice substance (kg/(m**3*s)). REAL, INTENT(INOUT) :: qmix(nx,ny,nz) ! Total mixing of water/ice substance (kg/(m**3*s)). REAL, INTENT(INOUT) :: tem1 (nx,ny,nz) ! Temporary work array REAL, INTENT(INOUT) :: tem2 (nx,ny,nz) ! Temporary work array REAL, INTENT(INOUT) :: tem3 (nx,ny,nz) ! Temporary work array REAL, INTENT(INOUT) :: tem4 (nx,ny,nz) ! Temporary work array REAL, INTENT(INOUT) :: tem5 (nx,ny,nz) ! Temporary work array REAL, INTENT(INOUT) :: tem6 (nx,ny,nz) ! Temporary work array REAL, INTENT(INOUT) :: tem7 (nx,ny,nz) ! Temporary work array REAL, INTENT(INOUT) :: tem8 (nx,ny,nz) ! Temporary work array REAL, INTENT(INOUT) :: tem9 (nx,ny,nz) ! Temporary work array REAL, INTENT(INOUT) :: tem10 (nx,ny,nz) ! Temporary work array REAL, INTENT(INOUT) :: tem11 (nx,ny,nz) ! Temporary work array REAL, INTENT(INOUT) :: tem1_0(0:nx,0:ny,0:nz) ! Non-standard size work array REAL, INTENT(INOUT) :: tem2_0(0:nx,0:ny,0:nz) ! Non-standard size work array REAL, INTENT(INOUT) :: tem3_0(0:nx,0:ny,0:nz) ! Non-standard size work array !----------------------------------------------------------------------- ! Misc. local variables: !----------------------------------------------------------------------- INTEGER :: i,j,k INTEGER :: tstrtlvl INTEGER :: ebc1,wbc1,nbc1,sbc1 INTEGER :: nq0exb,nqdtexb REAL :: deltat REAL :: cpu0, f_cputime ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! !----------------------------------------------------------------------- ! Compute the advection term for a general water substance ! q and store the result in qadv. !----------------------------------------------------------------------- cpu0 = f_cputime() IF(sadvopt==4) THEN ! Forward-based FCT scheme deltat = dtbig1 tstrtlvl = tpresent ELSE deltat = dtbig1*2 tstrtlvl = tpast END IF IF (sadvopt==1 .OR. sadvopt==2 .OR. sadvopt==3 ) THEN ! 2nd or 4th order centerd sc DO i=1,nx DO j=1,ny DO k=1,nz tem7(i,j,k) = 0.0 END DO END DO END DO CALL advq(nx,ny,nz,q,u,v,wcont,rhostr, tem7, mapfct, & qadv, & tem1,tem2,tem3,tem4,tem5,tem6) ELSE IF( sadvopt==4 .OR. sadvopt==5 ) THEN ! FCT advection CALL advqfct(nx,ny,nz,dtbig1,q,u,v,wcont, rhostr,mapfct,j3,qadv, & tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, & tem9,tem10,tem11,tem1_0,tem2_0,tem3_0) END IF advs_cpu_time = advs_cpu_time + f_cputime() - cpu0 !----------------------------------------------------------------------- ! Compute the mixing terms for the general water substance (q) ! equation, including both physical and computational mixing. ! Store the result in array qmix. !----------------------------------------------------------------------- CALL mixq(nx,ny,nz,q(1,1,1,tstrtlvl),rhostr,kmh,kmv,rprntl, & x,y,z,zp,mapfct, j1,j2,j3,j3inv, & qmix, & tem1,tem2,tem3,tem4,tem5,tem6) !----------------------------------------------------------------------- ! Call BRLXQ to added to qmix the additional boundary relaxation and ! spatial mixing in the boundary zone ! ! qmix = qmix + qexbc_term !----------------------------------------------------------------------- IF ( lbcopt==2 .AND. mgrid==1 ) THEN cpu0 = f_cputime() CALL brlxq(nx,ny,nz, deltat*0.5, qflag, q,rhostr, qmix, & exbcbuf(nqv0exb), exbcbuf(nqc0exb), & exbcbuf(nqr0exb), exbcbuf(nqi0exb), & exbcbuf(nqs0exb), exbcbuf(nqh0exb), & exbcbuf(nqvdtexb),exbcbuf(nqcdtexb), & exbcbuf(nqrdtexb),exbcbuf(nqidtexb), & exbcbuf(nqsdtexb),exbcbuf(nqhdtexb),bcrlx, & tem1,tem2,tem3,tem4) bc_cpu_time = bc_cpu_time + f_cputime() - cpu0 END IF !----------------------------------------------------------------------- ! Calculate qforce, store the result in tem1. !----------------------------------------------------------------------- DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 tem1(i,j,k)=-qadv(i,j,k)+qmix(i,j,k) IF(qflag == 2) THEN tem1(i,j,k)= tem1(i,j,k)+rhostr(i,j,k)*qccumsrc(i,j,k) ELSE IF(qflag == 3) THEN tem1(i,j,k)= tem1(i,j,k)+rhostr(i,j,k)*qrcumsrc(i,j,k) ELSE IF(qflag == 4) THEN tem1(i,j,k)= tem1(i,j,k)+rhostr(i,j,k)*qicumsrc(i,j,k) ELSE IF(qflag == 5) THEN tem1(i,j,k)= tem1(i,j,k)+rhostr(i,j,k)*qscumsrc(i,j,k) END IF END DO END DO END DO !----------------------------------------------------------------------- ! Treat the vertically implicit mixing term. !----------------------------------------------------------------------- IF (trbvimp==1) THEN ! Vertical implicit application cpu0 = f_cputime() DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 tem2(i,j,k)=rhostr(i,j,k)*kmv(i,j,k)*rprntl(i,j,k) & *j3inv(i,j,k)*j3inv(i,j,k) END DO END DO END DO CALL vmiximps(nx,ny,nz,deltat*0.5,q(1,1,1,tstrtlvl),rhostr, & tem2,tem1,tem3,tem4,tem5,tem6) tmix_cpu_time = tmix_cpu_time + f_cputime() - cpu0 END IF !----------------------------------------------------------------------- ! Integrate forward by one timestep the general water substance ! (q) equation, yielding q at time = tfuture. !----------------------------------------------------------------------- IF( nestgrd==1 .AND. mgrid/=1 ) THEN cpu0 = f_cputime() DO k=2,nz-2 DO j=2,ny-2 DO i=2,nx-2 q(i,j,k,tfuture)=q(i,j,k,tstrtlvl)+deltat*tem1(i,j,k) & /rhostr(i,j,k) END DO END DO END DO misc_cpu_time = misc_cpu_time + f_cputime() - cpu0 cpu0 = f_cputime() !----------------------------------------------------------------------- ! Set the boundary conditions on q for an adaptive (nested) ! grid run. If using only one grid, this portion of the code is ! skipped....proceed to next comment block. !----------------------------------------------------------------------- DO k=2,nz-2 DO i=1,nx-1 q(i, 1,k,tfuture)=2*q(i, 1,k,tpresent)-q(i, 1,k,tpast) q(i,ny-1,k,tfuture)=2*q(i,ny-1,k,tpresent)-q(i,ny-1,k,tpast) END DO END DO DO k=2,nz-2 DO j=2,ny-2 q( 1,j,k,tfuture)=2*q( 1,j,k,tpresent)-q( 1,j,k,tpast) q(nx-1,j,k,tfuture)=2*q(nx-1,j,k,tpresent)-q(nx-1,j,k,tpast) END DO END DO CALL bcq(nx,ny,nz,dtbig1,q,qdteb,qdtwb,qdtnb,qdtsb,0,0,0,0,tbc,bbc) bc_cpu_time = bc_cpu_time + f_cputime() - cpu0 ELSE cpu0 = f_cputime() DO k=2,nz-2 DO j=1,ny-1 DO i=1,nx-1 q(i,j,k,tfuture)=q(i,j,k,tstrtlvl)+deltat*tem1(i,j,k) & /rhostr(i,j,k) END DO END DO END DO !----------------------------------------------------------------------- ! Set the boundary conditions on q for a NON-adaptive (uniform) ! grid run. !----------------------------------------------------------------------- misc_cpu_time = misc_cpu_time + f_cputime() - cpu0 cpu0 = f_cputime() IF ( lbcopt==1 ) THEN ebc1=ebc wbc1=wbc sbc1=sbc nbc1=nbc IF( ebc==4 ) ebc1=0 IF( wbc==4 ) wbc1=0 IF( nbc==4 ) nbc1=0 IF( sbc==4 ) sbc1=0 CALL bcq(nx,ny,nz, dtbig1, q, qdteb,qdtwb,qdtnb,qdtsb, & ebc1,wbc1,nbc1,sbc1,tbc,bbc) ELSE CALL bcq(nx,ny,nz, dtbig1, q, qdteb,qdtwb,qdtnb,qdtsb, & 3,3,3,3,tbc,bbc) ! Zero-gradient condition will be ! reset if exbc data is available for q IF ( qflag==1 ) THEN nq0exb = nqv0exb nqdtexb = nqvdtexb ELSE IF ( qflag==2 ) THEN nq0exb = nqc0exb nqdtexb = nqcdtexb ELSE IF ( qflag==3 ) THEN nq0exb = nqr0exb nqdtexb = nqrdtexb ELSE IF ( qflag==4 ) THEN nq0exb = nqi0exb nqdtexb = nqidtexb ELSE IF ( qflag==5 ) THEN nq0exb = nqs0exb nqdtexb = nqsdtexb ELSE IF ( qflag==6 ) THEN nq0exb = nqh0exb nqdtexb = nqhdtexb END IF CALL exbcq(nx,ny,nz, qflag, curtim+dtbig,q(1,1,1,tfuture), & exbcbuf(nq0exb),exbcbuf(nqdtexb)) END IF bc_cpu_time = bc_cpu_time + f_cputime() - cpu0 END IF RETURN END SUBROUTINE solvq