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