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