! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE SOLVUV ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! SUBROUTINE SOLVUV(nx,ny,nz, exbcbufsz, dtsml1, curtsml, & & u,v,wcont,pprt, & & udteb,udtwb,udtnb,udtsb, & & vdteb,vdtwb,vdtnb,vdtsb, & & rhostr, uforce,vforce,ubar,vbar, & & x,y,z,zp, mapfct, j1,j2,j3,j3inv, & & exbcbuf,rhofct, & & rhostru,rhostrv,rhostrw,wpgrad, & & upgrad,vpgrad,tem1,tem2,tem3) ! !####################################################################### ! ! PURPOSE: ! ! Coordinate the time stepping of u and v momentum equations and the ! setting of boundary conditions. ! !####################################################################### ! ! ! AUTHOR: Ming Xue ! 10/10/92 ! ! MODIFICATION HISTORY: ! ! 5/20/92 (K. Droegemeier and M. Xue) ! Added full documentation. ! ! 6/07/92 ( M. Xue) ! Now only tfuture fields of u, v, w and pprt are passed in. ! ! 2/10/93 (K. Droegemeier) ! Cleaned up documentation. ! ! 4/10/93 (M. Xue & Hao Jin) ! Add the terrain. ! ! 9/1/94 (Y. Lu) ! Cleaned up documentation. ! ! 1/25/96 (Donghai Wang & Yuhe Liu) ! Added the map projection factor to ARPS governing equations. ! ! 10/16/97 (Donghai Wang) ! Using total density (rho) for the calculation of the pressure ! gradient force terms. ! ! 11/06/97 (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. ! !####################################################################### ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! dtsml1 Local value of small time step size ! curtsml Local curttim at a small time step ! ! u x component of velocity at tfuture (m/s) ! v y component of velocity at tfuture (m/s) ! wcont Contravariant vertical velocity (m/s) ! ! pprt Perturbation pressure at tfuture (Pascal) ! ! udteb Time tendency of the u field at the east boundary ! udtwb Time tendency of the u field at the west boundary ! udtnb Time tendency of the u field at the north boundary ! udtsb Time tendency of the u field at the south boundary ! ! vdteb Time tendency of the v field at the east boundary ! vdtwb Time tendency of the v field at the west boundary ! vdtnb Time tendency of the v field at the north boundary ! vdtsb Time tendency of the v field at the south boundary ! ! rhostr Base state density rhobar times j3 (kg/m**3) ! ! uforce Acoustically inactive forcing terms in u-eq. ! (kg/(m*s)**2) ! vforce Acoustically inactive forcing terms in v-eq. ! (kg/(m*s)**2) ! ! ubar Base state x velocity component (m/s) ! vbar Base state y velocity component (m/s) ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp Vertical coordinate of grid points in physical space(m) ! ! mapfct Map factors at scalar, u and v points ! ! j1 Coordinate transformation Jacobian -d(zp)/dx ! j2 Coordinate transformation Jacobian -d(zp)/dy ! j3 Coordinate transformation Jacobian d(zp)/dz ! j3inv Inverse of j3 ! ! rhostru Average rhostr at u points (kg/m**3). ! rhostrv Average rhostr at v points (kg/m**3). ! rhostrw Average rhostr at w points (kg/m**3). ! ! OUTPUT: ! ! u x component of velocity at time tfuture updated ! in time by dtsml1 (m/s) ! v y component of velocity at time tfuture updated ! in time by dtsml1 (m/s) ! w Vertical component of Cartesian velocity at tfuture ! updated in time by dtsml1 (m/s) ! wpgrad Pressure gradient term in w-eq. (kg/(m*s)**2) ! ! WORK ARRAYS: ! ! upgrad Pressure gradient term in u-eq. (kg/(m*s)**2) ! vpgrad Pressure gradient term in v-eq. (kg/(m*s)**2) ! tem1 Temporary work array ! tem2 Temporary work array ! tem3 Temporary work array ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations: ! !####################################################################### ! ! Force explicit declarations implicit none ! Number of grid points in 3 directions integer nx, ny, nz ! u-velocity at tfuture (m/s) real u (nx,ny,nz) ! v-velocity at tfuture (m/s) real v (nx,ny,nz) ! Contravariant vertical velocity (m/s) real wcont (nx,ny,nz) ! Perturbation pressure at tfuture real pprt (nx,ny,nz) ! (Pascal) ! Time tendency of u at east boundary real udteb (ny,nz) ! Time tendency of u at west boundary real udtwb (ny,nz) ! Time tendency of u at north boundary real udtnb (nx,nz) ! Time tendency of u at south boundary real udtsb (nx,nz) ! Time tendency of v at east boundary real vdteb (ny,nz) ! Time tendency of v at west boundary real vdtwb (ny,nz) ! Time tendency of v at north boundary real vdtnb (nx,nz) ! Time tendency of v at south boundary real vdtsb (nx,nz) ! Base state density rhobar times j3. real rhostr(nx,ny,nz) ! Acoustically inactive forcing terms real uforce(nx,ny,nz) ! in u-momentum equation (kg/(m*s)**2) ! Acoustically inactive forcing terms real vforce(nx,ny,nz) ! in v-momentum equation (kg/(m*s)**2) real ubar(nx,ny,nz) real vbar(nx,ny,nz) ! x-coord. of the physical and compu- real x (nx) ! tational grid. Defined at u-point. ! y-coord. of the physical and compui- real y (ny) ! tational grid. Defined at v-point. ! z-coord. of the computational grid. real z (nz) ! Defined at w-point on the staggered ! grid. ! Physical height coordinate defined at real zp (nx,ny,nz) ! w-point of the staggered grid. ! Map factors at scalar, u and v points real mapfct(nx,ny,6) ! Coordinate transformation Jacobian real j1 (nx,ny,nz) ! defined as - d( zp )/d( x ). ! Coordinate transformation Jacobian real j2 (nx,ny,nz) ! defined as - d( zp )/d( y ). ! Coordinate transformation Jacobian real j3 (nx,ny,nz) ! defined as d( zp )/d( z ). ! Inverse of j3 real j3inv (nx,ny,nz) ! EXBC buffer size integer exbcbufsz ! EXBC buffer array real exbcbuf( exbcbufsz ) ! rho-factor: rhobar/rho real rhofct(nx,ny,nz) ! Averaged rhostr at u points (kg/m**3). real rhostru(nx,ny,nz) ! Averaged rhostr at v points (kg/m**3). real rhostrv(nx,ny,nz) ! Averaged rhostr at w points (kg/m**3). real rhostrw(nx,ny,nz) ! Pressure gradient term in w-eq. real wpgrad(nx,ny,nz) ! !####################################################################### ! ! Temporary WORK ARRAYS: ! !####################################################################### ! ! Pressure gradient term in u-eq. real upgrad(nx,ny,nz) ! Pressure gradient term in v-eq. real vpgrad(nx,ny,nz) ! Temporary work array real tem1 (nx,ny,nz) ! Temporary work array real tem2 (nx,ny,nz) ! Temporary work array real tem3 (nx,ny,nz) ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! ! Local curttim at a small time step real curtsml ! Local small time step real dtsml1 integer i, j, k integer ebc1,wbc1,nbc1,sbc1 ! !####################################################################### ! ! Include files: ! !####################################################################### ! include 'globcst.inc' include 'bndry.inc' include 'exbc.inc' ! !####################################################################### ! ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! !####################################################################### ! ! Divergence damping is activated if divdmp = 1. ! Otherwise, this portion of the code is skipped to save ! computations. ! !####################################################################### IF( divdmp.eq.1 .or. divdmp.eq.2 ) THEN ! !####################################################################### ! ! The divergence damping terms are defined as ! ! d(cdvdmph*div)/dx for u eq., ! d(cdvdmph*div)/dy for v eq., ! d(cdvdmpv*div)/dz for w eq., ! ! where ! ! div = (difx(u*rhostr) + dify(v*rhostr) + difz(wcont*rhostr))/j3 ! !####################################################################### ! ! !####################################################################### ! ! Compute the momentum divergence term, defined as ! ! div = d(u*rhostr)/dx + d(v*rhostr)/dy + d(wcont*rhostr)/dz. ! ! and combine the pressure and divergence damping into one ! array so that the pressure gradient force and divergence damping ! can be calculated in a single step. ! !####################################################################### ! DO 10 k=1,nz-1 DO 10 j=1,ny-1 DO 10 i=1,nx upgrad(i,j,k)=u(i,j,k)*rhostru(i,j,k)*mapfct(i,j,5) 10 CONTINUE DO 15 k=1,nz-1 DO 15 j=1,ny DO 15 i=1,nx-1 vpgrad(i,j,k)=v(i,j,k)*rhostrv(i,j,k)*mapfct(i,j,6) 15 CONTINUE DO 20 k=1,nz DO 20 j=1,ny-1 DO 20 i=1,nx-1 tem1(i,j,k)=wcont(i,j,k)*rhostrw(i,j,k) 20 CONTINUE DO 150 k=1,nz-1 DO 150 j=1,ny-1 DO 150 i=1,nx-1 tem3(i,j,k) = j3inv(i,j,k) * ( mapfct(i,j,1)*mapfct(i,j,1)& * ((upgrad(i+1,j,k)-upgrad(i,j,k))*dxinv +(vpgrad(i,j+1,k)& -vpgrad(i,j,k))*dyinv) + (tem1(i,j,k+1)-tem1(i,j,k))* & dzinv ) 150 CONTINUE ELSE DO 350 k=1,nz-1 DO 350 j=1,ny-1 DO 350 i=1,nx-1 tem3(i,j,k)= 0.0 350 CONTINUE ENDIF ! !####################################################################### ! ! Compute the pressure gradient force terms for the three ! momentum equations and store the values in upgrad, vpgrad, ! and wpgrad. ! ! When divergence damping is activated (divdmp=1), these arrays also ! contain the divergence damping terms. ! !####################################################################### ! CALL pgrad(nx,ny,nz, pprt, tem3, & & j1,j2,j3, upgrad,vpgrad,wpgrad,tem1,tem2) ! !####################################################################### ! ! Average rhofct at u points, stored in tem1 ! v points, stored in tem2 ! !####################################################################### ! CALL avgsu(rhofct,nx,ny,nz, 1,ny-1, 1,nz-1, tem1) CALL avgsv(rhofct,nx,ny,nz, 1,nx-1, 1,nz-1, tem2) ! !####################################################################### ! ! Integrate the u, and v equations forward by one small ! timestep dtsml1 using a forward-in-time integration scheme ! (that is, forward relative to the pressure gradient terms). ! !####################################################################### ! IF( nestgrd.eq.1 .and. mgrid.ne.1 ) THEN DO 450 k=2,nz-2 DO 450 j=2,ny-2 DO 450 i=2,nx-1 u(i,j,k)=u(i,j,k) +dtsml1*(uforce(i,j,k) -mapfct(i,j,2)* & tem1(i,j,k)*upgrad(i,j,k)) /rhostru(i,j,k) 450 CONTINUE DO 500 k=2,nz-2 DO 500 j=2,ny-1 DO 500 i=2,nx-2 v(i,j,k)=v(i,j,k) +dtsml1*(vforce(i,j,k) -mapfct(i,j,3)* & tem2(i,j,k)*vpgrad(i,j,k)) /rhostrv(i,j,k) 500 CONTINUE ! !####################################################################### ! ! Set the lateral boundary conditions for u and v given the ! time tendencies in the case of nested inner grid. ! !####################################################################### ! DO 100 k=2,nz-2 DO 100 j=1,ny-1 u( 1,j,k)=u( 1,j,k)+dtsml1 * udtwb(j,k) u(nx,j,k)=u(nx,j,k)+dtsml1 * udteb(j,k) 100 CONTINUE DO 110 k=2,nz-2 DO 110 i=2,nx-1 u(i, 1,k)=u(i, 1,k)+dtsml1 * udtsb(i,k) u(i,ny-1,k)=u(i,ny-1,k)+dtsml1 * udtnb(i,k) 110 CONTINUE DO 200 k=2,nz-2 DO 200 j=1,ny v( 1,j,k)=v( 1,j,k)+dtsml1 * vdtwb(j,k) v(nx-1,j,k)=v(nx-1,j,k)+dtsml1 * vdteb(j,k) 200 CONTINUE DO 210 k=2,nz-2 DO 210 i=2,nx-2 v(i, 1,k)=v(i, 1,k)+dtsml1 * vdtsb(i,k) v(i,ny,k)=v(i,ny,k)+dtsml1 * vdtnb(i,k) 210 CONTINUE ! !####################################################################### ! ! Set the top and bottom boundary conditions for u and v. ! !####################################################################### ! CALL bcu(nx,ny,nz, dtsml1, u, udteb,udtwb,udtnb,udtsb, 0,0,0,0 ,& tbc,bbc) CALL bcv(nx,ny,nz, dtsml1, v, vdteb,vdtwb,vdtnb,vdtsb, 0,0,0,0 ,& tbc,bbc) ELSE DO 452 k=2,nz-2 DO 452 j=1,ny-1 DO 452 i=2,nx-1 u(i,j,k)=u(i,j,k) +dtsml1*(uforce(i,j,k) -mapfct(i,j,2)* & tem1(i,j,k)*upgrad(i,j,k)) /rhostru(i,j,k) 452 CONTINUE DO 502 k=2,nz-2 DO 502 j=2,ny-1 DO 502 i=1,nx-1 v(i,j,k)=v(i,j,k) +dtsml1*(vforce(i,j,k) -mapfct(i,j,3)* & tem2(i,j,k)*vpgrad(i,j,k)) /rhostrv(i,j,k) 502 CONTINUE ! !####################################################################### ! ! Set the boundary conditions for u and v. ! !####################################################################### ! 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 bcu(nx,ny,nz, dtsml1, u, udteb,udtwb,udtnb,udtsb, & ebc,wbc,nbc1,sbc1,tbc,bbc) CALL bcv(nx,ny,nz, dtsml1, v, vdteb,vdtwb,vdtnb,vdtsb, & ebc1,wbc1,nbc,sbc,tbc,bbc) ELSE CALL bcu(nx,ny,nz, dtsml1, u, udteb,udtwb,udtnb,udtsb, & 0,0,0,0,tbc,bbc) CALL bcv(nx,ny,nz, dtsml1, v, vdteb,vdtwb,vdtnb,vdtsb, & 0,0,0,0,tbc,bbc) CALL exbcuv(nx,ny,nz, curtsml, u,v, exbcbuf(nu0exb), exbcbuf( & nv0exb), exbcbuf(nudtexb),exbcbuf(nvdtexb)) ENDIF ENDIF RETURN END ! ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE SOLVWPEX ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! SUBROUTINE SOLVWPEX(nx,ny,nz,exbcbufsz, dtsml1, curtsml, & & u,v,w,wcont,ptprt,pprt,phydro, & & wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, & & rhostr,ptbar,pbar,csndsq, wforce,wpgrad,pforce, & & x,y,z,zp, mapfct, j1,j2,j3,j3inv, & & rhostru,rhostrv,rhostrw, & & exbcbuf,rhofct, & & div,pdiv,tem1,tem2,tem3) ! !####################################################################### ! ! PURPOSE: ! ! Coordinate the time stepping of the w momentum and pressure equati ! using explicit time integration scheme. Divergence damping and the ! acoustically inactive forcing terms in w and pressure equations ! have been evaluated prior to this subroutine and are passed into ! this routine. ! !####################################################################### ! ! AUTHOR: Ming Xue ! 10/10/91 ! ! MODIFICATION HISTORY: ! ! 5/20/92 (K. Droegemeier and M. Xue) ! Added full documentation. ! ! 6/07/92 ( M. Xue) ! Now only tfuture fields of pprt, u, v and w are passed in. ! ! 2/10/93 (K. Droegemeier) ! Cleaned up documentation. ! ! 4/10/93 (M. Xue & Hao Jin) ! Add the terrain. ! ! 9/6/94 (M.Xue) ! Pressure perturbation buoyancy and base state pressure ! advection terms moved into the small time steps. ! ! 9/1/94 (Y. Lu) ! Cleaned up documentation. ! ! 01/28/95 (G. Bassett) ! Option added to turn off buoyancy terms (for buoyopt=0). ! ! 1/25/96 (Donghai Wang & Yuhe Liu) ! Added the map projection factor to ARPS governing equations. ! ! 4/27/1996 (M. Xue) ! Added code for peqopt=2 case. ! ! 5/6/96 (M. Xue) ! Replaced csndsq in pressure buoyancy term by cpdcv*pbar/rhobar. ! ! 10/16/97 (Donghai Wang) ! Using total density (rho) for the calculation of the pressure ! gradient force terms. ! ! 10/16/97 (Donghai Wang) ! Added the second order terms in the linerized buoyancy terms. ! ! 11/05/97 (D. Weber) ! Added phydro array for use in the bottom boundary condition for ! perturbation pressure (hydrostatic). ! ! 11/06/97 (D. Weber) ! Added three additional levels to the mapfct array. The three ! levels (4,5,6) represent the inverse of the first three in order. ! The inverse map factors are computed to improve efficiency. ! !####################################################################### ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! ! dtsml1 Local value of small time step size ! curtsml Local curttim at a small time step ! ! u x component of velocity at tfuture (m/s) ! v y component of velocity at tfuture (m/s) ! w Vertical component of Cartesian velocity at tfuture ! (m/s) ! wcont Contravariant vertical velocity (m/s) ! ptprt Perturbation potential temperature at all time levels ! (K) ! pprt Perturbation pressure at tfuture (Pascal) ! ! phydro Big time step forcing term for use in computing the ! hydrostatic pressure at k=1. ! ! wdteb Time tendency of the w field at the east boundary ! wdtwb Time tendency of the w field at the west boundary ! wdtnb Time tendency of the w field at the north boundary ! wdtsb Time tendency of the w field at the south boundary ! ! pdteb Time tendency of the pressure field at the east ! boundary ! pdtwb Time tendency of the pressure field at the west ! boundary ! pdtnb Time tendency of the pressure field at the north ! boundary ! pdtsb Time tendency of the pressure field at the south ! boundary ! ! rhostr Base state density rhobar times j3 (kg/m**3) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! csndsq Sound wave speed squared. ! ! wforce Acoustically inactive forcing terms in w-eq. ! (kg/(m*s)**2) ! pforce Acoustically inactive terms in pressure eq. (Pascal/s) ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp Vertical coordinate of grid points in physical space(m) ! ! mapfct Map factors at scalar, u and v points ! ! j1 Coordinate transformation Jacobian -d(zp)/dx ! j2 Coordinate transformation Jacobian -d(zp)/dy ! j3 Coordinate transformation Jacobian d(zp)/dz ! ! rhostru Average rhostr at u points (kg/m**3). ! rhostrv Average rhostr at v points (kg/m**3). ! rhostrw Average rhostr at w points (kg/m**3). ! ! OUTPUT: ! ! pprt Perturbation pressure at tfuture updated in time ! by stsml1 (Pascal) ! w Vertical component of Cartesian velocity at tfuture ! updated in time by dtsml1 (m/s) ! ! WORK ARRAYS: ! ! wpgrad Pressure gradient term in w-eq. (kg/(m*s)**2) ! div Velocity divergence (1/s), a local array ! pdiv Divergence term in pressure eq. ! (Pascal/s), a local array ! tem1 Temporary work array ! tem2 Temporary work array ! tem3 Temporary work array ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations: ! !####################################################################### ! ! Force explicit declarations implicit none ! Number of grid points in 3 directions integer nx, ny, nz ! Local curttim at a small time step real curtsml ! Local small time step size (s) real dtsml1 ! u-velocity at tfuture (m/s) real u (nx,ny,nz) ! v-velocity at tfuture (m/s) real v (nx,ny,nz) ! w-velocity at tfuture (m/s) real w (nx,ny,nz) ! Contravariant vertical velocity (m/s) real wcont (nx,ny,nz) ! Perturbation potential temperature (K) real ptprt (nx,ny,nz) ! Perturbation pressure at tfuture real pprt (nx,ny,nz) ! (Pascal) ! Big time step forcing for computing real phydro(nx,ny) ! hydrostatic pprt at k=1. ! Time tendency of w at east boundary real wdteb (ny,nz) ! Time tendency of w at west boundary real wdtwb (ny,nz) ! Time tendency of w at north boundary real wdtnb (nx,nz) ! Time tendency of w at south boundary real wdtsb (nx,nz) ! Time tendency of pressure at east real pdteb (ny,nz) ! boundary ! Time tendency of pressure at west real pdtwb (ny,nz) ! boundary ! Time tendency of pressure at north real pdtnb (nx,nz) ! boundary ! Time tendency of pressure at south real pdtsb (nx,nz) ! boundary ! Base state potential temperature (K) real ptbar (nx,ny,nz) ! Base state pressure (Pascal). real pbar (nx,ny,nz) ! Sound wave speed squared. real csndsq(nx,ny,nz) ! Base state density rhobar times j3. real rhostr(nx,ny,nz) ! Acoustically inactive forcing terms real wforce(nx,ny,nz) ! in w-momentum equation (kg/(m*s)**2) ! Acoustically inactive forcing terms real pforce(nx,ny,nz) ! in the pressure equation (Pascal/s) ! x-coord. of the physical and compu- real x (nx) ! tational grid. Defined at u-point. ! y-coord. of the physical and compu- real y (ny) ! tational grid. Defined at v-point. ! z-coord. of the computational grid. real z (nz) ! Defined at w-point on the staggered ! grid. ! Physical height coordinate defined at real zp (nx,ny,nz) ! w-point of the staggered grid. ! Map factors at scalar, u and v points real mapfct(nx,ny,6) ! Coordinate transformation Jacobian real j1 (nx,ny,nz) ! defined as - d( zp )/d( x ). ! Coordinate transformation Jacobian real j2 (nx,ny,nz) ! defined as - d( zp )/d( y ). ! Coordinate transformation Jacobian real j3 (nx,ny,nz) ! defined as d( zp )/d( z ). ! Coordinate transformation Jacobian real j3inv (nx,ny,nz) ! defined as d( zp )/d( z ). ! EXBC buffer size integer exbcbufsz ! EXBC buffer array real exbcbuf( exbcbufsz ) ! rho-factor:rhobar/rho real rhofct(nx,ny,nz) ! Averaged rhostr at u points (kg/m**3). real rhostru(nx,ny,nz) ! Averaged rhostr at v points (kg/m**3). real rhostrv(nx,ny,nz) ! Averaged rhostr at w points (kg/m**3). real rhostrw(nx,ny,nz) ! !####################################################################### ! ! Temporary work arrays: ! !####################################################################### ! ! Pressure gradient term in w-eq. real wpgrad(nx,ny,nz) ! Velocity divergence (1/s) real div (nx,ny,nz) ! Divergence term in pressure eq. real pdiv (nx,ny,nz) ! (Pascal/s) ! Temporary work array real tem1 (nx,ny,nz) ! Temporary work array real tem2 (nx,ny,nz) ! Temporary work array real tem3 (nx,ny,nz) ! !####################################################################### ! ! Include files: ! !####################################################################### ! include 'bndry.inc' include 'exbc.inc' include 'globcst.inc' include 'phycst.inc' ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i, j, k integer ebc1,wbc1,nbc1,sbc1 ! Switch for pprt/(rhobar*csndsq) term in w-eq. integer wpprt ! Switch for rhobar*g*w term in w-eq. integer prgw real prgwg5, g05 real pttem,ptem,tem ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! g05 = g*0.5 IF( bsnesq .eq.1 ) THEN ! Switch for pprt/(rhobar*csndsq) term in w-eq. wpprt = 0 ! Switch for rhobar*g*w term in w-eq. prgw = 0 ELSE ! Switch for pprt/(rhobar*csndsq) term in w-eq. wpprt = 1 ! Switch for rhobar*g*w term in w-eq. prgw = 1 ENDIF ! !####################################################################### ! ! The contribution from ptprt to the buoyancy term is calculated ! inside the small time steps if the potential temperature equation ! is solved inside small time steps, i.e., if ptsmlstp=1. ! ! The contribution from pressure perturbation to the buoyancy is ! always calculated inside the small time steps for better ! computational stability ! ! If buoyopt = 0 then turn off all buoyancy. ! !####################################################################### ! tem = 0.5*(1.0-cpdcv) IF ( buoyopt .ne. 0 ) THEN IF ( ptsmlstp.eq.1 ) THEN !1st-order IF ( buoy2nd .eq. 0) THEN DO 101 k=1,nz-1 DO 101 j=1,ny-1 DO 101 i=1,nx-1 tem2(i,j,k)=rhostr(i,j,k)*( ptprt(i,j,k)/ptbar(i,j,k) & -wpprt*pprt(i,j,k)/(cpdcv*pbar(i,j,k)) ) 101 CONTINUE !2nd-order ELSE DO 106 k=1,nz-1 DO 106 j=1,ny-1 DO 106 i=1,nx-1 pttem = ptprt(i,j,k)/ptbar(i,j,k) ptem = wpprt*pprt(i,j,k)/(cpdcv*pbar(i,j,k)) tem2(i,j,k)=rhostr(i,j,k)*( pttem-pttem*pttem-ptem- & tem*ptem*ptem + 0.5*pttem*ptem) 106 CONTINUE ENDIF ELSE !1st-order IF (buoy2nd .eq. 0) THEN DO 102 k=1,nz-1 DO 102 j=1,ny-1 DO 102 i=1,nx-1 tem2(i,j,k)= -wpprt*pprt(i,j,k)*rhostr(i,j,k)/(cpdcv* & pbar(i,j,k)) 102 CONTINUE !2nd-order ELSE DO 107 k=1,nz-1 DO 107 j=1,ny-1 DO 107 i=1,nx-1 ptem = wpprt*pprt(i,j,k)/(cpdcv*pbar(i,j,k)) tem2(i,j,k)=-rhostr(i,j,k)*(ptem+tem*ptem*ptem ) 107 CONTINUE ENDIF ENDIF DO 103 k=2,nz-1 DO 103 j=1,ny-1 DO 103 i=1,nx-1 tem1(i,j,k)=(tem2(i,j,k)+tem2(i,j,k-1))*g05 103 CONTINUE ELSE DO 104 k=1,nz-1 DO 104 j=1,ny-1 DO 104 i=1,nx-1 tem1(i,j,k) = 0.0 104 CONTINUE ENDIF ! !####################################################################### ! ! To assume mirror symmetry about the top and bottom boundary, ! the density/temperature has zero gradient across the boundary ! but g is assumed to change sign, so that the buoyancy is zero ! at the boundary. c.f. subroutine BUOYCY. ! !####################################################################### ! DO 115 j=1,ny-1 DO 115 i=1,nx-1 tem1(i,j,2)=0.0 tem1(i,j,nz-1)=0.0 115 CONTINUE ! !####################################################################### ! ! Average rhofct at w points, stored in tem2 ! !####################################################################### ! CALL avgsw(rhofct,nx,ny,nz, 1,nx-1, 1,ny-1,tem2) ! DO 120 k=2,nz-1 DO 120 j=1,ny-1 DO 120 i=1,nx-1 wpgrad(i,j,k)=tem2(i,j,k)*(wpgrad(i,j,k)-tem1(i,j,k) ) 120 CONTINUE ! !####################################################################### ! ! Integrate w equations forward by one small timestep dtsml1 using ! a forward-in-time integration scheme (that is, forward relative ! to the pressure gradient terms). ! !####################################################################### ! IF( nestgrd.eq.1 .and. mgrid.ne.1 ) THEN DO 500 k=2,nz-1 DO 500 j=2,ny-2 DO 500 i=2,nx-2 w(i,j,k)=w(i,j,k) +dtsml1 *(wforce(i,j,k)-wpgrad(i,j,k))/ & rhostrw(i,j,k) 500 CONTINUE ! !####################################################################### ! ! Set the lateral boundary conditions for w given the time tendencie ! in the case of nested inner grid. ! !####################################################################### ! DO 100 k=3,nz-2 DO 100 j=1,ny-1 w( 1,j,k)=w( 1,j,k)+dtsml1 * wdtwb(j,k) w(nx-1,j,k)=w(nx-1,j,k)+dtsml1 * wdteb(j,k) 100 CONTINUE DO 110 k=3,nz-2 DO 110 i=2,nx-2 w(i, 1,k)=w(i, 1,k)+dtsml1 * wdtsb(i,k) w(i,ny-1,k)=w(i,ny-1,k)+dtsml1 * wdtnb(i,k) 110 CONTINUE ELSE DO 501 k=2,nz-1 DO 501 j=1,ny-1 DO 501 i=1,nx-1 w(i,j,k)=w(i,j,k) +dtsml1 *(wforce(i,j,k)-wpgrad(i,j,k))/ & rhostrw(i,j,k) 501 CONTINUE ! !####################################################################### ! ! Apply the lateral boundary conditions for w. ! ! For the open boundary case, w at the lateral boundarues are solved ! from w equation directly therefore does not need to be reset. ! !####################################################################### ! ! Internal boundary conditions 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 lbcw(nx,ny,nz,dtsml,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, & ebc1,wbc1,nbc1,sbc1) ! Apply zero gradient condition ELSE CALL lbcw(nx,ny,nz,dtsml,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, & 3,3,3,3) ENDIF ENDIF ! !####################################################################### ! ! Calculate wcont at time tfuture, including the boundary ! points. Wcont at the lateral boundaries is calculated ! from boundary u, v and w. Wcont at the top and bottom ! depends on the boundary condition option. ! !####################################################################### ! CALL wcontra(nx,ny,nz,u,v,w,mapfct,j1,j2,j3, & & rhostr,rhostru,rhostrv,rhostrw,wcont,tem1,tem2) ! !####################################################################### ! ! Set the top and bottom boundary conditions for w based on u, v and ! wcont at the top and bottom boundaries. ! !####################################################################### ! CALL vbcw(nx,ny,nz,w,wcont,tbc,bbc,u,v, & & rhostr,rhostru,rhostrv,rhostrw, & & j1,j2,j3) IF( peqopt.eq.1) THEN ! !####################################################################### ! ! Average j3 to velocity points. ! !####################################################################### ! CALL avgsw(j3,nx,ny,nz, 1,nx-1, 1,ny-1, tem3) CALL avgsv(j3,nx,ny,nz, 1,nx-1, 1,nz-1, tem2) CALL avgsu(j3,nx,ny,nz, 1,ny-1, 1,nz-1, tem1) !####################################################################### ! ! Calculate the velocity divergence using newly updated velocity. ! !####################################################################### ! DO 260 k=2,nz-2 DO 260 j=1,ny-1 DO 260 i=1,nx-1 div(i,j,k) = mapfct(i,j,1)*mapfct(i,j,1) *(( u(i+1,j,k)* & tem1(i+1,j,k)*mapfct(i+1,j,5) - u(i ,j,k)*tem1(i ,j,k)* & mapfct(i ,j,5))*dxinv +( v(i,j+1,k)*tem2(i,j+1,k)*mapfct( & i,j+1,6) - v(i,j ,k)*tem2(i,j ,k)*mapfct(i,j ,6))*dyinv) & +( wcont(i,j,k+1)*tem3(i,j,k+1) - wcont(i,j,k)*tem3(i,j,k)& ) * dzinv 260 CONTINUE ! !####################################################################### ! ! Compute the divergence term and the base state pressure advection ! in the pressure equation. ! ! The pressure divergence term is: -rhostr*csndsq*div/j3 ! The base state pressure advection is -j3*w*d(pbar)/dzp = w*rhostr* ! ! These two terms are saved in pdiv. ! !####################################################################### ! prgwg5=prgw * g05 DO 265 k=2,nz-2 DO 265 j=1,ny-1 DO 265 i=1,nx-1 pdiv(i,j,k)= -csndsq(i,j,k)*rhostr(i,j,k)*div(i,j,k)* & j3inv(i,j,k) +prgwg5*(w(i,j,k)+w(i,j,k+1))*rhostr(i,j,k) 265 CONTINUE ELSE !####################################################################### ! ! Calculate the density weighted mass divergence using newly ! updated velocity. Then calculate the pressure divergence term: ! -csndsq*div ! !####################################################################### ! DO 10 k= 2,nz-2 DO 10 j= 1,ny-1 DO 10 i= 1,nx tem1(i,j,k)=u(i,j,k)*rhostru(i,j,k)*mapfct(i,j,5) 10 CONTINUE DO 15 k= 2,nz-2 DO 15 j= 1,ny DO 15 i= 1,nx-1 tem2(i,j,k)=v(i,j,k)*rhostrv(i,j,k)*mapfct(i,j,6) 15 CONTINUE DO 16 k= 2,nz-1 DO 16 j= 1,ny-1 DO 16 i= 1,nx-1 tem3(i,j,k)=wcont(i,j,k)*rhostrw(i,j,k) 16 CONTINUE DO 262 k=2,nz-2 DO 262 j=1,ny-1 DO 262 i=1,nx-1 pdiv(i,j,k)= -csndsq(i,j,k)*( mapfct(i,j,1)*mapfct(i,j,1) & *((tem1(i+1,j,k)-tem1(i,j,k))*dxinv + (tem2(i,j+1,k)-tem2(& i,j,k))*dyinv) + (tem3(i,j,k+1)-tem3(i,j,k))*dzinv ) 262 CONTINUE ENDIF ! !####################################################################### ! ! Integrate forward the pressure equation by one small timestep, ! using a backward (relative to the divergence term) integration ! scheme. ! ! And set lateral boundary conditions for the pressure equation. ! !####################################################################### ! IF( nestgrd.eq.1 .and. mgrid.ne.1 ) THEN DO 600 k=2,nz-2 DO 600 j=2,ny-2 DO 600 i=2,nx-2 pprt(i,j,k)=pprt(i,j,k) +dtsml1 *(pforce(i,j,k)+pdiv(i,j, & k))*j3inv(i,j,k) 600 CONTINUE DO 300 k=2,nz-2 DO 300 j=1,ny-1 pprt( 1,j,k)=pprt( 1,j,k)+dtsml1 * pdtwb(j,k) pprt(nx-1,j,k)=pprt(nx-1,j,k)+dtsml1 * pdteb(j,k) 300 CONTINUE DO 310 k=2,nz-2 DO 310 i=2,nx-2 pprt(i, 1,k)=pprt(i, 1,k)+dtsml1 * pdtsb(i,k) pprt(i,ny-1,k)=pprt(i,ny-1,k)+dtsml1 * pdtnb(i,k) 310 CONTINUE ! ! Call the pprt bottom boundary condition subroutine to ! compute the hydrostatic pprt at k=1. ! CALL pprtbbc(nx,ny,nz,g05,buoy2nd,rhostr,pprt,ptprt, pbar,ptbar,& phydro, tem1) ! new pprt at k=1. ! !####################################################################### ! ! And top and bottom boundary conditions for the pressure equation. ! !####################################################################### ! CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, & tem1(1,1,1),0,0,0,0 ,tbc,bbc) ! Not nested grid ELSE DO 710 k=2,nz-2 DO 710 j=1,ny-1 DO 710 i=1,nx-1 pprt(i,j,k)=pprt(i,j,k) +dtsml1 *(pforce(i,j,k)+pdiv(i,j, & k))*j3inv(i,j,k) 710 CONTINUE ! ! Call the pprt bottom boundary condition subroutine to ! compute the hydrostatic pprt at k=1. ! CALL pprtbbc(nx,ny,nz,g05,buoy2nd,rhostr,pprt,ptprt, pbar,ptbar,& phydro, tem1) ! new pprt at k=1. ! !####################################################################### ! ! Apply the boundary conditions for the pressure equation. ! ! For the open boundary case, the boundary value is solved from the ! equation. ! !####################################################################### ! ! Internal boundary conditions 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 bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, & tem1(1,1,1), ebc1,wbc1,nbc1,sbc1,tbc,bbc) ! External boundary condition ELSE CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, & tem1(1,1,1), 0,0,0,0,tbc,bbc) CALL exbcp(nx,ny,nz, curtsml, pprt, exbcbuf(npr0exb),exbcbuf( & nprdtexb)) ENDIF ENDIF RETURN END ! ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE SOLVWPIM ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! SUBROUTINE SOLVWPIM(nx,ny,nz,exbcbufsz, dtsml1, curtsml, & & u,v,w,wcont,ptprt,pprt,phydro, & & wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, & & rhostr,ptbar,pbar,csndsq, wforce,wpgrad,pforce, & & x,y,z,zp, mapfct, j1,j2,j3,j3inv, & & trigs1,trigs2,ifax1,ifax2, & & wsave1,wsave2,vwork1,vwork2, & & rhostru,rhostrv,rhostrw, & & exbcbuf,rhofct, & & fw,fp,wcontuv,tem1,tem2,tem3,tem4) ! !####################################################################### ! ! PURPOSE: ! ! Perform the time integration of w and pressure equations using ! implicit schem in vertical direction. The acoustically inactive ! forcing terms in these equations have been evaluated prior to this ! subroutine and are stored in wforce and pforce. ! !####################################################################### ! ! ! AUTHOR: M. Xue & H. Jin ! 8/20/93. ! ! 9/6/94 (M.Xue) ! A new version of vertically implicit small time step solver, ! The perturbation pressure contribution to buoyancy and the ! base state pressure advection terms are also treated implicitly ! in the vertical direction inside the small time steps. ! ! 9/1/94 (Y. Lu) ! Cleaned up documentation. ! ! 3/2/1995 (M. Xue and A. Shapiro) ! Significant bug fix in loop 600. rhostru there was mistakenly ! written as rhostrw. ! ! 10/31/95 (D. Weber) ! Added linear hydrostatic upper radiation condition for w-p. ! References are Klemp and Durran (JAS, 1983) and Chen (1991). ! Includes the use of trigs1,trigs2,ifax1,ifax2. ! ! 1/25/96 (Donghai Wang & Yuhe Liu) ! Added the map projection factor to ARPS governing equations. ! ! 4/27/1996 (M. Xue) ! Corrected the formulations of the coefficients of the tridiagonal ! equations. The oringal formulation was slightly inacurate. ! ! 5/6/96 (M. Xue) ! Replaced csndsq in pressure buoyancy term by cpdcv*pbar/rhobar. ! ! 3/11/1997 (M. Xue and D. Weber) ! Corrected a minor bug related to the radiation top boundary ! condition. ! ! 07/22/97 (D. Weber) ! Added even fft for linear hydrostatic w-p upper radiation conditio ! References are Klemp and Durran (JAS, 1983) and Chen (1991). ! Includes the use of wsave1,wsave2,vwork1,vwork2 (fftopt=2). ! ! 10/21/97 (Donghai Wang) ! Using total density (rho) in the calculation of the pressure ! gradient force terms. ! ! 10/21/97 (Donghai Wang) ! Added the second order terms in the linerized buoyancy terms. ! ! 11/05/97 (D. Weber) ! Added phydro array for use in the bottom boundary condition for ! perturbation pressure (hydrostatic). ! ! 11/06/97 (D. Weber) ! Added three additional levels to the mapfct array. The three ! levels (4,5,6) represent the inverse of the first three in order. ! The inverse map factors are computed to improve efficiency. ! Computed constants outside loops in selected loops. ! !####################################################################### ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! ! dtsml1 Local value of small time step size ! curtsml Local curttim at a small time step ! ! u x component of velocity at tfuture (m/s) ! v y component of velocity at tfuture (m/s) ! w Vertical component of Cartesian velocity at tfuture ! (m/s) ! wcont Contravariant vertical velocity (m/s) ! ptprt Perturbation potential temperature at time tpresent (K) ! pprt Perturbation pressure at time tpresent (Pascal) ! ! wdteb Time tendency of the w field at the east boundary ! wdtwb Time tendency of the w field at the west boundary ! wdtnb Time tendency of the w field at the north boundary ! wdtsb Time tendency of the w field at the south boundary ! ! pdteb Time tendency of the pressure field at the east ! boundary ! pdtwb Time tendency of the pressure field at the west ! boundary ! pdtnb Time tendency of the pressure field at the north ! boundary ! pdtsb Time tendency of the pressure field at the south ! boundary ! ! phydro Big time step forcing term for use in computing the ! hydrostatic pressure at k=1. ! ! rhostr Base state density rhobar times j3 (kg/m**3) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! csndsq Sound wave speed squared. ! ! wforce Acoustically inactive forcing terms in w-eq. ! (kg/(m*s)**2) ! wpgrad Pressure gradient term in w-eq. (kg/(m*s)**2) ! pforce Acoustically inactive terms in pressure eq. (Pascal/s) ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp Vertical coordinate of grid points in physical space(m) ! ! mapfct Map factors at scalar, u and v points ! ! j1 Coordinate transformation Jacobian -d(zp)/dx ! j2 Coordinate transformation Jacobian -d(zp)/dy ! j3 Coordinate transformation Jacobian d(zp)/dz ! trigs1 Array containing pre-computed trig function for ! fftopt=1. ! trigs2 Array containing pre-computed trig function for ! fftopt=1. ! ifax1 Array containing the factors of nx for fftopt=1. ! ifax2 Array containing the factors of ny for fftopt=1. ! ! vwork1 2-D work array for fftopt = 2. ! vwork2 2-D work array for fftopt = 2. ! wsave1 Work array for fftopt = 2. ! wsave2 Work array for fftopt = 2. ! ! rhostru Average rhostr at u points (kg/m**3). ! rhostrv Average rhostr at v points (kg/m**3). ! rhostrw Average rhostr at w points (kg/m**3). ! ! OUTPUT: ! ! pprt Perturbation pressure at tfuture updated in time ! by dtsml1 (Pascal) ! w Vertical component of Cartesian velocity at tfuture ! updated in time by dtsml1 (m/s) ! ! WORK ARRAYS: ! ! fw Work array to carry force terms in w-eqation. ! fp Work array to carry force terms in p-eqation. ! wcontuv Work array to carry the contributions of u & v to wcont ! tem1 Work array ! tem2 Work array ! tem3 Work array ! tem4 Work array ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations: ! !####################################################################### ! ! Force explicit declarations implicit none ! Number of grid points in 3 directions integer nx, ny, nz ! Local small time step size (s) real dtsml1 ! Local curttim at a small time step real curtsml ! u-velocity at tfuture (m/s) real u (nx,ny,nz) ! v-velocity at tfuture (m/s) real v (nx,ny,nz) ! w-velocity at tfuture (m/s) real w (nx,ny,nz) ! Contravariant vertical velocity (m/s) real wcont (nx,ny,nz) ! Perturbation potential temperature (K) real ptprt (nx,ny,nz) ! Perturbation pressure at tfuture real pprt (nx,ny,nz) ! (Pascal) ! Time tendency of w at east boundary real wdteb (ny,nz) ! Time tendency of w at west boundary real wdtwb (ny,nz) ! Time tendency of w at north boundary real wdtnb (nx,nz) ! Time tendency of w at south boundary real wdtsb (nx,nz) ! Time tendency of pressure at east real pdteb (ny,nz) ! boundary ! Time tendency of pressure at west real pdtwb (ny,nz) ! boundary ! Time tendency of pressure at north real pdtnb (nx,nz) ! boundary ! Time tendency of pressure at south real pdtsb (nx,nz) ! boundary ! Big time step forcing for computing real phydro(nx,ny) ! hydrostatic pprt at k=1. ! Base state density rhobar times j3. real rhostr(nx,ny,nz) ! Base state potential temperature (K) real ptbar (nx,ny,nz) ! Base state pressure (Pascal). real pbar (nx,ny,nz) ! Sound wave speed squared. real csndsq(nx,ny,nz) ! Acoustically inactive forcing terms real wforce(nx,ny,nz) ! in w-momentum equation (kg/(m*s)**2) ! Pressure gradient term in w-eq. real wpgrad(nx,ny,nz) ! (kg/(m*s)**2) ! Acoustically inactive forcing terms real pforce(nx,ny,nz) ! in the pressure equation (Pascal/s) ! x-coord. of the physical and compu- real x (nx) ! tational grid. Defined at u-point. ! y-coord. of the physical and compu- real y (ny) ! tational grid. Defined at v-point. ! z-coord. of the computational grid. real z (nz) ! Defined at w-point on the staggered ! grid. ! Physical height coordinate defined at real zp (nx,ny,nz) ! w-point of the staggered grid. ! Map factors at scalar, u and v points real mapfct(nx,ny,6) ! Coordinate transformation Jacobian real j1 (nx,ny,nz) ! defined as - d( zp )/d( x ). ! Coordinate transformation Jacobian real j2 (nx,ny,nz) ! defined as - d( zp )/d( y ). ! Coordinate transformation Jacobian real j3 (nx,ny,nz) ! defined as d( zp )/d( z ). ! Coordinate transformation Jacobian real j3inv (nx,ny,nz) ! defined as d( zp )/d( z ). ! EXBC buffer size integer exbcbufsz ! EXBC buffer array real exbcbuf( exbcbufsz ) ! rho-factor: rhobar/rho real rhofct(nx,ny,nz) ! Array containing pre-computed trig real trigs1(3*(nx-1)/2+1) ! function for fftopt=1. ! Array containing pre-computed trig real trigs2(3*(ny-1)/2+1) ! function for fftopt=1. ! Array containing the factors of nx for integer ifax1(13) ! fftopt=1. ! Array containing the factors of ny for integer ifax2(13) ! fftopt=1. ! 2-D work array for fftopt = 2. real vwork1 (nx+1,ny+1) ! 2-D work array for fftopt = 2. real vwork2 (ny,nx+1) ! Work array for fftopt = 2. real wsave1 (3*(ny-1)+15) ! Work array for fftopt = 2. real wsave2 (3*(nx-1)+15) ! Averaged rhostr at u points (kg/m**3). real rhostru(nx,ny,nz) ! Averaged rhostr at v points (kg/m**3). real rhostrv(nx,ny,nz) ! Averaged rhostr at w points (kg/m**3). real rhostrw(nx,ny,nz) ! !####################################################################### ! ! Temporary WORK ARRAYS: ! !####################################################################### ! ! Pressure gradient term in w-eq. real fp (nx,ny,nz) ! Force in w equation.(kg/(m*s)**2) real fw (nx,ny,nz) ! Contributions of u & v to wcont real wcontuv(nx,ny,nz) ! Temporary work 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) ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i, j, k integer ebc1,wbc1,nbc1,sbc1 real g05,pk,nk,g05wp,g05pr integer wpprt, prgw, itema real tema,temb,w1,w2,nrho !Switch for 1st-order or 2nd-order in buoyancy integer buoy2swt real ptemk,ptemk1,pttemk,pttemk1 ! !####################################################################### ! ! Include files: ! !####################################################################### ! include 'globcst.inc' include 'bndry.inc' include 'exbc.inc' ! Physical constants include 'phycst.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! IF( bsnesq .eq.1 ) THEN ! Switch for pprt/(rhobar*csndsq) term in w-eq. wpprt = 0 ! Switch for rhobar*g*w term in w-eq. prgw = 0 ELSE ! Switch for pprt/(rhobar*csndsq) term in w-eq. wpprt = 1 ! Switch for rhobar*g*w term in w-eq. prgw = 1 ENDIF g05 = g*0.5 g05wp = g05*wpprt g05pr = g05*prgw ! !####################################################################### ! ! Calculate the horizontal velocity divergence using newly updated ! u and v velocity plus half vertical divergence from wcont, and ! store the result in tem3. ! ! Namely, tem3 = difx(u*avgx(j3))+dify(v*avgy(j3))+ ! (1-tacoef)*difz(wcont*avgz(j3)) ! !####################################################################### ! CALL avgsu(j3,nx,ny,nz, 1,ny-1, 1,nz-1, tem1) CALL avgsv(j3,nx,ny,nz, 1,nx-1, 1,nz-1, tem2) DO 100 k=2,nz-2 DO 100 j=1,ny-1 DO 100 i=1,nx-1 tem3(i,j,k) = mapfct(i,j,1)*mapfct(i,j,1) * ( (u(i+1,j,k)* & tem1(i+1,j,k)*mapfct(i+1,j,5) -u(i,j,k)*tem1(i,j,k)*mapfct( & i,j,5))*dxinv + (v(i,j+1,k)*tem2(i,j+1,k)*mapfct(i,j+1,6) & -v(i,j,k)*tem2(i,j,k)*mapfct(i,j,6))*dyinv ) 100 CONTINUE CALL avgsw(j3,nx,ny,nz, 1,nx-1, 1,ny-1, tem1) tema = (1.0-tacoef)*dzinv DO 105 k=2,nz-2 DO 105 j=1,ny-1 DO 105 i=1,nx-1 tem3(i,j,k)= tem3(i,j,k) +(wcont(i,j,k+1)*tem1(i,j,k+1)- & wcont(i,j,k)*tem1(i,j,k)) *tema 105 CONTINUE ! !####################################################################### ! ! Compute the forcing terms in pressure equation ! ! fp = dtsml/j3 *(pforce-rhobar*csndsq*tem3+(1-tacoef)*rhostr*g*w) ! ! rhostr*g*w is the base state pressure advection term. ! !####################################################################### ! tema = g05pr*(1.0-tacoef) DO 170 k=2,nz-2 DO 170 j=1,ny-1 DO 170 i=1,nx-1 fp(i,j,k)=dtsml1*j3inv(i,j,k) *(pforce(i,j,k) -csndsq(i,j,k)& *rhostr(i,j,k)*j3inv(i,j,k)*tem3(i,j,k) +tema*(w(i,j,k)+w(i,& j,k+1))*rhostr(i,j,k) ) 170 CONTINUE ! !####################################################################### ! ! Compute the first two terms in the coutravariant vertical ! velocity formualation: ! ! wcontuv = (avgx(avgz(ustr)*j1)+avgy(avgz(vstr)*j2))/rhostrw ! !####################################################################### ! IF( ternopt.eq.0 ) THEN DO 240 k=2,nz-1 DO 240 j=1,ny-1 DO 240 i=1,nx-1 wcontuv(i,j,k)=0.0 240 CONTINUE ELSE DO 30 k= 1,nz-1 DO 30 j= 1,ny-1 DO 30 i= 1,nx tem1(i,j,k)=u(i,j,k)*rhostru(i,j,k) 30 CONTINUE DO 35 k= 1,nz-1 DO 35 j= 1,ny DO 35 i= 1,nx-1 tem2(i,j,k)=v(i,j,k)*rhostrv(i,j,k) 35 CONTINUE DO 40 k= 2,nz-1 DO 40 j= 1,ny-1 DO 40 i= 1,nx-1 wcontuv(i,j,k)=((tem1(i ,j,k)+tem1(i ,j,k-1))*j1(i & ,j,k) +(tem1(i+1,j,k)+tem1(i+1,j,k-1))*j1(i+1,j,k) & +(tem2(i ,j,k)+tem2(i ,j,k-1))*j2(i ,j,k) +(tem2(i,j+1,k)+& tem2(i,j+1,k-1))*j2(i,j+1,k)) *0.25 * mapfct(i,j,1) / & rhostrw(i,j,k) 40 CONTINUE DO 244 j=1,ny-1 DO 244 i=1,nx-1 wcontuv(i,j,nz-1)=0.0 wcontuv(i,j,2) = ((u(i ,j,2)+u(i ,j,1))*j1(i,j,2) +(u(i+1,j,& 2)+u(i+1,j,1))*j1(i+1,j,2) +(v(i,j ,2)+v(i,j ,1))*j2(i,j,2) & +(v(i,j+1,2)+v(i,j+1,1))*j2(i,j+1,2)) *0.25*mapfct(i,j,1) 244 CONTINUE ENDIF ! !####################################################################### ! ! Average rhofct at w points, stored in tem4 ! !####################################################################### ! CALL avgsw(rhofct,nx,ny,nz, 1,nx-1, 1,ny-1,tem4) ! !####################################################################### ! ! Compute the right-hand-side forcing term in tridiagonal linear ! equation. Array w is used to store this forcing term. ! ! First, we add the contribution of pressure perturbation to the ! buoyancy to wforce and wpgrad. ! ! This term has a form of - rhostr*g*pprt/(cpdcv*pbar), ! !####################################################################### ! tema = g05wp/cpdcv !1st-order IF ( buoy2nd .eq. 0) THEN !Switch for 1st-order or 2nd-order buoy2swt = 0 !2nd-order ELSE buoy2swt = 1 ENDIF temb = 0.5*buoy2swt*(1.0-cpdcv)/cpdcv DO 350 k=3,nz-2 DO 350 j=1,ny-1 DO 350 i=1,nx-1 ptemk = pprt(i,j,k )*rhostr(i,j,k )/pbar(i,j,k ) ptemk1= pprt(i,j,k-1)*rhostr(i,j,k-1)/pbar(i,j,k-1) fw(i,j,k) = wforce(i,j,k)-tem4(i,j,k)*(wpgrad(i,j,k) & +(ptemk+ptemk1+ temb*(ptemk*ptemk+ptemk1*ptemk1)) *tema) 350 CONTINUE ! !####################################################################### ! ! When potential temperature equation is solved within small time ! steps, the contribution from ptprt to buoyancy term is calculated ! here. ! !####################################################################### ! IF( ptsmlstp.eq.1 ) THEN tema = 1.0/(2.0*cpdcv) DO 352 k=3,nz-2 DO 352 j=1,ny-1 DO 352 i=1,nx-1 ptemk = pprt(i,j,k )*rhostr(i,j,k )/pbar(i,j,k ) ptemk1= pprt(i,j,k-1)*rhostr(i,j,k-1)/pbar(i,j,k-1) pttemk = ptprt(i,j,k )*rhostr(i,j,k )/ptbar(i,j,k ) pttemk1= ptprt(i,j,k-1)*rhostr(i,j,k-1)/ptbar(i,j,k-1) fw(i,j,k)=fw(i,j,k)+ (pttemk+pttemk1 -buoy2swt*(pttemk* & pttemk+pttemk1*pttemk1 +tema*(ptemk*pttemk + ptemk1* & pttemk1))) *g05*tem4(i,j,k) 352 CONTINUE ENDIF ! !####################################################################### ! ! tem3=fp-dtsml*rhostr*csndsq*tacoef* d(wcontuv)/dz /(j3*j3) ! !####################################################################### ! tema = tacoef*dtsml1*dzinv DO 400 k=2,nz-2 DO 400 j=1,ny-1 DO 400 i=1,nx-1 tem3(i,j,k)=fp(i,j,k) -tema*csndsq(i,j,k)*rhostr(i,j,k) & *(wcontuv(i,j,k+1)-wcontuv(i,j,k)) *(j3inv(i,j,k)*j3inv(i,j,& k)) 400 CONTINUE ! !####################################################################### ! ! fw=fw-tacoef*d(tem3)/dz-g*tacoef*avgz(rhostr*tem3/(cpdcv*pbar))) ! !####################################################################### ! tema = g05wp/cpdcv DO 420 k=3,nz-2 DO 420 j=1,ny-1 DO 420 i=1,nx-1 fw(i,j,k)=fw(i,j,k)-tacoef*tem4(i,j,k)*( (tem3(i,j,k)-tem3( & i,j,k-1))*dzinv + (rhostr(i,j,k )*tem3(i,j,k )/pbar(i,j,k & ) +rhostr(i,j,k-1)*tem3(i,j,k-1)/pbar(i,j,k-1))*tema) 420 CONTINUE ! !####################################################################### ! ! fw=fw*dtsml/rhostr + w ! !####################################################################### ! DO 425 k=3,nz-2 DO 425 j=1,ny-1 DO 425 i=1,nx-1 fw(i,j,k)=fw(i,j,k)*dtsml1/rhostrw(i,j,k) + w(i,j,k) 425 CONTINUE ! !####################################################################### ! ! Prepare the left-hand-side coefficents for tridiagonal ! equation set: ! ! A(k)*w(k-1)+B(k)*w(k)+C(k)*w(k+1)=D(k) for k=3,nz-2 ! ! w(2) and w(nz-1) are used as the boundary conditions. ! ! Due to the lack of work arrays, we are storing the coefficients ! A in rhostru, B in rhostrv, and C in rhostrw. D is stored in fw. ! ! rhostru, rhostrv and rhostrw are re-calculated from rhostr after ! they have been used. ! !####################################################################### ! DO 505 k=2,nz-2 DO 505 j=1,ny-1 DO 505 i=1,nx-1 tem1(i,j,k)= rhostr(i,j,k)*g05pr*j3inv(i,j,k) tem2(i,j,k)= rhostr(i,j,k)*csndsq(i,j,k)*j3inv(i,j,k)**2* & dzinv 505 CONTINUE tema = (dtsml1*tacoef)**2 * wpprt*g/cpdcv temb = (dtsml1*tacoef)**2 * dzinv DO 510 k=3,nz-2 DO 510 j=1,ny-1 DO 510 i=1,nx-1 pk = tem4(i,j,k)*tema/(pbar(i,j,k)+pbar(i,j,k-1)) nk = tem4(i,j,k)*temb/rhostrw(i,j,k) rhostru(i,j,k)= (-nk+pk)*(tem1(i,j,k-1)+tem2(i,j,k-1)) rhostrw(i,j,k)= ( nk+pk)*(tem1(i,j,k )-tem2(i,j,k )) rhostrv(i,j,k)= 1 +nk*(tem1(i,j,k)+tem2(i,j,k)-tem1(i,j,k-1)& +tem2(i,j,k-1)) +pk*(tem1(i,j,k)+tem2(i,j,k)+tem1(i,j,k-1)- & tem2(i,j,k-1)) 510 CONTINUE ! !####################################################################### ! ! Reset fw on the boundaries using the top and bottom boundary ! conditions of w. ! ! w=0.0 at k=nz-1. At the lower boundary, wcont=0.0, therefore ! w=-wcontuv, where wcontuv was calculated in loop 240. ! !####################################################################### ! DO 600 j=1,ny-1 DO 600 i=1,nx-1 fw(i,j,3) =fw(i,j,3)+wcontuv(i,j,2)*rhostru(i,j,3) fw(i,j,nz-2)=fw(i,j,nz-2)+0.0*rhostrw(i,j,nz-2) 600 CONTINUE ! !####################################################################### ! ! Call the tridiagonal solver for either a rigid (tbc=1) or open upp ! boundary condition for w (tbc=4). ! ! For tbc=4, we have a choice of fft's for the upper boundary: ! fftopt =1, periodic fft for linearized hydrostatic radiation condi ! fftopt =2, even fft for linearized hydrostatic radiation condition ! ! The references are Klemp and Durran Jas (1983) and Chen (MWR) 1991 ! ! w1 is the coef for w(nz-1) and w2 is the coef. for w(nz-2) in the ! pressure equation. The pressure equation is solved at p(nz-2) ! for w(nz-1). tem2(i,j,nz-2) is the summation of the known terms ! in the pressure equation. The only unknowns are p(nz-2), w(nz-1) ! and w(nz-1) at the new time level. In the tridiagonal solver ! the elimination step is performed and w(nz-2) is given in termsx ! of w(nz-1) and known terms. THe next step is to use the rad. ! condition (in fourier space): ! ! ! p = N * rhobar * w(nz-1) / abs(kx+ky) ! ! The end result is a relation for w(nz-1) given known quantities. ! Trigs1 and trigs2 are the predetermined trig functions used in th ! fft program in tridiag. ! !####################################################################### ! apply linear radiation condition to w at nz-1. IF(tbc.eq.4)THEN tema = rhostr(nx/2,ny/2,nz-2)*csndsq(nx/2,ny/2,nz-2)* j3inv(nx/ & 2,ny/2,nz-2) temb = dtsml1*j3inv(nx/2,ny/2,nz-2) w2=temb*tacoef*(tema*dzinv+g05*prgw*rhostr(nx/2,ny/2,nz-2)) w1=temb*tacoef*(-tema*dzinv+g05*prgw*rhostr(nx/2,ny/2,nz-2)) DO 605 j=1,ny-1 DO 605 i=1,nx-1 tem2(i,j,2)=pprt(i,j,nz-2)+tem3(i,j,nz-2) 605 CONTINUE nrho=sqrt(g/ptbar(nx/2,ny/2,nz-2)*(ptbar(nx/2,ny/2,nz-1) & -ptbar(nx/2,ny/2,nz-3))*dzinv*0.5) *rhostr(nx/2,ny/2,nz-2)* & j3inv(nx/2,ny/2,nz-2) END IF ! !####################################################################### ! ! NOTE: tem2(1,1,4) is passed into subroutine tridiag and becomes ! a 2-D array of size (nx+1,ny+1). ! !####################################################################### ! CALL tridiag(nx,ny,nz,rhostru,rhostrv,rhostrw,fw,tem1, & & tem2(1,1,1),tem2(1,1,2),w1,w2,nrho,tem2(1,1,3),trigs1,trigs2, & & ifax1,ifax2,wsave1,wsave2,vwork1,vwork2,tem2(1,1,4)) ! !####################################################################### ! ! Restore rhostru, rhostrv and rhostrw after they have been used ! are work arrays. ! !####################################################################### ! CALL rhouvw(nx,ny,nz,rhostr,rhostru,rhostrv,rhostrw) ! !####################################################################### ! ! On exit of tridiag, the interior solution of w is stored in fw. ! !####################################################################### ! set itema = nz-1 IF(tbc.eq.4)THEN itema = nz-1 ! set itema = nz-2 ELSE itema = nz-2 END IF IF( nestgrd.ne.1 .or. mgrid.eq.1 ) THEN ! For non-nesting case or the ! base-grid of the nested case DO 620 k=3,itema DO 620 j=1,ny-1 DO 620 i=1,nx-1 w(i,j,k) = fw(i,j,k) 620 CONTINUE ! Internal boundary condition 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 lbcw(nx,ny,nz,dtsml1,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, & ebc1,wbc1,nbc1,sbc1) ! External boundary condition ELSE CALL lbcw(nx,ny,nz,dtsml1,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, & 3,3,3,3) ENDIF ! For nested interior grid ELSE DO 640 k=3,itema DO 640 j=2,ny-2 DO 640 i=2,nx-2 w(i,j,k) = fw(i,j,k) 640 CONTINUE DO 700 k=3,nz-2 DO 700 j=1,ny-1 w( 1,j,k)=w( 1,j,k)+dtsml1 * wdtwb(j,k) w(nx-1,j,k)=w(nx-1,j,k)+dtsml1 * wdteb(j,k) 700 CONTINUE DO 710 k=3,nz-2 DO 710 i=2,nx-2 w(i, 1,k)=w(i, 1,k)+dtsml1 * wdtsb(i,k) w(i,ny-1,k)=w(i,ny-1,k)+dtsml1 * wdtnb(i,k) 710 CONTINUE ENDIF ! !####################################################################### ! ! Calculate wcont at time tfuture, including the boundary ! points. Wcont at the lateral boundaries is calculated ! from boundary u, v and w. Wcont at the top and bottom ! depends on the boundary condition option. ! !####################################################################### ! CALL wcontra(nx,ny,nz,u,v,w,mapfct,j1,j2,j3, & & rhostr,rhostru,rhostrv,rhostrw,wcont,tem1,tem2) ! !####################################################################### ! ! Set the top and bottom boundary conditions for w based on u, v and ! wcont at the top and bottom boundaries. ! !####################################################################### ! CALL vbcw(nx,ny,nz,w,wcont,tbc,bbc,u,v, & & rhostr,rhostru,rhostrv,rhostrw, & & j1,j2,j3) ! !####################################################################### ! ! Calculate the new pressure ! ! pprt(new) = pprt(old)+fp+dtsml*tacoef/j3* ! (rhostr*g*avg(w)-rhostr/j3*csndsq*difz(j3*wcont)) ! !####################################################################### ! DO 780 k=2,nz-1 DO 780 j=1,ny-1 DO 780 i=1,nx-1 tem1(i,j,k)=0.5*(j3(i,j,k)+j3(i,j,k-1)) 780 CONTINUE tema = tacoef*dtsml1 DO 790 k=2,nz-2 DO 790 j=1,ny-1 DO 790 i=1,nx-1 fp(i,j,k)=fp(i,j,k)+tema*j3inv(i,j,k)* (rhostr(i,j,k)*g05pr*& (w(i,j,k+1)+w(i,j,k)) -rhostr(i,j,k)*csndsq(i,j,k)* (tem1(i,& j,k+1)*wcont(i,j,k+1)-tem1(i,j,k)*wcont(i,j,k)) *(j3inv(i,j,& k)*dzinv) ) 790 CONTINUE IF( nestgrd.eq.1 .and. mgrid.ne.1 ) THEN ! For nested interior grid DO 800 k=2,nz-2 DO 800 j=2,ny-2 DO 800 i=2,nx-2 pprt(i,j,k)=pprt(i,j,k)+fp(i,j,k) 800 CONTINUE DO 900 k=2,nz-2 DO 900 j=1,ny-1 pprt( 1,j,k)=pprt( 1,j,k)+dtsml1 * pdtwb(j,k) pprt(nx-1,j,k)=pprt(nx-1,j,k)+dtsml1 * pdteb(j,k) 900 CONTINUE DO 910 k=2,nz-2 DO 910 i=2,nx-2 pprt(i, 1,k)=pprt(i, 1,k)+dtsml1 * pdtsb(i,k) pprt(i,ny-1,k)=pprt(i,ny-1,k)+dtsml1 * pdtnb(i,k) 910 CONTINUE ! ! Call the pprt bottom boundary condition subroutine to ! compute the hydrostatic pprt at k=1. ! CALL pprtbbc(nx,ny,nz,g05,buoy2swt,rhostr,pprt,ptprt, pbar, & ptbar,phydro, tem1) ! new pprt at k=1. CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, & tem1(1,1,1), 0,0,0,0,tbc,bbc) ! For non-nesting case or the ELSE ! base-grid of the nested case DO 810 k=2,nz-2 DO 810 j=1,ny-1 DO 810 i=1,nx-1 pprt(i,j,k)=pprt(i,j,k)+fp(i,j,k) 810 CONTINUE ! ! Call the pprt bottom boundary condition subroutine to ! compute the hydrostatic pprt at k=1. ! CALL pprtbbc(nx,ny,nz,g05,buoy2swt,rhostr,pprt,ptprt, pbar, & ptbar,phydro, tem1) ! new pprt at k=1. ! !####################################################################### ! ! Apply the boundary conditions for the pressure equation. ! ! For the open boundary case, the boundary value of pprt is ! predicted by the pressure equation. ! !####################################################################### ! ! Internal boundary conditions 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 bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, & tem1(1,1,1), ebc1,wbc1,nbc1,sbc1,tbc,bbc) ! External boundary condition ELSE CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, & tem1(1,1,1), 0,0,0,0,tbc,bbc) CALL exbcp(nx,ny,nz, curtsml, pprt, exbcbuf(npr0exb),exbcbuf( & nprdtexb)) ENDIF ENDIF RETURN END ! ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE SOLVWPIM1 ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! SUBROUTINE SOLVWPIM1(nx,ny,nz,exbcbufsz, dtsml1, curtsml, & & u,v,w,wcont,ptprt,pprt,phydro, & & wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, & & rhostr,ptbar,pbar,csndsq, wforce,wpgrad,pforce, & & x,y,z,zp, mapfct, j1,j2,j3,j3inv, & & trigs1,trigs2,ifax1,ifax2, & & wsave1,wsave2,vwork1,vwork2, & & rhostru,rhostrv,rhostrw, & & exbcbuf,rhofct, & & fw,fp,tem1,tem2,tem3,tem4,tem5) ! !####################################################################### ! ! PURPOSE: ! ! Perform the time integration of w and pressure equations using ! implicit scheme in vertical direction. The acoustically inactive ! forcing terms in these equations have been evaluated prior to this ! subroutine and are stored in wforce and pforce. ! ! This routine is for peqopt=2. ! !####################################################################### ! ! ! AUTHOR: M. Xue ! 4/26/1996. ! Written for peqopt=2 based on SOLVWPIM. ! ! MODIFICATION HISTORY: ! ! 07/22/97 (D. Weber) ! Added wsave1,wsave2,vwork1,vwork2 for use in the even fft version ! of the upper w-p radiation condition (fftopt=2). ! ! 10/21/97 (Donghai Wang) ! Using total density (rho) in the calculation of the pressure ! gradient force terms. ! ! 10/21/97 (Donghai Wang) ! Added the second order terms in the linerized buoyancy terms. ! ! 11/05/97 (D. Weber) ! Added phydro array for use in the bottom boundary condition for ! perturbation pressure (hydrostatic). ! ! 11/06/97 (D. Weber) ! Added three additional levels to the mapfct array. The three ! levels (4,5,6) represent the inverse of the first three in order. ! The inverse map factors are computed to improve efficiency. ! Removed multiplication of constants from loops. ! !####################################################################### ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! ! dtsml1 Local value of small time step size ! curtsml Local curttim at a small time step ! ! u x component of velocity at tfuture (m/s) ! v y component of velocity at tfuture (m/s) ! w Vertical component of Cartesian velocity at tfuture ! (m/s) ! wcont Contravariant vertical velocity (m/s) ! ptprt Perturbation potential temperature at time tpresent (K) ! pprt Perturbation pressure at time tpresent (Pascal) ! ! wdteb Time tendency of the w field at the east boundary ! wdtwb Time tendency of the w field at the west boundary ! wdtnb Time tendency of the w field at the north boundary ! wdtsb Time tendency of the w field at the south boundary ! ! pdteb Time tendency of the pressure field at the east ! boundary ! pdtwb Time tendency of the pressure field at the west ! boundary ! pdtnb Time tendency of the pressure field at the north ! boundary ! pdtsb Time tendency of the pressure field at the south ! boundary ! ! phydro Big time step forcing term for use in computing the ! hydrostatic pressure at k=1. ! ! rhostr Base state density rhobar times j3 (kg/m**3) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! csndsq Sound wave speed squared. ! ! wforce Acoustically inactive forcing terms in w-eq. ! (kg/(m*s)**2) ! wpgrad Pressure gradient term in w-eq. (kg/(m*s)**2) ! pforce Acoustically inactive terms in pressure eq. (Pascal/s) ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp Vertical coordinate of grid points in physical space(m) ! ! mapfct Map factors at scalar, u and v points ! ! j1 Coordinate transformation Jacobian -d(zp)/dx ! j2 Coordinate transformation Jacobian -d(zp)/dy ! j3 Coordinate transformation Jacobian d(zp)/dz ! trigs1 Array containing pre-computed trig function for ! fftopt=1. ! trigs2 Array containing pre-computed trig function for ! fftopt=1. ! ifax1 Array containing the factors of nx for fftopt=1. ! ifax2 Array containing the factors of ny for fftopt=1. ! ! vwork1 2-D work array for fftopt=2. ! vwork2 2-D work array for fftopt=2. ! wsave1 Work array for fftopt=2. ! wsave2 Work array for fftopt=2. ! ! rhostru Average rhostr at u points (kg/m**3). ! rhostrv Average rhostr at v points (kg/m**3). ! rhostrw Average rhostr at w points (kg/m**3). ! ! OUTPUT: ! ! pprt Perturbation pressure at tfuture updated in time ! by dtsml1 (Pascal) ! w Vertical component of Cartesian velocity at tfuture ! updated in time by dtsml1 (m/s) ! ! WORK ARRAYS: ! ! fw Work array to carry force terms in w-eqation. ! fp Work array to carry force terms in p-eqation. ! tem1 Work array ! tem2 Work array ! tem3 Work array ! tem4 Work array ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations: ! !####################################################################### ! ! Force explicit declarations implicit none ! Number of grid points in 3 directions integer nx, ny, nz ! Local small time step size (s) real dtsml1 ! Local curttim at a small time step real curtsml ! u-velocity at tfuture (m/s) real u (nx,ny,nz) ! v-velocity at tfuture (m/s) real v (nx,ny,nz) ! w-velocity at tfuture (m/s) real w (nx,ny,nz) ! Contravariant vertical velocity (m/s) real wcont (nx,ny,nz) ! Perturbation potential temperature (K) real ptprt (nx,ny,nz) ! Perturbation pressure at tfuture real pprt (nx,ny,nz) ! (Pascal) ! Time tendency of w at east boundary real wdteb (ny,nz) ! Time tendency of w at west boundary real wdtwb (ny,nz) ! Time tendency of w at north boundary real wdtnb (nx,nz) ! Time tendency of w at south boundary real wdtsb (nx,nz) ! Time tendency of pressure at east real pdteb (ny,nz) ! boundary ! Time tendency of pressure at west real pdtwb (ny,nz) ! boundary ! Time tendency of pressure at north real pdtnb (nx,nz) ! boundary ! Time tendency of pressure at south real pdtsb (nx,nz) ! boundary ! Big time step forcing for computing real phydro(nx,ny) ! hydrostatic pprt at k=1. ! Base state density rhobar times j3. real rhostr(nx,ny,nz) ! Base state potential temperature (K) real ptbar (nx,ny,nz) ! Base state pressure (Pascal). real pbar (nx,ny,nz) ! Sound wave speed squared. real csndsq(nx,ny,nz) ! Acoustically inactive forcing terms real wforce(nx,ny,nz) ! in w-momentum equation (kg/(m*s)**2) ! Pressure gradient term in w-eq. real wpgrad(nx,ny,nz) ! (kg/(m*s)**2) ! Acoustically inactive forcing terms real pforce(nx,ny,nz) ! in the pressure equation (Pascal/s) ! x-coord. of the physical and compu- real x (nx) ! tational grid. Defined at u-point. ! y-coord. of the physical and compu- real y (ny) ! tational grid. Defined at v-point. ! z-coord. of the computational grid. real z (nz) ! Defined at w-point on the staggered ! grid. ! Physical height coordinate defined at real zp (nx,ny,nz) ! w-point of the staggered grid. ! Map factors at scalar, u and v points real mapfct(nx,ny,6) ! Coordinate transformation Jacobian real j1 (nx,ny,nz) ! defined as - d( zp )/d( x ). ! Coordinate transformation Jacobian real j2 (nx,ny,nz) ! defined as - d( zp )/d( y ). ! Coordinate transformation Jacobian real j3 (nx,ny,nz) ! defined as d( zp )/d( z ). ! Coordinate transformation Jacobian real j3inv (nx,ny,nz) ! defined as d( zp )/d( z ). ! EXBC buffer size integer exbcbufsz ! EXBC buffer array real exbcbuf( exbcbufsz ) ! rho-factor: rhobar/rho real rhofct(nx,ny,nz) ! Array containing pre-computed trig real trigs1(3*(nx-1)/2+1) ! function for fftopt=1. ! Array containing pre-computed trig real trigs2(3*(ny-1)/2+1) ! function for fftopt=1. ! Array containing the factors of nx integer ifax1(13) ! for fftopt=1. ! Array containing the factors of ny integer ifax2(13) ! for fftopt=1. ! 2-D work array for fftopt=2 option. real vwork1 (nx+1,ny+1) ! 2-D work array for fftopt=2 option. real vwork2 (ny,nx+1) ! Work array for fftopt=2 option. real wsave1 (3*(ny-1)+15) ! Work array for fftopt=2 option. real wsave2 (3*(nx-1)+15) ! Averaged rhostr at u points (kg/m**3). real rhostru(nx,ny,nz) ! Averaged rhostr at v points (kg/m**3). real rhostrv(nx,ny,nz) ! Averaged rhostr at w points (kg/m**3). real rhostrw(nx,ny,nz) ! !####################################################################### ! ! Temporary WORK ARRAYS: ! !####################################################################### ! ! Pressure gradient term in w-eq. real fp (nx,ny,nz) ! Force in w equation.(kg/(m*s)**2) real fw (nx,ny,nz) ! Temporary work 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) ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i, j, k integer ebc1,wbc1,nbc1,sbc1 real g05,pk,nk,g05wp integer wpprt, itema real tema,temb,w1,w2,nrho integer prgw !Switch for 1st-order or 2nd-order in buoyancy integer buoy2swt ! !####################################################################### ! ! Include files: ! !####################################################################### ! include 'globcst.inc' include 'bndry.inc' include 'exbc.inc' ! Physical constants include 'phycst.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! IF( bsnesq .eq.1 ) THEN ! Switch for pprt/(rhobar*csndsq) term in w-eq. wpprt = 0 ELSE ! Switch for pprt/(rhobar*csndsq) term in w-eq. wpprt = 1 ENDIF prgw = 0 g05 = g*0.5 g05wp = g05*wpprt ! !####################################################################### ! ! Calculate the horizontal velocity divergence using newly updated ! u and v velocity plus half vertical divergence from wcont, and ! store the result in tem3. ! ! Namely, tem3 = difx(u*avgx(j3))+dify(v*avgy(j3))+ ! (1-tacoef)*difz(wcont*avgz(j3)) ! !####################################################################### ! DO 10 k= 2,nz-2 DO 10 j= 1,ny-1 DO 10 i= 1,nx tem1(i,j,k)=u(i,j,k)*rhostru(i,j,k)*mapfct(i,j,5) 10 CONTINUE DO 15 k= 2,nz-2 DO 15 j= 1,ny DO 15 i= 1,nx-1 tem2(i,j,k)=v(i,j,k)*rhostrv(i,j,k)*mapfct(i,j,6) 15 CONTINUE DO 16 k= 2,nz-1 DO 16 j= 1,ny-1 DO 16 i= 1,nx-1 tem3(i,j,k)=wcont(i,j,k)*rhostrw(i,j,k) 16 CONTINUE DO 100 k=2,nz-2 DO 100 j=1,ny-1 DO 100 i=1,nx-1 tem4(i,j,k)=mapfct(i,j,1)*mapfct(i,j,1) *((tem1(i+1,j,k)- & tem1(i,j,k))*dxinv +(tem2(i,j+1,k)-tem2(i,j,k))*dyinv) & +(tem3(i,j,k+1)-tem3(i,j,k))*dzinv*(1.0-tacoef) 100 CONTINUE ! !####################################################################### ! ! Compute the forcing terms in pressure equation ! ! fp = dtsml/j3 *(pforce-csndsq*tem3) ! ! rhostr*g*w is the base state pressure advection term. ! !####################################################################### ! DO 170 k=2,nz-2 DO 170 j=1,ny-1 DO 170 i=1,nx-1 fp(i,j,k)=dtsml1*j3inv(i,j,k) * (pforce(i,j,k)-csndsq(i,j,k)& *tem4(i,j,k)) 170 CONTINUE !####################################################################### ! ! Compute two more terms in fp related to coordinate transformation: ! ! tem3 = (avgx(avgz(ustr)*j1)+avgy(avgz(vstr)*j2))/avgz(j3) ! !####################################################################### ! IF( ternopt.ne.0 ) THEN DO 30 k= 1,nz-1 DO 30 j= 1,ny-1 DO 30 i= 1,nx tem1(i,j,k)=u(i,j,k)*rhostru(i,j,k) 30 CONTINUE DO 35 k= 1,nz-1 DO 35 j= 1,ny DO 35 i= 1,nx-1 tem2(i,j,k)=v(i,j,k)*rhostrv(i,j,k) 35 CONTINUE DO 40 k= 2,nz-1 DO 40 j= 1,ny-1 DO 40 i= 1,nx-1 tem3(i,j,k)=((tem1(i ,j,k)+tem1(i ,j,k-1))*j1(i ,j,k) & +(tem1(i+1,j,k)+tem1(i+1,j,k-1))*j1(i+1,j,k) +(tem2(i & ,j,k)+tem2(i ,j,k-1))*j2(i ,j,k) +(tem2(i,j+1,k)+tem2(i,j+& 1,k-1))*j2(i,j+1,k)) *0.5 * mapfct(i,j,1)/(j3(i,j,k) +j3( & i,j,k-1)) 40 CONTINUE ! !####################################################################### ! ! fp=fp-(dtsml/j3)*csndsq*tacoef* d(tem3)/dz ! !####################################################################### ! DO 400 k=2,nz-2 DO 400 j=1,ny-1 DO 400 i=1,nx-1 fp(i,j,k)=fp(i,j,k) -dtsml1*j3inv(i,j,k)*csndsq(i,j,k) & *tacoef*(tem3(i,j,k+1)-tem3(i,j,k))*dzinv 400 CONTINUE ENDIF ! !####################################################################### ! ! Average rhofct at w points, stored in tem5 ! !####################################################################### ! CALL avgsw(rhofct,nx,ny,nz, 1,nx-1, 1,ny-1,tem5) ! !####################################################################### ! ! Compute the right-hand-side forcing term in tridiagonal linear ! equation. Array w is used to store this forcing term. ! !####################################################################### !####################################################################### ! ! fw=wforce-wpgrad-tacoef*d(fp)/dz ! -g*avgz(rhostr*(pprt+taceof*pfrc)/(gamma*pbar)) ! !####################################################################### ! !1st-order IF ( buoy2nd .eq. 0) THEN !Switch for 1st-order or 2nd-order buoy2swt = 0 !2nd-order ELSE buoy2swt = 1 ENDIF tema = 0.5*buoy2swt*(1.0-cpdcv)/cpdcv temb = 1.0/cpdcv DO 345 k=2,nz-2 DO 345 j=1,ny-1 DO 345 i=1,nx-1 tem1(i,j,k)=rhostr(i,j,k)*((pprt(i,j,k)+tacoef*fp(i,j,k)) & /pbar(i,j,k) +tema*pprt(i,j,k)*pprt(i,j,k) /(pbar(i,j,k)* & pbar(i,j,k)))*temb 345 CONTINUE ! !####################################################################### ! ! When potential temperature equation is solved within small time ! steps, the contribution from ptprt to buoyancy term is calculated ! here. ! !####################################################################### ! IF( ptsmlstp.eq.1 ) THEN DO 352 k=2,nz-2 DO 352 j=1,ny-1 DO 352 i=1,nx-1 tem2(i,j,k)= ptprt(i,j,k)*rhostr(i,j,k)/ptbar(i,j,k) & *(1.0-buoy2swt*ptprt(i,j,k)/ptbar(i,j,k)) 352 CONTINUE ELSE DO 353 k=2,nz-2 DO 353 j=1,ny-1 DO 353 i=1,nx-1 tem2(i,j,k)= 0.0 353 CONTINUE ENDIF tema = tacoef*dzinv DO 350 k=3,nz-2 DO 350 j=1,ny-1 DO 350 i=1,nx-1 fw(i,j,k) = wforce(i,j,k)-tem5(i,j,k)*(wpgrad(i,j,k) & +tema*(fp(i,j,k)-fp(i,j,k-1)) +g05wp*(tem1(i,j,k)+tem1(i,j, & k-1)) -g05 *(tem2(i,j,k)+tem2(i,j,k-1))) 350 CONTINUE ! !####################################################################### ! ! fw=fw*dtsml/rhostr + w ! !####################################################################### ! DO 425 k=3,nz-2 DO 425 j=1,ny-1 DO 425 i=1,nx-1 fw(i,j,k)=fw(i,j,k)*dtsml1/rhostrw(i,j,k) + w(i,j,k) 425 CONTINUE ! !####################################################################### ! ! Prepare the left-hand-side coefficents for tridiagonal ! equation set: ! ! A(k)*w(k-1)+B(k)*w(k)+C(k)*w(k+1)=D(k) for k=3,nz-2 ! ! w(2) and w(nz-1) are used as the boundary conditions. ! ! Due to the lack of work arrays, we are storing the coefficients ! A in rhostru, B in rhostrv, and C in rhostrw. D is stored in fw. ! ! rhostru, rhostrv and rhostrw are re-calculated from rhostr after ! they have been used. ! !####################################################################### ! DO 515 k=2,nz-1 DO 515 j=1,ny-1 DO 515 i=1,nx-1 tem4(i,j,k)=0.5*(rhostr(i,j,k)*j3inv(i,j,k) +rhostr(i,j,k-1)& *j3inv(i,j,k-1)) 515 CONTINUE tema = (dtsml1*tacoef)**2 * wpprt*g*dzinv /cpdcv temb = (dtsml1*tacoef*dzinv)**2 DO 510 k=3,nz-2 DO 510 j=1,ny-1 DO 510 i=1,nx-1 pk = tem5(i,j,k)*tema/(pbar(i,j,k)+pbar(i,j,k-1)) nk = tem5(i,j,k)*temb/rhostrw(i,j,k) tem1(i,j,k)= (-nk+pk)*csndsq(i,j,k-1)*j3inv(i,j,k-1) tem3(i,j,k)=-( nk+pk)*csndsq(i,j,k )*j3inv(i,j,k ) tem2(i,j,k)= 1-tem4(i,j,k)*(tem1(i,j,k)+tem3(i,j,k)) tem1(i,j,k)=tem1(i,j,k)*tem4(i,j,k-1) tem3(i,j,k)=tem3(i,j,k)*tem4(i,j,k+1) 510 CONTINUE ! !####################################################################### ! ! Reset fw on the boundaries using the top and bottom boundary ! conditions of w. ! !####################################################################### ! DO 244 j=1,ny-1 DO 244 i=1,nx-1 w(i,j,nz-1)=0.0 w(i,j,2) =-((u(i ,j,2)+u(i ,j,1))*j1(i,j,2) +(u(i+1,j,2)+u(i+ & 1,j,1))*j1(i+1,j,2) +(v(i,j ,2)+v(i,j ,1))*j2(i,j,2) +(v(i,j+ & 1,2)+v(i,j+1,1))*j2(i,j+1,2)) *0.25*mapfct(i,j,1) 244 CONTINUE DO 600 j=1,ny-1 DO 600 i=1,nx-1 fw(i,j,3) =fw(i,j, 3)-w(i,j, 2)*tem1(i,j, 3) fw(i,j,nz-2)=fw(i,j,nz-2)-w(i,j,nz-1)*tem3(i,j,nz-2) 600 CONTINUE ! !####################################################################### ! ! Call the tridiagonal solver for either a rigid or open upper ! boundary condition for w. ! ! tbc = 4, periodic fft for linearized hydrostatic radiation conditi ! tbc = 6, even fft for linearized hydrostatic radiation condition. ! The references are Klemp and Durran Jas (1983) and Chen (MWR) 1991 ! ! w1 is the coef for w(nz-1) and w2 is the coef. for w(nz-2 in the ! pressure equation. The pressure equation is solved at p(nz-2 ! for w(nz-1). tem4(i,j,nz-2) is the summation of the known terms ! in the pressure equation. The only unknowns are p(nz-2), w(nz-1) ! and w(nz-1) at the new time level. In the tridiagonal solver ! the elimination step is performed and w(nz-2) is given in termsx ! of w(nz-1) and known terms. THe next step is to use the rad. ! condition (in fourier space): ! ! ! p = N * rhobar * w(nz-1) / abs(kx+ky) ! ! The end result is a relation for w(nz-1) given known quantities. ! Trigs1 and trigs2 are the predetermined trig functions used in th ! fft program in tridiag. ! !####################################################################### ! apply linear radiation condition to w at nz-1. IF(tbc.eq.4)THEN tema = rhostr(nx/2,ny/2,nz-2)*csndsq(nx/2,ny/2,nz-2)* j3inv(nx/ & 2,ny/2,nz-2) temb = dtsml1*j3inv(nx/2,ny/2,nz-2) w2=temb*tacoef*(tema*dzinv+g05*prgw*rhostr(nx/2,ny/2,nz-2)) w1=temb*tacoef*(-tema*dzinv+g05*prgw*rhostr(nx/2,ny/2,nz-2)) DO 605 j=1,ny-1 DO 605 i=1,nx-1 tem4(i,j,2)=pprt(i,j,nz-2)+fp(i,j,nz-2) 605 CONTINUE nrho=sqrt(g/ptbar(nx/2,ny/2,nz-2)*(ptbar(nx/2,ny/2,nz-1) & -ptbar(nx/2,ny/2,nz-3))*dzinv*0.5) *rhostr(nx/2,ny/2,nz-2)* & j3inv(nx/2,ny/2,nz-2) END IF ! !####################################################################### ! ! NOTE: tem4(1,1,4) is passed into subroutine tridiag and becomes ! a 2-D array of size (nx+1,ny+1). ! ! rhostrw is used as a work array in tridiag. ! !####################################################################### ! CALL tridiag(nx,ny,nz,tem1,tem2,tem3,fw,rhostrw, & & tem4(1,1,1),tem4(1,1,2),w1,w2,nrho,tem4(1,1,3), & & trigs1,trigs2,ifax1,ifax2,wsave1,wsave2,vwork1,vwork2, & & tem4(1,1,4)) CALL avgsw(rhostr,nx,ny,nz, 1,nx-1, 1,ny-1, rhostrw) ! call TRIDIAG_old(nx,ny,nz,tem1,tem2,tem3,fw, tem4) ! !####################################################################### ! ! On exit of tridiag, the interior solution of w is stored in fw. ! !####################################################################### ! set itema = nz-1 IF(tbc.eq.4)THEN itema = nz-1 ! set itema = nz-2 ELSE itema = nz-2 END IF IF( nestgrd.ne.1 .or. mgrid.eq.1 ) THEN ! For non-nesting case or the ! base-grid of the nested case DO 620 k=3,itema DO 620 j=1,ny-1 DO 620 i=1,nx-1 w(i,j,k) = fw(i,j,k) 620 CONTINUE ! Internal boundary condition 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 lbcw(nx,ny,nz,dtsml1,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, & ebc1,wbc1,nbc1,sbc1) ! External boundary condition ELSE CALL lbcw(nx,ny,nz,dtsml1,w,wcont,wdteb,wdtwb,wdtnb,wdtsb, & 3,3,3,3) ENDIF ! For nested interior grid ELSE DO 640 k=3,itema DO 640 j=2,ny-2 DO 640 i=2,nx-2 w(i,j,k) = fw(i,j,k) 640 CONTINUE DO 700 k=3,nz-2 DO 700 j=1,ny-1 w( 1,j,k)=w( 1,j,k)+dtsml1 * wdtwb(j,k) w(nx-1,j,k)=w(nx-1,j,k)+dtsml1 * wdteb(j,k) 700 CONTINUE DO 710 k=3,nz-2 DO 710 i=2,nx-2 w(i, 1,k)=w(i, 1,k)+dtsml1 * wdtsb(i,k) w(i,ny-1,k)=w(i,ny-1,k)+dtsml1 * wdtnb(i,k) 710 CONTINUE ENDIF ! !####################################################################### ! ! Calculate wcont at time tfuture, including the boundary ! points. Wcont at the lateral boundaries is calculated ! from boundary u, v and w. Wcont at the top and bottom ! depends on the boundary condition option. ! !####################################################################### ! CALL wcontra(nx,ny,nz,u,v,w,mapfct,j1,j2,j3, & & rhostr,rhostru,rhostrv,rhostrw,wcont,tem1,tem2) ! !####################################################################### ! ! Set the top and bottom boundary conditions for w based on u, v and ! wcont at the top and bottom boundaries. ! !####################################################################### ! CALL vbcw(nx,ny,nz,w,wcont,tbc,bbc,u,v, & & rhostr,rhostru,rhostrv,rhostrw, & & j1,j2,j3) ! !####################################################################### ! ! Calculate the new pressure ! ! pprt(new)=pprt(old)+fp-dtsml/j3*tacoef*csndsq*difz(rhobar*w)) ! !####################################################################### ! tema = tacoef*dtsml1*dzinv DO 795 k=2,nz-1 DO 795 j=1,ny-1 DO 795 i=1,nx-1 tem1(i,j,k)= 0.5*(rhostr(i,j,k )*j3inv(i,j,k ) +rhostr(i,j, & k-1)*j3inv(i,j,k-1)) 795 CONTINUE DO 790 k=2,nz-2 DO 790 j=1,ny-1 DO 790 i=1,nx-1 fp(i,j,k)=fp(i,j,k)-tema*j3inv(i,j,k)*csndsq(i,j,k)* & (tem1(i,j,k+1)*w(i,j,k+1)-tem1(i,j,k)*w(i,j,k)) 790 CONTINUE IF( nestgrd.eq.1 .and. mgrid.ne.1 ) THEN ! For nested interior grid DO 800 k=2,nz-2 DO 800 j=2,ny-2 DO 800 i=2,nx-2 pprt(i,j,k)=pprt(i,j,k)+fp(i,j,k) 800 CONTINUE DO 900 k=2,nz-2 DO 900 j=1,ny-1 pprt( 1,j,k)=pprt( 1,j,k)+dtsml1 * pdtwb(j,k) pprt(nx-1,j,k)=pprt(nx-1,j,k)+dtsml1 * pdteb(j,k) 900 CONTINUE DO 910 k=2,nz-2 DO 910 i=2,nx-2 pprt(i, 1,k)=pprt(i, 1,k)+dtsml1 * pdtsb(i,k) pprt(i,ny-1,k)=pprt(i,ny-1,k)+dtsml1 * pdtnb(i,k) 910 CONTINUE ! ! Call the pprt bottom boundary condition subroutine to ! compute the hydrostatic pprt at k=1. ! CALL pprtbbc(nx,ny,nz,g05,buoy2swt,rhostr,pprt,ptprt, pbar, & ptbar,phydro, tem1) ! new pprt at k=1. CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, & tem1(1,1,1), 0,0,0,0 ,tbc,bbc) ! For non-nesting case or the ELSE ! base-grid of the nested case DO 810 k=2,nz-2 DO 810 j=1,ny-1 DO 810 i=1,nx-1 pprt(i,j,k)=pprt(i,j,k)+fp(i,j,k) 810 CONTINUE ! ! Call the pprt bottom boundary condition subroutine to ! compute the hydrostatic pprt at k=1. ! CALL pprtbbc(nx,ny,nz,g05,buoy2swt,rhostr,pprt,ptprt, pbar, & ptbar,phydro, tem1) ! new pprt at k=1. ! !####################################################################### ! ! Apply the boundary conditions for the pressure equation. ! ! For the open boundary case, the boundary value of pprt is ! predicted by the pressure equation. ! !####################################################################### ! ! Internal boundary conditions 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 bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, & tem1(1,1,1), ebc1,wbc1,nbc1,sbc1,tbc,bbc) ! External boundary condition ELSE CALL bcp(nx,ny,nz, dtsml1, pprt, pdteb, pdtwb, pdtnb, pdtsb, & tem1(1,1,1), 0,0,0,0,tbc,bbc) CALL exbcp(nx,ny,nz, curtsml, pprt, exbcbuf(npr0exb),exbcbuf( & nprdtexb)) ENDIF ENDIF RETURN END ! ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE TRIDIAG ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! SUBROUTINE TRIDIAG(nx,ny,nz,a,b,c,d, tem2, tem1, temxy2,w1,w2, & & nrho,work,trigs1,trigs2,ifax1,ifax2, & & wsave1,wsave2,vwork1,vwork2, & & temxy1) ! !####################################################################### ! ! PURPOSE: ! ! Solve a tridiagonal linear equation. ! ! The tridiagonal equation set to be solved is ! ! a(k)*w(k-1)+b(k)*w(k)+c(k)*w(k+1)=d(k) for k=3,nz-2 ! ! given w(2) and w(nz-1) as the boundary conditions. ! ! The solution w(k) is stored directly into d(k). ! ! Reference: Numerical Recipes: FORTRAN Version, 1989. Page 40. ! !####################################################################### ! ! AUTHOR: Ming Xue ! 9/9/94 ! ! MODIFICATION HISTORY: ! ! 10/31/95 (D. Weber) ! Added linear hydrostatic upper radiation condition for w-p. ! References are Klemp and Durran (JAS, 1983) and Chen (1991). ! Includes the use of trigs1,trigs2,ifax1,ifax2. ! ! 07/22/97 (D. Weber) ! Added linear hydrostatic upper radiation condition for w-p using ! an even fft (fftopt=2). ! !####################################################################### ! ! INPUT : ! ! nx Number of grid points in the x-direction ! ny Number of grid points in the y-direction ! nz Number of grid points in the vertical ! ! a lhs coefficient of tridigonal eqations ! b lhs coefficient of tridigonal eqations ! c lhs coefficient of tridigonal eqations ! d rhs coefficient of tridigonal eqations ! temxy2 Pforce array at nz-2. ! trigs1 Array containing pre-computed trig function for ! fftopt=1. ! trigs2 Array containing pre-computed trig function for ! fftopt=1. ! ifax1 Array containing the factors of nx for fftopt=1. ! ifax2 Array containing the factors of ny for fftopt=1. ! ! vwork1 2-D work array for fftopt=2 option. ! vwork2 2-D work array for fftopt=2 option. ! wsave1 Work array for fftopt=2. ! wsave2 Work array for fftopt=2. ! ! OUTPUT: ! ! d The solution w. ! ! WORK ARRAYS: ! ! temxy1 Work array at nz-2. ! tem1 Work array ! tem2 Work array ! work Work array for fft. ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations: ! !####################################################################### ! ! Force explicit declarations implicit none ! ! Include files: ! include 'bndry.inc' include 'globcst.inc' ! Number of grid points in 3 directions integer nx, ny, nz ! lhs coefficient of tridigonal eqations real a (nx,ny,nz) ! lhs coefficient of tridigonal eqations real b (nx,ny,nz) ! lhs coefficient of tridigonal eqations real c (nx,ny,nz) ! rhs coefficient of tridigonal eqations real d (nx,ny,nz) ! The contains solution w on exit. ! Pforce array at nz-2. real temxy2(nx,ny) ! Array containing pre-computed trig real trigs1(3*(nx-1)/2+1) ! function for fftopt=1. ! Array containing pre-computed trig real trigs2(3*(ny-1)/2+1) ! function for fftopt=1. ! Array containing the factors of nx for integer ifax1(13) ! fftopt=1. ! Array containing the factors of ny for integer ifax2(13) ! fftopt=1. ! 2-D work array for fftopt=2. real vwork1 (nx+1,ny+1) ! 2-D work array for fftopt=2. real vwork2 (ny,nx+1) ! Work array for fftopt=2. real wsave1 (3*(ny-1)+15) ! Work array for fftopt=2. real wsave2 (3*(nx-1)+15) ! Work array at nz-2. real temxy1(nx+1,ny+1) ! 2-D work array. real tem1 (nx,ny) ! 3-D work array. real tem2 (nx,ny,nz) ! 2-D work array for fft code. real work (ny,nx) ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i, j, k real nrho real w1,w2,wc ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! DO 100 j=1,ny-1 DO 100 i=1,nx-1 d(i,j,3)=d(i,j,3)/b(i,j,3) tem1(i,j)=b(i,j,3) 100 CONTINUE DO 200 k=4,nz-2 DO 200 j=1,ny-1 DO 200 i=1,nx-1 tem2(i,j,k) = c(i,j,k-1)/tem1(i,j) tem1(i,j)=b(i,j,k)-a(i,j,k)*tem2(i,j,k) d(i,j,k)=(d(i,j,k)-a(i,j,k)*d(i,j,k-1))/tem1(i,j) 200 CONTINUE ! apply upper radiation boundary condition. IF(tbc.eq.4)THEN ! ! Computing the new c(nz-2) (which is stored in nz-1).... ! DO 210 j=1,ny-1 DO 210 i=1,nx-1 tem2(i,j,nz-1) = c(i,j,nz-2)/tem1(i,j) 210 CONTINUE ! coef. for w(nz-1) in tri. elim. pha wc = tem2(nx/2,ny/2,nz-1) ! ! Combining the pforce array and d(nz-2) from the elimination phase. ! Note: these arrays are a function of x and/or y. All others varia ! are determined from base state which are not a function ! of x,y. (note: zflat must be at or below the height of scalar nz- ! DO 220 j=1,ny-1 DO 220 i=1,nx-1 temxy1(i,j) = temxy2(i,j) + w2*d(i,j,nz-2) 220 CONTINUE !####################################################################### ! ! Call the upper radiation subroutine to update w at the top (nz-2) ! !####################################################################### ! call the even vfftpack fft... IF(fftopt.eq.2)THEN CALL uprad3(nx,ny,w1,w2,wc,nrho,wsave1,wsave2,vwork1,vwork2, & temxy1,work) ! call the periodic fft.... ELSEIF(fftopt.eq.1)THEN CALL uprad1(nx,ny,w1,w2,wc,nrho,ifax1,ifax2,trigs1,trigs2, & temxy1,work) ! end of fftopt.... END IF DO 230 j=1,ny-1 DO 230 i=1,nx-1 d(i,j,nz-1) = temxy1(i,j) 230 CONTINUE !####################################################################### ! ! Perform the back substitution phase of the tridiagonal solver. ! !####################################################################### DO 300 k=nz-2,3,-1 DO 300 j=1,ny-1 DO 300 i=1,nx-1 d(i,j,k)=d(i,j,k)-tem2(i,j,k+1)*d(i,j,k+1) 300 CONTINUE ! perform backsubtitution for other bc's. ELSEIF(tbc.ne.4)THEN DO 310 k=nz-3,3,-1 DO 310 j=1,ny-1 DO 310 i=1,nx-1 d(i,j,k)=d(i,j,k)-tem2(i,j,k+1)*d(i,j,k+1) 310 CONTINUE ! end of tbc loop... ENDIF RETURN END ! ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE WCONTRA ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for the Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! SUBROUTINE WCONTRA(nx,ny,nz,u,v,w,mapfct,j1,j2,j3, & & rhostr,rhostru,rhostrv,rhostrw,wcont,ustr,vstr) ! !####################################################################### ! ! PURPOSE: ! ! Calculate wcont, the contravariant vertical velocity (m/s) ! !####################################################################### ! ! AUTHOR: Ming Xue & Hao Jin ! 1/4/1993. ! ! Modification history: ! 8/29/94 (A. Shapiro) ! Bug fix. Call to vbcwcont moved outside IF block. ! ! 9/9/94 (M. Xue) ! Optimized. ! ! 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. ! !####################################################################### ! ! 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 ! ! 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 Cartesian velocity ! at all time levels (m/s) ! ! mapfct Map factors at scalar, u and v points ! ! j1 Coordinate transform Jacobian -d(zp)/dx ! j2 Coordinate transform Jacobian -d(zp)/dy ! j3 Coordinate transform Jacobian d(zp)/dz ! ! rhostr j3 times base state density rhobar(kg/m**3). ! rhostru Average rhostr at u points (kg/m**3). ! rhostrv Average rhostr at v points (kg/m**3). ! rhostrw Average rhostr at w points (kg/m**3). ! ! OUTPUT: ! ! wcont Vertical component of contravariant velocity in ! computational coordinates (m/s) ! ! WORK ARRAYS: ! ! ustr Work array ! vstr Work array ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations. ! !####################################################################### ! implicit none ! The number of grid points in 3 integer nx,ny,nz ! directions ! Total u-velocity (m/s) real u (nx,ny,nz) ! Total v-velocity (m/s) real v (nx,ny,nz) ! Total w-velocity (m/s) real w (nx,ny,nz) ! Map factors at scalar, u and v points real mapfct(nx,ny,6) ! Coordinate transform Jacobian real j1 (nx,ny,nz) ! defined as - d( zp )/d( x ). ! Coordinate transform Jacobian real j2 (nx,ny,nz) ! defined as - d( zp )/d( y ). ! Coordinate transform Jacobian real j3 (nx,ny,nz) ! defined as d( zp )/d( z ). ! j3 times base state density rhobar real rhostr(nx,ny,nz) ! (kg/m**3). ! Average rhostr at u points (kg/m**3). real rhostru(nx,ny,nz) ! Average rhostr at v points (kg/m**3). real rhostrv(nx,ny,nz) ! Average rhostr at w points (kg/m**3). real rhostrw(nx,ny,nz) ! Vertical velocity in computational real wcont (nx,ny,nz) ! coordinates (m/s) ! temporary work array real ustr (nx,ny,nz) ! temporary work array real vstr (nx,ny,nz) !####################################################################### ! ! Include files: ! !####################################################################### ! include 'globcst.inc' ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i,j,k ! !####################################################################### !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! No coord. transformation case. IF( crdtrns .eq. 0 ) THEN DO 10 k= 2,nz-1 DO 10 j= 1,ny-1 DO 10 i= 1,nx-1 wcont(i,j,k)=w(i,j,k) 10 CONTINUE ELSEIF( ternopt.eq.0) THEN DO 20 k= 2,nz-1 DO 20 j= 1,ny-1 DO 20 i= 1,nx-1 wcont(i,j,k)=w(i,j,k)*2.0/(j3(i,j,k)+j3(i,j,k-1)) 20 CONTINUE ELSE DO 30 k= 1,nz-1 DO 30 j= 1,ny-1 DO 30 i= 1,nx ustr(i,j,k)=u(i,j,k)*rhostru(i,j,k) 30 CONTINUE DO 35 k= 1,nz-1 DO 35 j= 1,ny DO 35 i= 1,nx-1 vstr(i,j,k)=v(i,j,k)*rhostrv(i,j,k) 35 CONTINUE DO 40 k= 2,nz-1 DO 40 j= 1,ny-1 DO 40 i= 1,nx-1 wcont(i,j,k)= ( ((ustr(i ,j,k)+ustr(i ,j,k-1))*j1(i & ,j,k) +(ustr(i+1,j,k)+ustr(i+1,j,k-1))*j1(i+1,j,k) & +(vstr(i ,j,k)+vstr(i ,j,k-1))*j2(i ,j,k) +(vstr(i,j+1,k)+& vstr(i,j+1,k-1))*j2(i,j+1,k)) * 0.25 * mapfct(i,j,1) & / rhostrw(i,j,k) + w(i,j,k) ) * 2.0/(j3(i,j,k)+j3(i,j,k-1)& ) 40 CONTINUE ENDIF CALL vbcwcont(nx,ny,nz,wcont) RETURN END ! ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE SOLVPT_LRG ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! SUBROUTINE SOLVPT_LRG(nx,ny,nz, exbcbufsz, dtbig1, & & ptprt, ptdteb,ptdtwb,ptdtnb,ptdtsb, & & rhostr,ptforce, exbcbuf) ! !####################################################################### ! ! PURPOSE: ! ! Coordinate the time integration of the potential temperature ! equation in large time steps. ! !####################################################################### ! ! AUTHOR: Ming Xue ! 10/10/91 ! ! MODIFICATION HISTORY: ! ! 9/17/94 (M. Xue) ! Rewritten for small time step integration of ptprt equation. ! ! 11/6/1995 (M. Xue) ! Added option for fourth order vertical advection for ptbar. ! ! 4/17/96 (M. Xue) ! Removed the block for 4th order advection of ptbar. ! ! 4/24/1997 (M. Xue) ! Rewrote as SOLVPT_LRG based on SOLVPT. ! !####################################################################### ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! ! dtbig1 The big time step size for this call (s) ! ! ptprt Perturbation potential temperature at times tpast and ! tpresent (K) ! ! ptdteb Time tendency of the ptprt field at the east boundary ! ptdtwb Time tendency of the ptprt field at the west boundary ! ptdtnb Time tendency of the ptprt field at the north boundary ! ptdtsb Time tendency of the ptprt field at the south boundary ! ! rhostr Base state density rhobar times j3 (kg/m**3) ! ! ptforce Gravity wave inactive forcing terms in pt-eq. ! (K*kg/(m**3*s)) ! ! OUTPUT: ! ! ptprt Perturbation potential temperature at time tfuture (K) ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations: ! !####################################################################### ! ! Force explicit declarations implicit none ! The no. of t-levels of integer nt ! t-dependent arrays. ! Index of time level for the integer tpast ! past time. ! Index of time level for the integer tpresent ! present time. ! Index of time level for the integer tfuture ! future time. parameter (nt=3, tpast=1, tpresent=2, tfuture=3) ! Number of grid points in 3 directions integer nx, ny, nz ! The big time step size for this call real dtbig1 ! (s) ! Perturbation potential temperature (K) real ptprt (nx,ny,nz,nt) ! Time tendency of ptprt at east real ptdteb(ny,nz) ! boundary ! Time tendency of ptprt at west real ptdtwb(ny,nz) ! boundary ! Time tendency of ptprt at north real ptdtnb(nx,nz) ! boundary ! Time tendency of ptprt at south real ptdtsb(nx,nz) ! boundary ! Base state density rhobar times j3. real rhostr(nx,ny,nz) ! Gravity wave inactive forcing terms real ptforce(nx,ny,nz) ! in potential temperature eq. ! (K*kg/(m**3*s)) ! EXBC buffer size integer exbcbufsz ! EXBC buffer array real exbcbuf( exbcbufsz ) ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i, j, k, tstrtlvl integer ebc1,wbc1,nbc1,sbc1 real deltat ! !####################################################################### ! ! Include files: ! !####################################################################### ! include 'globcst.inc' include 'bndry.inc' include 'exbc.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! !####################################################################### ! ! Integrate forward by one timestep the potential temperature ! equation. When PT-eq. is integrated inside small time steps, ! it is stepped forward by a small time step, otherwise, it is ! stepped forward by a large time step using leapfrog scheme ! (i.e. 2*dtbig1 from ptprt at tpast). ! !####################################################################### ! ! Forward-based FCT scheme IF(sadvopt.eq.4) THEN deltat = dtbig1 tstrtlvl = tpresent ELSE deltat = dtbig1*2 tstrtlvl = tpast ENDIF IF( nestgrd.eq.1 .and. mgrid.ne.1 ) THEN DO 100 k=2,nz-2 DO 100 j=2,ny-2 DO 100 i=2,nx-2 ptprt(i,j,k,tfuture)=ptprt(i,j,k,tstrtlvl) +deltat* & ptforce(i,j,k)/rhostr(i,j,k) 100 CONTINUE ! !####################################################################### ! ! Integrate PT equation and the boundary conditions for a nested gri ! !####################################################################### ! DO 110 k=2,nz-2 DO 110 i=1,nx-1 ptprt(i, 1,k,tfuture)= 2*ptprt(i, 1,k,tpresent)-ptprt(i, & 1,k,tpast) ptprt(i,ny-1,k,tfuture)= 2*ptprt(i,ny-1,k,tpresent)-ptprt(i,& ny-1,k,tpast) 110 CONTINUE DO 120 k=2,nz-2 DO 120 j=2,ny-2 ptprt( 1,j,k,tfuture)= 2*ptprt( 1,j,k,tpresent)-ptprt( & 1,j,k,tpast) ptprt(nx-1,j,k,tfuture)= 2*ptprt(nx-1,j,k,tpresent)-ptprt( & nx-1,j,k,tpast) 120 CONTINUE CALL bcsclr(nx,ny,nz, deltat*0.5 , ptprt(1,1,1,tstrtlvl),ptprt( & 1,1,1,tpresent), ptprt(1,1,1,tfuture),ptdteb,ptdtwb,ptdtnb, & ptdtsb, 0,0,0,0 ,tbc,bbc) ELSE ! !####################################################################### ! ! Integrate PT equation and the boundary conditions for a base grid. ! !####################################################################### ! DO 200 k=2,nz-2 DO 200 j=1,ny-1 DO 200 i=1,nx-1 ptprt(i,j,k,tfuture)=ptprt(i,j,k,tstrtlvl) +deltat* & ptforce(i,j,k)/rhostr(i,j,k) 200 CONTINUE ! Internal boundary conditions 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 bcsclr(nx,ny,nz, deltat*0.5 , ptprt(1,1,1,tstrtlvl), & ptprt(1,1,1,tpresent), ptprt(1,1,1,tfuture),ptdteb,ptdtwb, & ptdtnb,ptdtsb, ebc1,wbc1,nbc1,sbc1,tbc,bbc) ! External boundary condition ELSE CALL bcsclr(nx,ny,nz, deltat*0.5 , ptprt(1,1,1,tstrtlvl), & ptprt(1,1,1,tpresent), ptprt(1,1,1,tfuture),ptdteb,ptdtwb, & ptdtnb,ptdtsb, 0,0,0,0,tbc,bbc) CALL exbcpt(nx,ny,nz, curtim+dtbig, ptprt(1,1,1,tfuture), & exbcbuf(npt0exb),exbcbuf(nptdtexb)) ENDIF ENDIF RETURN END ! ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE SOLVPT_SML ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! SUBROUTINE SOLVPT_SML(nx,ny,nz, exbcbufsz, dtbig1,dtsml1,curtsml, & & ptprt,w, ptdteb,ptdtwb,ptdtnb,ptdtsb, & & ptbar,rhostr,rhostrw,j3,ptforce, & & exbcbuf, & & ptadv,tem1) ! !####################################################################### ! ! PURPOSE: ! ! Coordinate the time integration of the potential temperature ! equation in small time steps. ! !####################################################################### ! ! AUTHOR: Ming Xue ! 10/10/91 ! ! MODIFICATION HISTORY: ! ! 9/17/94 (M. Xue) ! Rewritten for small time step integration of ptprt equation. ! ! 11/6/1995 (M. Xue) ! Added option for fourth order vertical advection for ptbar. ! ! 4/17/96 (M. Xue) ! Removed the block for 4th order advection of ptbar. ! ! 4/24/1997 (M. Xue) ! Rewrote as SOLVPT_SML based on SOLVPT. ! !####################################################################### ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! ! dtbig1 The big time step size for this call (s) ! dtsml1 The big time step size for this call (s) ! curtsml Current time during small time step integration. ! ! ptprt Perturbation potential temperature at times tpast and ! tpresent (K) ! w Vertical component of Cartesian velocity ! ! ptdteb Time tendency of the ptprt field at the east boundary ! ptdtwb Time tendency of the ptprt field at the west boundary ! ptdtnb Time tendency of the ptprt field at the north boundary ! ptdtsb Time tendency of the ptprt field at the south boundary ! ! ptbar Base state potential temperature (K) ! rhostr Base state density rhobar times j3 (kg/m**3) ! rhostrw rhostr averaged to w-point. ! ! j3 Coordinate transformation Jacobian d(zp)/dz ! ! ptforce Gravity wave inactive forcing terms in pt-eq. ! (K*kg/(m**3*s)) ! ! OUTPUT: ! ! ptprt Perturbation potential temperature at time tfuture (K) ! ! WORK ARRAYS: ! ! ptadv Advection of base state potential temperature ! tem1 Temporary work array. ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations: ! !####################################################################### ! ! Force explicit declarations implicit none ! The no. of t-levels of integer nt ! t-dependent arrays. ! Index of time level for the integer tpast ! past time. ! Index of time level for the integer tpresent ! present time. ! Index of time level for the integer tfuture ! future time. parameter (nt=3, tpast=1, tpresent=2, tfuture=3) ! Number of grid points in 3 directions integer nx, ny, nz ! The big time step size for this call real dtbig1 ! (s) ! The big time step size for this call real dtsml1 ! (s) ! Current time during small time step real curtsml ! integration. ! Perturbation potential temperature (K) real ptprt (nx,ny,nz,nt) ! Total w-velocity (m/s) real w (nx,ny,nz,nt) ! Time tendency of ptprt at east real ptdteb(ny,nz) ! boundary ! Time tendency of ptprt at west real ptdtwb(ny,nz) ! boundary ! Time tendency of ptprt at north real ptdtnb(nx,nz) ! boundary ! Time tendency of ptprt at south real ptdtsb(nx,nz) ! boundary ! Base state potential temperature (K) real ptbar (nx,ny,nz) ! Base state density rhobar times j3. real rhostr(nx,ny,nz) ! rhostr averaged to w-point. real rhostrw(nx,ny,nz) ! Coordinate transformation Jacobian real j3 (nx,ny,nz) ! defined as d( zp )/d( z ). ! Gravity wave inactive forcing terms real ptforce(nx,ny,nz) ! in potential temperature eq. ! (K*kg/(m**3*s)) ! EXBC buffer size integer exbcbufsz ! EXBC buffer array real exbcbuf( exbcbufsz ) ! Temporary array to store base state real ptadv (nx,ny,nz) ! potential temperature advection. ! Temporary work array real tem1 (nx,ny,nz) ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i, j, k, onvf integer ebc1,wbc1,nbc1,sbc1 real deltat, dz05 , time ! !####################################################################### ! ! Include files: ! !####################################################################### ! include 'globcst.inc' include 'bndry.inc' include 'exbc.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! !####################################################################### ! ! Integrate forward by one timestep the potential temperature ! equation. When PT-eq. is integrated inside small time steps, ! it is stepped forward by a small time step, otherwise, it is ! stepped forward by a large time step using leapfrog scheme ! (i.e. 2*dtbig1 from ptprt at tpast). ! !####################################################################### ! deltat = dtsml1 time = curtsml ! !####################################################################### ! ! Base state potential temperature advection. This term is added to ! the array ptadv to yield the total potential temperature advection ! ! ptbar is assumed to be independent of physical x and y, therefore ! d(ptbar)/dx and d(ptbar)/dy for constant z are zero, the base stat ! advection is -w*d(ptbar)/dzp = -w/j3*d(ptbar)/dz. ! ! This term is responsible for internal gravity waves, when potentia ! temperature equation is solved inside the small time steps, this ! term is evaluated inside the small time steps. ! !####################################################################### ! dz05 = dz*0.5 DO 360 k=2,nz-1 DO 360 j=1,ny-1 DO 360 i=1,nx-1 tem1(i,j,k)=rhostrw(i,j,k)*w(i,j,k,tfuture) *(ptbar(i,j,k)- & ptbar(i,j,k-1)) /((j3(i,j,k)+j3(i,j,k-1))*dz05) 360 CONTINUE onvf = 0 CALL avgz(tem1 , onvf, & & nx,ny,nz, 1,nx-1, 1,ny-1, 2,nz-2, ptadv) IF( nestgrd.eq.1 .and. mgrid.ne.1 ) THEN DO 100 k=2,nz-2 DO 100 j=2,ny-2 DO 100 i=2,nx-2 ptprt(i,j,k,tfuture)=ptprt(i,j,k,tfuture) +deltat*( & ptforce(i,j,k)-ptadv(i,j,k))/rhostr(i,j,k) 100 CONTINUE ! !####################################################################### ! ! Integrate PT equation and the boundary conditions for a nested gri ! !####################################################################### ! DO 110 k=2,nz-2 DO 110 j=1,ny-1 ptprt( 1,j,k,tfuture)=ptprt( 1,j,k,tfuture) +deltat* ptdtwb(& j,k) ptprt(nx-1,j,k,tfuture)=ptprt(nx-1,j,k,tfuture) +deltat* & ptdteb(j,k) 110 CONTINUE DO 120 k=2,nz-2 DO 120 i=2,nx-2 ptprt(i, 1,k,tfuture)=ptprt(i, 1,k,tfuture) +deltat* ptdtsb(& i,k) ptprt(i,ny-1,k,tfuture)=ptprt(i,ny-1,k,tfuture) +deltat* & ptdtnb(i,k) 120 CONTINUE CALL bcsclr(nx,ny,nz, deltat*0.5 , ptprt(1,1,1,tfuture),ptprt(1,& 1,1,tpresent), ptprt(1,1,1,tfuture),ptdteb,ptdtwb,ptdtnb,ptdtsb,& 0,0,0,0 ,tbc,bbc) ELSE ! !####################################################################### ! ! Integrate PT equation and the boundary conditions for a base grid. ! !####################################################################### ! DO 200 k=2,nz-2 DO 200 j=1,ny-1 DO 200 i=1,nx-1 ptprt(i,j,k,tfuture)=ptprt(i,j,k,tfuture) +deltat*( & ptforce(i,j,k)-ptadv(i,j,k))/rhostr(i,j,k) 200 CONTINUE ! Internal boundary conditions 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 bcsclr(nx,ny,nz, deltat*0.5 , ptprt(1,1,1,tfuture),ptprt(& 1,1,1,tpresent), ptprt(1,1,1,tfuture),ptdteb,ptdtwb,ptdtnb, & ptdtsb, ebc1,wbc1,nbc1,sbc1,tbc,bbc) ! External boundary condition ELSE CALL bcsclr(nx,ny,nz, deltat*0.5 , ptprt(1,1,1,tfuture),ptprt(& 1,1,1,tpresent), ptprt(1,1,1,tfuture),ptdteb,ptdtwb,ptdtnb, & ptdtsb, 0,0,0,0,tbc,bbc) CALL exbcpt(nx,ny,nz, time , ptprt(1,1,1,tfuture), exbcbuf( & npt0exb),exbcbuf(nptdtexb)) ENDIF ENDIF RETURN END ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE SOLVQV ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! ! SUBROUTINE SOLVQV(nx,ny,nz, exbcbufsz, dtbig1, & & qv,u,v,wcont, qvdteb,qvdtwb,qvdtnb,qvdtsb, & & rhostr,qvbar,kmh,kmv,rprntl,qvsflx,pbldpth, & & x,y,z,zp, mapfct, j1,j2,j3,j3inv,qvcumsrc, & & usflx,vsflx,ptsflx,ptsfc,qvsfc,ptbar,ptprt, & & 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 equation for water vapor ! specific humidity qv, and the setting of boundary conditions. ! !####################################################################### ! ! AUTHOR: Ming Xue ! 10/10/91 ! ! MODIFICATION HISTORY: ! ! 5/20/92 (M. Xue) ! Added full documentation. ! ! 2/10/93 (K. Droegemeier) ! Cleaned up documentation. ! ! 4/10/93 (M. Xue & Hao Jin) ! Add the terrain. ! ! 8/22/95 (M. Xue) ! Added ptcumsrc term to the right hand side of qv equation. ! It was omitted before. ! ! 08/30/95 (M. Xue and Y. Liu) ! Changed EXBC for water variables to zero-gradient when the ! variables are missing in the boundary file. ! ! 08/30/95 (Yuhe Liu) ! Fixed a bug which double computes the qv advection and mixing. ! ! 1/25/96 (Donghai Wang & Yuhe Liu) ! Added the map projection factor to ARPS governing equations. ! ! 4/1/96 (Donghai Wang, X. Song and M. Xue) ! Added the implicit treatment for the vertical mixing. ! ! 3/24/1997 (M. Xue) ! Code to handle the case of forward time integration added. ! ! 7/10/1997 (Fanyou Kong -- CMRP) ! Added the positive definite advection scheme option (sadvopt = 5) ! ! 7/17/1998 (M. Xue) ! Changed call to ADVQFCT. ! !####################################################################### ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! ! dtbig1 The big time step size (s) ! ! qv Water vapor specific humidity at tpast and tpresent ! (kg/kg) ! ! u x component of velocity at all time levels (m/s) ! v y component of velocity at all time levels (m/s) ! wcont Contravariant vertical velocity (m/s) ! ! qvdteb Time tendency of the qv field at the east boundary ! qvdtwb Time tendency of the qv field at the west boundary ! qvdtnb Time tendency of the qv field at the north boundary ! qvdtsb Time tendency of the qv field at the south boundary ! ! rhostr Base state density rhobar times j3 (kg/m**3) ! qvbar Base state water vapor specific humidity (kg/kg) ! ! kmh Horizontal turb. mixing coef. for momentum ( m**2/s ) ! kmv Vertical turb. mixing coef. for momentum ( m**2/s ) ! rprntl Reciprocal of Prandtl number ! ! qvsflx Surface flux of moisture (kg/(m**2*s)). ! pbldpth Planetary boundary layer depth (m) ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp Vertical coordinate of grid points in physical space(m) ! ! mapfct Map factors at scalar points ! ! j1 Coordinate transformation Jacobian -d(zp)/dx ! j2 Coordinate transformation Jacobian -d(zp)/dy ! j3 Coordinate transformation Jacobian d(zp)/dz ! qvcumsrc Source term in qv-equation due to cumulus parameterizat ! ! OUTPUT: ! ! qv Water vapor specific humidity at time tfuture (kg/kg) ! ! WORK ARRAYS: ! ! qadv Advection term of water vapor eq. (kg/(m**3*s)). ! A local array. ! qmix Total mixing in water vapor eq. (kg/(m**3*s)). ! A local array. ! tem1 Temporary work array. ! tem2 Temporary work array. ! tem3 Temporary work array. ! tem4 Temporary work array. ! tem5 Temporary work array. ! tem6 Temporary work array. ! tem6 Temporary work array. ! tem7 Temporary work array. ! tem8 Temporary work array. ! tem9 Temporary work array. ! tem10 Temporary work array. ! tem11 Temporary work array. ! ! tem1_0 Temporary work array. ! tem2_0 Temporary work array. ! tem3_0 Temporary work array. ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations: ! !####################################################################### ! ! Force explicit declarations implicit none ! The no. of t-levels of t-dependent integer nt ! arrays. ! Index of time level for the past time. integer tpast ! Index of time level for the present integer tpresent ! time. ! Index of time level for the future integer tfuture ! time. parameter (nt=3, tpast=1, tpresent=2, tfuture=3) ! Number of grid points in 3 directions integer nx, ny, nz ! Local big time step real dtbig1 ! Water vapor specific humidity (kg/kg) real qv (nx,ny,nz,nt) ! Total u-velocity (m/s) real u (nx,ny,nz,nt) ! Total v-velocity (m/s) real v (nx,ny,nz,nt) ! Contravariant vertical velocity (m/s) real wcont (nx,ny,nz) ! Time tendency of qv at east boundary real qvdteb(ny,nz) ! Time tendency of qv at west boundary real qvdtwb(ny,nz) ! Time tendency of qv at north boundary real qvdtnb(nx,nz) ! Time tendency of qv at south boundary real qvdtsb(nx,nz) ! Base state density rhobar times j3. real rhostr(nx,ny,nz) ! Base state water vapor specific real qvbar (nx,ny,nz) ! humidity (kg/kg) ! Horizontal turb. mixing coef. for real kmh (nx,ny,nz) ! momentum. ( m**2/s ) ! Vertical turb. mixing coef. for real kmv (nx,ny,nz) ! momentum. ( m**2/s ) ! Reciprocal of Prandtl number real rprntl(nx,ny,nz) ! Surface flux of u-momentum real usflx(nx,ny) ! Surface flux of v-momentum real vsflx(nx,ny) ! Surface heat flux (K*kg/(m*s**2)) real ptsflx(nx,ny) ! Surface flux of moisture (kg/(m**2*s)) real qvsflx(nx,ny) ! Planetary boundary layer depth (m) real pbldpth(nx,ny,nt) ! x-coord. of the physical and compu- real x (nx) ! tational grid. Defined at u-point. ! y-coord. of the physical and compu- real y (ny) ! tational grid. Defined at v-point. ! z-coord. of the computational grid. real z (nz) ! Defined at w-point on the staggered ! grid. ! Physical height coordinate defined at real zp (nx,ny,nz) ! w-point of the staggered grid. ! Map factors at scalar points real mapfct(nx,ny) ! Coordinate transformation Jacobian real j1 (nx,ny,nz) ! defined as - d( zp )/d( x ). ! Coordinate transformation Jacobian real j2 (nx,ny,nz) ! defined as - d( zp )/d( y ). ! Coordinate transformation Jacobian real j3 (nx,ny,nz) ! defined as d( zp )/d( z ). ! Coordinate transformation Jacobian real j3inv (nx,ny,nz) ! defined as d( zp )/d( z ). ! Source in qv-equation due to cumulus real qvcumsrc(nx,ny,nz) ! parameterization ! Base state potential temperature (K) real ptbar (nx,ny,nz) ! Perturbation potential temperature real ptprt (nx,ny,nz,nt) ! Potential temperature at the ground le real ptsfc(nx,ny) ! Effective qv at the surface (kg/kg) real qvsfc(nx,ny) ! EXBC buffer size integer exbcbufsz ! EXBC buffer array real exbcbuf( exbcbufsz ) ! EXBC relaxation coefficients real bcrlx (nx,ny) ! Advection of water vapor (kg/(m**3*s)) real qadv(nx,ny,nz) ! A local array. ! Total mixing of water vapor real qmix(nx,ny,nz) ! (kg/(m**3*s)) ! A local array. ! Temporary work 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) ! Temporary work array. real tem2_0(0:nx,0:ny,0:nz) ! Temporary work array. real tem3_0(0:nx,0:ny,0:nz) ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i,j,k ! Starting level of time integration integer tstrtlvl real deltat integer ebc1,wbc1,nbc1,sbc1, qflag real cpu0, f_cputime ! !####################################################################### ! ! Include files: ! !####################################################################### ! include 'bndry.inc' include 'exbc.inc' include 'globcst.inc' ! !####################################################################### ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! !####################################################################### ! ! Compute the advection term of the water vapor specific humidity ! equation and store the result in array qadv. ! !####################################################################### ! cpu0 = f_cputime() ! Forward-based FCT scheme IF(sadvopt.eq.4) THEN 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 advection CALL advq(nx,ny,nz,qv,u,v,wcont,rhostr,qvbar, mapfct, qadv, & tem1,tem2,tem3,tem4,tem5,tem6) ! FCT advection ELSEIF( sadvopt .eq. 4.or.sadvopt.eq.5) THEN CALL advqfct(nx,ny,nz,dtbig1,qv,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 ! !####################################################################### ! ! Compute the mixing terms for the water vapor specific humidity ! equation. This includes both physical and computational ! mixing. Store the result in array qmix. ! !####################################################################### ! CALL mixqv(nx,ny,nz, & & qv(1,1,1,tstrtlvl),rhostr,qvbar,kmh,kmv,rprntl,qvsflx, & & pbldpth(1,1,tpresent), & & x,y,z,zp, mapfct, j1,j2,j3,j3inv, & & usflx,vsflx,ptsflx,ptsfc,qvsfc,ptbar,ptprt, & & qmix, & & tem1,tem2,tem3,tem4,tem5,tem6) ! !####################################################################### ! ! Call BRLXQ to calculate the boundary relaxation and computation ! mixing for qv ! ! qmix = qmix + qvexbc_term ! !####################################################################### ! IF ( lbcopt.eq.2 .and. mgrid.eq.1 ) THEN cpu0 = f_cputime() qflag = 1 CALL brlxq(nx,ny,nz, deltat*0.5, qflag, qv,rhostr, qmix, & exbcbuf(nqv0exb), exbcbuf(nqc0exb), exbcbuf(nqr0exb), exbcbuf( & nqi0exb), exbcbuf(nqs0exb), exbcbuf(nqh0exb), exbcbuf(nqvdtexb),& exbcbuf(nqcdtexb), exbcbuf(nqrdtexb),exbcbuf(nqidtexb), exbcbuf(& nqsdtexb),exbcbuf(nqhdtexb),bcrlx, tem1,tem2,tem3,tem4) bc_cpu_time = bc_cpu_time + f_cputime() - cpu0 ENDIF ! !####################################################################### ! ! Calculate qvforce, store the result in tem1. ! !####################################################################### ! 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)+rhostr(i,j,k)*qvcumsrc(i,j,k) & +qmix(i,j,k) 101 CONTINUE ! !####################################################################### ! ! Treat the vertically implicit mixing term. ! !####################################################################### ! ! Vertical implicit application IF (trbvimp.eq.1) THEN 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,qv(1,1,1,tstrtlvl),rhostr, & tem2,tem1,tem3,tem4,tem5,tem6) tmix_cpu_time = tmix_cpu_time + f_cputime() - cpu0 ENDIF ! !####################################################################### ! ! Integrate forward by one timestep the water vapor specific humidit ! equation, yielding qv at time = tfuture. ! !####################################################################### ! 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 qv(i,j,k,tfuture)=qv(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() ! !####################################################################### ! ! Set the boundary conditions on qv for an adaptive (nested) ! grid run. If using only one grid, this portion of the code is ! skipped....proceed to next comment block. ! !####################################################################### ! DO 100 k=2,nz-2 DO 100 i=1,nx-1 qv(i, 1,k,tfuture)= 2*qv(i, 1,k,tpresent)-qv(i, 1,k,tpast) qv(i,ny-1,k,tfuture)= 2*qv(i,ny-1,k,tpresent)-qv(i,ny-1,k, & tpast) 100 CONTINUE DO 105 k=2,nz-2 DO 105 j=2,ny-2 qv( 1,j,k,tfuture)= 2*qv( 1,j,k,tpresent)-qv( 1,j,k,tpast) qv(nx-1,j,k,tfuture)= 2*qv(nx-1,j,k,tpresent)-qv(nx-1,j,k, & tpast) 105 CONTINUE CALL bcqv(nx,ny,nz, deltat*0.5, qv,qvbar, qvdteb,qvdtwb,qvdtnb, & qvdtsb, 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 qv(i,j,k,tfuture)=qv(i,j,k,tstrtlvl) +deltat*tem1(i,j,k)/ & rhostr(i,j,k) 502 CONTINUE misc_cpu_time = misc_cpu_time + f_cputime() - cpu0 ! !####################################################################### ! ! Set the boundary conditions on qv for a NON-adaptive (uniform) ! grid run. ! !####################################################################### ! 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 bcqv(nx,ny,nz, dtbig1, qv,qvbar, qvdteb,qvdtwb,qvdtnb, & qvdtsb, ebc1,wbc1,nbc1,sbc1,tbc,bbc) ELSE CALL bcqv(nx,ny,nz, dtbig1, qv,qvbar, qvdteb,qvdtwb,qvdtnb, & qvdtsb, 3,3,3,3,tbc,bbc) qflag = 1 CALL exbcq(nx,ny,nz, qflag, curtim+dtbig,qv(1,1,1,tfuture), & exbcbuf(nqv0exb), exbcbuf(nqvdtexb)) ENDIF bc_cpu_time = bc_cpu_time + f_cputime() - cpu0 ENDIF RETURN END ! ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### 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. ! !####################################################################### ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! ! dtbig1 The big time step size (s) ! qflag Indicator for water/ice species ! ! q One of the liquid or ice variables at tpast and ! tpresent (kg/kg) ! ! u x component of velocity at all time levels (m/s) ! v y component of velocity at all time levels (m/s) ! wcont Contravariant vertical velocity (m/s) ! ! qdteb Time tendency of liquid/ice variable q at east boundary ! qdtwb Time tendency of liquid/ice variable q at west boundary ! qdtnb Time tendency of liquid/ice variable q at north ! boundary ! qdtsb Time tendency of liquid/ice variable q at south ! boundary ! ! rhostr Base state density rhobar times j3 (kg/m**3) ! ! kmh Horizontal turb. mixing coef. for momentum ( m**2/s ) ! kmv Vertical turb. mixing coef. for momentum ( m**2/s ) ! rprntl Reciprocal of Prandtl number ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp Vertical coordinate of grid points in physical space(m) ! ! mapfct Map factors at scalar points ! ! j1 Coordinate transformation Jacobian -d(zp)/dx ! j2 Coordinate transformation Jacobian -d(zp)/dy ! j3 Coordinate transformation Jacobian d(zp)/dz ! qccumsrc Source term in qc-equation due to cumulus parameterizat ! qrcumsrc Source term in qr-equation due to cumulus parameterizat ! qicumsrc Source term in qi-equation due to cumulus parameterizat ! qscumsrc Source term in qs-equation due to cumulus parameterizat ! ! ! OUTPUT: ! ! q Liquid/ice variable q at tfuture (kg/kg) ! ! WORK ARRAYS: ! ! qadv Advection term of water variable eq. (kg/(m**3*s)). ! A local array. ! qmix Total mixing in water variable eq. (kg/(m**3*s)). ! A local array. ! tem1 Temporary work array. ! tem2 Temporary work array. ! tem3 Temporary work array. ! tem4 Temporary work array. ! tem5 Temporary work array. ! tem6 Temporary work array. ! tem7 Temporary work array. ! tem8 Temporary work array. ! tem9 Temporary work array. ! tem10 Temporary work array. ! tem11 Temporary work array. ! ! tem1_0 Temporary work array. ! tem2_0 Temporary work array. ! tem3_0 Temporary work array. ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations: ! !####################################################################### ! ! Force explicit declarations implicit none ! The no. of t-levels of t-dependent integer nt ! arrays. ! Index of time level for the past time. integer tpast ! Index of time level for the present integer tpresent ! time. ! Index of time level for the future integer tfuture ! time. parameter (nt=3, tpast=1, tpresent=2, tfuture=3) ! Number of grid points in 3 directions integer nx, ny, nz ! Local big time step real dtbig1 ! Indicator for water/ice species integer qflag ! One of the water/ice variables (kg/kg) real q (nx,ny,nz,nt) ! Total u-velocity (m/s) real u (nx,ny,nz,nt) ! Total v-velocity (m/s) real v (nx,ny,nz,nt) ! Contravariant vertical velocity (m/s) real wcont (nx,ny,nz) ! Time tendency of liquid/ice at real qdteb(ny,nz) ! e-boundary ! Time tendency of liquid/ice at real qdtwb(ny,nz) ! w-boundary ! Time tendency of liquid/ice at real qdtnb(nx,nz) ! n-boundary ! Time tendency of liquid/ice at real qdtsb(nx,nz) ! s-boundary ! Base state density rhobar times j3. real rhostr(nx,ny,nz) ! Horizontal turb. mixing coef. for real kmh (nx,ny,nz) ! momentum. ( m**2/s ) ! Vertical turb. mixing coef. for real kmv (nx,ny,nz) ! momentum. ( m**2/s ) ! Reciprocal of Prandtl number real rprntl(nx,ny,nz) ! x-coord. of the physical and compu- real x (nx) ! tational grid. Defined at u-point. ! y-coord. of the physical and compu- real y (ny) ! tational grid. Defined at v-point. ! z-coord. of the computational grid. real z (nz) ! Defined at w-point on the staggered ! grid. ! Physical height coordinate defined at real zp (nx,ny,nz) ! w-point of the staggered grid. ! Map factors at scalar points real mapfct(nx,ny) ! Coordinate transformation Jacobian real j1 (nx,ny,nz) ! defined as - d( zp )/d( x ). ! Coordinate transformation Jacobian real j2 (nx,ny,nz) ! defined as - d( zp )/d( y ). ! Coordinate transformation Jacobian real j3 (nx,ny,nz) ! defined as d( zp )/d( z ). ! Coordinate transformation Jacobian real j3inv (nx,ny,nz) ! defined as d( zp )/d( z ). ! Source in qc-equation due to cumulus real qccumsrc(nx,ny,nz) ! parameterization ! Source in qr-equation due to cumulus real qrcumsrc(nx,ny,nz) ! parameterization ! Source in qi-equation due to cumulus real qicumsrc(nx,ny,nz) ! parameterization ! Source in qs-equation due to cumulus real qscumsrc(nx,ny,nz) ! parameterization ! EXBC buffer size integer exbcbufsz ! EXBC buffer array real exbcbuf( exbcbufsz ) ! EXBC relaxation coefficients real bcrlx (nx,ny) ! Advection of water/ice substance real qadv(nx,ny,nz) ! (kg/(m**3*s)). A local array. ! Total mixing of water/ice substance real qmix(nx,ny,nz) ! (kg/(m**3*s)). A local array. ! Temporary work 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) ! automatic 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) ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i,j,k integer tstrtlvl real deltat integer ebc1,wbc1,nbc1,sbc1 integer nq0exb,nqdtexb real cpu0, f_cputime ! !####################################################################### ! ! Include files: ! !####################################################################### ! include 'bndry.inc' include 'exbc.inc' include 'globcst.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! !####################################################################### ! ! Compute the advection term for a general water substance ! q and store the result in qadv. ! !####################################################################### ! cpu0 = f_cputime() ! Forward-based FCT scheme IF(sadvopt.eq.4) THEN 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) ! FCT advection ELSEIF( sadvopt.eq.4 .or. sadvopt.eq.5 ) THEN 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 ! !####################################################################### ! ! 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.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 ! !####################################################################### ! ! Calculate qforce, store the result in tem1. ! !####################################################################### ! 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 ! !####################################################################### ! ! Treat the vertically implicit mixing term. ! !####################################################################### ! ! Vertical implicit application IF (trbvimp.eq.1) THEN 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 ! !####################################################################### ! ! Integrate forward by one timestep the general water substance ! (q) equation, yielding q at time = tfuture. ! !####################################################################### ! 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() ! !####################################################################### ! ! 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 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 ! !####################################################################### ! ! 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.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 fo 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 ! ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE UPRAD1 ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! SUBROUTINE UPRAD1(nx,ny,w1,w2,wc,nrho,ifax1,ifax2,trigs1,trigs2, & & temxy1,work) ! !####################################################################### ! ! PURPOSE: ! ! To apply the upper radiation condition using a periodic 1 or 2-D ! fft. The variable temxy1 contains input information and is ! overwritten with the final w (at nz-1) on exit. ! ! In the case of the open boundary, a linear hydrostatic relation ! between vertical velocity and pressure is used in fourier space. ! The use of fast fourier transform is required for this condition. ! ! The fft routine used is FFT991 obtained from NCAR. It is a ! version of the one used by ECMWF. Before fft991 is called, set99 ! is called to setup the variables needed for the transform routine. ! (init3d.f) ! !####################################################################### ! ! AUTHOR: Dan Weber ! 07/22/1997. ! ! Modification History: ! !####################################################################### ! ! INPUT: ! ! nx Number of grid points in the x-direction. ! ny Number of grid points in the y-direction. ! ! w1 Pre-fft w(nz-1) coefficient. ! w2 Pre-fft w(nz-2) coefficient. ! wc Pre-fft w(nz-1) coefficient from the reduction ! phase of the tridiaginal solver. ! ! trigs1 Array containing pre-computed trig function for ! fftopt=1. ! trigs2 Array containing pre-computed trig function for ! fftopt=1. ! ifax1 Array containing the factors of nx for fftopt=1. ! ifax2 Array containing the factors of ny for fftopt=1. ! ! OUTPUT: ! ! ! temxy1 On input known forcing term, on output final w at nz-1. ! ! ! WORK ARRAYS: ! ! ! work 2-D work array for fft code. ! ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations. ! !####################################################################### ! implicit none ! Number of grid points in the x-directi integer nx ! Number of grid points in the y-directi integer ny ! pre-fft w(nz-1) coefficient. real w1 ! pre-fft w(nz-2) coefficient. real w2 ! pre-fft w(nz-1) coefficient from the real wc ! reduction phase of the tridiagonal sol ! Array containing pre-computed trig real trigs1(3*(nx-1)/2+1) ! function for fftopt=1. ! Array containing pre-computed trig real trigs2(3*(ny-1)/2+1) ! function for fftopt=1. ! Array containing the factors of nx for integer ifax1(13) ! fftopt=1. ! Array containing the factors of ny for integer ifax2(13) ! fftopt=1. ! On input known forcing term, real temxy1 (nx+1,ny+1) ! on output final w at nz-1. ! 2-D work array for fft code. real work (ny,nx) ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i,j,itema,itemb real eps,kx,ky,pi,nrho,tema,temb,temc,temd,teme,temf,temg ! !####################################################################### ! ! Include files: ! !####################################################################### ! include 'globcst.inc' include 'bndry.inc' !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ !####################################################################### ! ! Perform the fft on temxy1 ! !####################################################################### ! 1-d transform in x, note nx must be ODD! IF(ny.eq.4)THEN itema = 1 itemb = nx CALL fft991(temxy1,work,trigs1,ifax1,1,nx+1,nx-1,1,-1) ! 1-d transform in y, note ny must be ODD! ELSEIF(nx.eq.4)THEN itema = ny itemb = 1 CALL fft991(temxy1,work,trigs2,ifax2,nx+1,1,ny-1,1,-1) ! Do 2-d transform, note nx,ny must be ODD! ELSE itema = ny itemb = nx CALL fft991(temxy1,work,trigs1,ifax1,1,nx+1,nx-1,ny-1,-1) CALL fft991(temxy1,work,trigs2,ifax2,nx+1,1,ny-1,nx,-1) ! end of run type if block. END IF !####################################################################### ! ! Compute the wave space w at nz-1. ! !####################################################################### eps= 1.0e-8 pi = 4.0*atan(1.0) tema = pi/(nx-2) temb = pi/(ny-2) temc = 2.0*dxinv temd = 2.0*dyinv temf = w1-w2*wc temg = eps - nrho ! note: kx and ky are in finite difference form. DO 230 j=1,itema ky = temd*sin(int((j-1)*0.5+0.001)*temb) ky = ky*ky DO 230 i=1,itemb kx = temc*sin(int((i-1)*0.5+0.001)*tema) teme =sqrt(kx*kx+ky) temxy1(i,j)=-teme*temxy1(i,j)/(teme*temf+temg) 230 CONTINUE !####################################################################### ! ! Perform the reverse fft on temxy1 to obtain w(nz-1) in real space. ! !####################################################################### ! 1-d transform in x, note nx must be ODD! IF(ny.eq.4)THEN CALL fft991(temxy1,work,trigs1,ifax1,1,nx+1,nx-1,1,1) DO 240 j=2,ny-1 DO 240 i=1,nx-1 temxy1(i,j) = temxy1(i,1) 240 CONTINUE ! 1-d transform in y, note ny must be ODD! ELSEIF(nx.eq.4)THEN CALL fft991(temxy1,work,trigs2,ifax2,nx+1,1,ny-1,1,1) DO 250 j=1,ny-1 DO 250 i=2,nx-1 temxy1(i,j) = temxy1(1,j) 250 CONTINUE ! Do 2-d transform, note nx,ny must be ODD! ELSE CALL fft991(temxy1,work,trigs2,ifax2,nx+1,1,ny-1,nx, 1) CALL fft991(temxy1,work,trigs1,ifax1,1,nx+1,nx-1,ny-1, 1) ! end of run type if block. END IF RETURN END ! ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE UPRAD3 ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! SUBROUTINE UPRAD3(nx,ny,w1,w2,wc,nrho,wsave1,wsave2,vwork1,vwork2,& & temxy1,work) ! !####################################################################### ! ! PURPOSE: ! ! To apply the upper radiation condition using an even (vfftpack) ! sequence fft. The variable temxy1 contains input information and ! is overwritten with the updated w (at nz-1) on exit. ! ! In the case of the open boundary, a linear, hydrostatic relation ! between vertical velocity to pressure is used in fourier space. ! The use of fast fourier transform is NOT required for this ! option but is recommended. If nx and ny are not of special ! character, then a basic fourier transform is used. ! ! The fft routine used is vcost from the vfftpack software suite ! available at PSC. ! ! Before vcost is called, the arrays wsave1,wsave2,vwork1,vwork2 ! must be initialized (init3d.f). ! !####################################################################### ! ! AUTHOR: Dan Weber ! 07/22/1997. ! ! ! Modification History: ! !####################################################################### ! ! INPUT: ! ! nx Number of grid points in the x-direction. ! ny Number of grid points in the y-direction. ! ! w1 Pre-fft w(nz-1) coefficient. ! w2 Pre-fft w(nz-2) coefficient. ! wc Pre-fft w(nz-1) coefficient from the reduction ! phase of the tridiaginal solver. ! ! vwork1 2-D work array for fftopt=2 option. ! vwork2 2-D work array for fftopt=2 option. ! wsave1 Work array for fftopt=2 option. ! wsave2 Work array for fftopt=2 option. ! ! OUTPUT: ! ! ! temxy1 On input known forcing term, on output final w at nz-1. ! ! ! WORK ARRAYS: ! ! ! work 2-D work array for fft code. ! ! !####################################################################### ! ! !####################################################################### ! ! Variable Declarations. ! !####################################################################### implicit none ! Number of grid points in the x integer nx ! direction. ! Number of grid points in the y integer ny ! direction. ! Pre-fft w(nz-1) coefficient. real w1 ! Pre-fft w(nz-2) coefficient. real w2 ! Pre-fft w(nz-1) coefficient from the real wc ! reduction phase of the tridiagonal ! solver. ! 2-D work array for fftopt=2. real vwork1 (nx+1,ny+1) ! 2-D work array for fftopt=2. real vwork2 (ny,nx+1) ! Work array for fftopt=2. real wsave1 (3*(ny-1)+15) ! Work array for fftopt=2. real wsave2 (3*(nx-1)+15) ! On input known forcing term, real temxy1 (nx+1,ny+1) ! on output final w at nz-1. ! 2-D work array (global) real work (ny,nx) ! !####################################################################### ! ! Include files: ! !####################################################################### ! include 'globcst.inc' include 'bndry.inc' ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i,j real kx,ky,pi,tema,temb,temc,temd,teme,temf real eps,nrho ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ !####################################################################### ! ! Perform the fft on temxy1. ! !####################################################################### ! 1-d transform in x. IF(ny.eq.4)THEN ! swap temxy1(i,j) array into work(j,i) DO 10 i=1,nx-1 work(1,i) = temxy1(i,1) 10 CONTINUE CALL vcost(1,nx-1,work,vwork2,ny,wsave2) ! swap back for computations.... DO 20 i=1,nx-1 temxy1(i,1) = work(1,i) 20 CONTINUE ! 1-d transform in y. ELSEIF(nx.eq.4)THEN CALL vcost(1,ny-1,temxy1,vwork1,nx+1,wsave1) ! Do 2-d transform. ELSE CALL vcost(nx-1,ny-1,temxy1,vwork1,nx+1,wsave1) ! swapping arrays. DO 50 j=1,ny-1 DO 50 i=1,nx-1 work(j,i) = temxy1(i,j) 50 CONTINUE CALL vcost(ny-1,nx-1,work,vwork2,ny,wsave2) ! swapping arrays for computations. DO 60 j=1,ny-1 DO 60 i=1,nx-1 temxy1(i,j) = work(j,i) 60 CONTINUE ! end of run type if block. END IF !####################################################################### ! ! Compute the wave space w at nz-1. ! !####################################################################### eps = 1.0e-13 pi = 4.0*atan(1.0) ! finite difference form. tema = 0.5*pi/(nx-2) ! finite difference form. temb = 0.5*pi/(ny-2) temc = 2.0*dxinv temd = 2.0*dyinv teme = w1-w2*wc temf = eps - nrho ! first part of 2-d case. temxy1(1,1) = -temxy1(1,1)/(teme+temf) ! perform computations in x direction. IF(ny.eq.4)THEN ! compute from i=2,nx-1 and j=1... DO 70 i=2,nx-1 kx = temc*sin(int(i-1)*tema) kx = sqrt(kx*kx) temxy1(i,1)=-kx*temxy1(i,1)/(kx*teme+temf) 70 CONTINUE ! perform computations in y direction. ELSEIF(nx.eq.4)THEN ! compute from j=2,ny-1 and i=1. DO 80 j=2,ny-1 ky = temd*sin(int(j-1)*temb) ky = sqrt(ky*ky) temxy1(1,j)=-ky*temxy1(1,j)/(ky*teme+temf) 80 CONTINUE ! compute the 2-D version. ELSE ! compute from i=2,nx-1 and j=1. DO 85 i=2,nx-1 kx = temc*sin(int(i-1)*tema) kx = sqrt(kx*kx) temxy1(i,1)=-kx*temxy1(i,1)/(kx*teme+temf) 85 CONTINUE ! compute from i=1,j=2,ny-1. DO 90 j=2,ny-1 ky = temd*sin(int(j-1)*temb) ky = sqrt(ky*ky) temxy1(1,j)=-ky*temxy1(1,j)/(ky*teme+temf) 90 CONTINUE ! compute from i=2,nx-1,j=2,ny-1. DO 100 j=2,ny-1 ky = temd*sin(int(j-1)*temb) ky = ky*ky DO 100 i=2,nx-1 kx = temc*sin(int(i-1)*tema) kx =sqrt(kx*kx+ky) temxy1(i,j)=-kx*temxy1(i,j)/(kx*teme+temf) 100 CONTINUE ! end of the run type if block. END IF !####################################################################### ! ! Perform the reverse fft on temxy1 to obtain w in real space. ! !####################################################################### ! 1-d transform in x, note nx must be EVEN! IF(ny.eq.4)THEN DO 110 i=1,nx-1 work(1,i) = temxy1(i,1) 110 CONTINUE CALL vcost(1,nx-1,work,vwork2,ny,wsave2) DO 120 j=1,ny-1 DO 120 i=1,nx-1 temxy1(i,j) = work(1,i) 120 CONTINUE ! 1-d transform in y, note ny must be EVEN! ELSEIF(nx.eq.4)THEN CALL vcost(1,ny-1,temxy1,vwork1,nx+1,wsave1) DO 140 j=1,ny-1 DO 140 i=1,nx-1 temxy1(i,j) = temxy1(1,j) 140 CONTINUE ! Do 2-d transform, note nx,ny must be EVEN! ELSE ! swapping arrays. DO 150 j=1,ny-1 DO 150 i=1,nx-1 work(j,i) = temxy1(i,j) 150 CONTINUE CALL vcost(ny-1,nx-1,work,vwork2,ny,wsave2) ! swapping arrays. DO 160 j=1,ny-1 DO 160 i=1,nx-1 temxy1(i,j) = work(j,i) 160 CONTINUE CALL vcost(nx-1,ny-1,temxy1,vwork1,nx+1,wsave1) ! end of runopt if block, w (nz-1) is updated. END IF RETURN END ! ! ! ################################################################## ! ################################################################## ! ###### ###### ! ###### SUBROUTINE PPRTBBC ###### ! ###### ###### ! ###### Developed by ###### ! ###### Center for Analysis and Prediction of Storms ###### ! ###### University of Oklahoma ###### ! ###### ###### ! ################################################################## ! ################################################################## ! SUBROUTINE PPRTBBC(nx,ny,nz,g05,buoy2swt,rhostr,pprt,ptprt, & & pbar,ptbar,phydro, & & tem1) ! !####################################################################### ! ! PURPOSE: ! ! To compute the perturbation pressure below the ground (scalar k=1) ! The technique uses the perturbation hydrostatic relation obtained ! from the model vertical momentum equation. The technique is ! summarized here. ! ! The calculation of the hydrostatic pressure bottom boundary ! condition involves two steps. The first step is used to ! approximate the second order pprt term in the buoyancy term ! and the second step obtains the final pprt at k=1. ! ! Step A: pprt(1)=D=(pprt(2) - F - dz*phydro + 2.0*A*B)/(1-C) ! ! where, A = dz*g*avgz(rhobar*J3)/(2*cpdcv) ! B = pprt(i,j,2)/pbar(i,j,2) ! C = A/pbar(i,j,1) ! F = alpha*div*(1)-alpha*div*(2) not included.... ! note: ! D = the intermediate result pprt(1) for use in the ! second order pprt term ! ! Step B: pprt(1) = D + A*(1/cpdcv -1)*(B*B + E*E)/(1.0-C) ! ! where, E = D/pbar(i,j,1) ! ! The result is stored in phydro(i,j) and passed into bcp to set ! pprt(i,j,1). ! ! !####################################################################### ! ! AUTHOR: Dan Weber ! 11/06/1997. ! ! ! Modification History: ! !####################################################################### ! ! !####################################################################### ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! ! ptprt Perturbation potential temperature at time tpresent (K) ! pprt Perturbation pressure at time tpresent (Pascal) ! ! phydro Big time step forcing term for use in computing the ! hydrostatic pressure at k=1. ! ! rhostr Base state density rhobar times j3 (kg/m**3) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! ! OUTPUT: ! ! tem1 Holds pprt(i,j,1) via hydrostatic balance. ! !####################################################################### ! ! Variable Declarations. ! !####################################################################### ! ! Force explicit declarations implicit none ! Number of grid points in 3 directions integer nx, ny, nz ! Perturbation potential temperature (K) real ptprt (nx,ny,nz) ! Perturbation pressure at tfuture real pprt (nx,ny,nz) ! (Pascal) ! Big time step forcing for computing real phydro(nx,ny) ! hydrostatic pprt at k=1. ! Base state density rhobar times j3. real rhostr(nx,ny,nz) ! Base state potential temperature (K) real ptbar (nx,ny,nz) ! Base state pressure (Pascal). real pbar (nx,ny,nz) ! !####################################################################### ! ! Temporary WORK ARRAYS: ! !####################################################################### ! ! Temporary work array real tem1 (nx,ny,nz) ! !####################################################################### ! ! Misc. local variables: ! !####################################################################### ! integer i, j real g05 real tema,temb,a,b,c,d,e !Switch for 1st-order or 2nd-order in buoyancy integer buoy2swt real ptemk,ptemk1,pttemk,pttemk1 ! !####################################################################### ! ! Include files: ! !####################################################################### ! include 'bndry.inc' include 'globcst.inc' ! Physical constants include 'phycst.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! add in the 1st and 2nd order pt terms. IF( ptsmlstp.eq.1 ) THEN tema = 0.5/cpdcv DO 100 j=1,ny-1 DO 100 i=1,nx-1 ptemk = pprt(i,j,2)*rhostr(i,j,2)/pbar(i,j,2) ptemk1= pprt(i,j,1)*rhostr(i,j,1)/pbar(i,j,1) pttemk = ptprt(i,j,2)*rhostr(i,j,2)/ptbar(i,j,2) pttemk1= ptprt(i,j,1)*rhostr(i,j,1)/ptbar(i,j,1) tem1(i,j,1)=phydro(i,j)+g05*(pttemk+pttemk1 -buoy2swt*(tema*& (ptemk*pttemk+ptemk1*pttemk1) +pttemk*pttemk+pttemk1* & pttemk1)) 100 CONTINUE ELSE DO 110 j=1,ny-1 DO 110 i=1,nx-1 tem1(i,j,1)= phydro(i,j) 110 CONTINUE ENDIF tema =buoy2swt*0.5*(1.0/cpdcv -1.0) temb = dz*g*0.25/cpdcv DO 120 j=1,ny-1 DO 120 i=1,nx-1 ! compute the intermediate result... A = temb*(rhostr(i,j,2)+rhostr(i,j,1)) B = pprt(i,j,2)/pbar(i,j,2) C = A/pbar(i,j,1) D = (pprt(i,j,2) - dz*tem1(i,j,1) + A*B) / (1.0-C) ! compute the final pprt(i,j,1).... E = D/pbar(i,j,1) tem1(i,j,1) = D + A*tema*(B*B + E*E) / (1.0-C) 120 CONTINUE CALL bcs2d(nx,ny,tem1(1,1,1), ebc,wbc,nbc,sbc) RETURN END