subroutine hadvppm(nn,dt,dx,con,vel,area,areav,flxarr,mynn) c c----CAMx v7Beta6 190902 c c HADVPPM performs advection using the one-dimensional implementation c of the piecewise parabolic method of Colella and Woodward (1984). c A piecewise continuous parabola is used as the intepolation polynomial. c The slope of the parabola at cell edges is computed from a cumulative c function of the advected quantity. These slopes are further modified c so that the interpolation function is monotone. c c This version based on CMAQ HPPM.F, v 1.1.1.1 9/14/98, written by c M.T. Odman (10/5/93), NCSC. This version assumes constant grid cell c size. c c The following definitions are used: c c |-----------> Positive direction c c |Boundary|<-----------------Domain----------------->|Boundary| c c | CON(1) | CON(2) | ... | CON(I) | ... |CON(N-1)| CON(N) | c c VEL(1)-->| VEL(I-1)-->| |-->VEL(I) |-->VEL(N-1) c c FP(1)-->| FP(I-1)-->| |-->FP(I) |-->FP(N-1) c c FM(2)<--| FM(I)<--| |<--FM(I+1) |<--FM(N) c c -->| DX |<-- c c Copyright 1996 - 2018 c Ramboll c c Modifications: c 5/17/00 small modification to flux1,2 to improve mass accounting c 10/13/03 area weighting applied to winds rather than conc c 11/04/03 only single FLXARR vector passed back for use in c ZRATES and probling tools c 9/1/09 Revised area weighting, removed from wind vector c c Input arguments: c nn Number of cells c dt Time step (s) c dx Length of cell (m) c con Concentration vector (umol/m3) c vel Wind speed vector (m/s) c area Cell area adjustment vector (1/m2) c areav Interfacial area adjustment vector (m2) c c Output arguments: c con Concentration vector (umol/m3) c flxarr Conc change from interfacial mass flux (umol/m3) c c Routines called: c none c c Called by: c XYADVEC c ZRATES c include "camx.prm" c integer mynn real con(mynn),vel(mynn),area(mynn),areav(mynn),flxarr(mynn) c c-----Local parameters c STEEPEN is a flag for discontinuty capturing (steepening) c This is disabled in this version as marked by c*** c c*** logical STEEPEN c*** parameter (STEEPEN=.false.) c*** parameter (ETA1=20.0, ETA2=0.05, EPS=0.01) parameter (TWO3RDS=2./3.) c c*** real fm(MX1D),fp(MX1D),cm(MX1D),cl(MX1D),cr(MX1D),dc(MX1D), c*** & c6(MX1D),d2c(MX1D),eta(MX1D),etabar(MX1D),cld(MX1D),crd(MX1D) c real fm(MXCELLS) real fp(MXCELLS) real cm(MXCELLS) real cl(MXCELLS) real cr(MXCELLS) real dc(MXCELLS) real c6(MXCELLS) c c-----Entry point c c-----Set all fluxes to zero. Either positive or negative flux will c remain zero depending on the sign of the velocity c c*** zeta = dx*dx do i = 1,nn fm(i) = 0. fp(i) = 0. enddo c c-----Zero order polynomial at the boundary cells c cm(2) = con(2) cm(nn) = con(nn-1) c c-----First order polynomial at the next cells, no monotonicity constraint c needed c cm(3) = (con(3) + con(2))/2. cm(nn-1) = (con(nn-1) + con(nn-2))/2. c c-----Second order polynomial inside the domain c do i = 3,nn-2 c c-----Compute average slope in the i'th cell c dc(i) = 0.5*(con(i+1) - con(i-1)) c c-----Guarantee that CM lies between CON(I) and CON(I+1) c monotonicity constraint if ((con(i+1) - con(i))*(con(i) - con(i-1)).gt.0.) then dc(i) = sign(1.,dc(i))*min( & abs(dc(i)), & 2.*abs(con(i+1) - con(i)), & 2.*abs(con(i) - con(i-1))) else dc(i) = 0. endif enddo c do i = 3,nn-3 cm(i+1) = con(i) + & 0.5*(con(i+1) - con(i)) + (dc(i) - dc(i+1))/6. enddo c do i = 2,nn-1 cr(i) = cm(i+1) cl(i) = cm(i) enddo c c-----Optional discontinuty capturing c This is disbaled completely in this version c c*** if (STEEPEN) then c*** do i = 2,nn-1 c*** eta(i) = 0. c*** cld(i) = con(i) c*** crd(i) = con(i) c*** enddo c***c c***c-----Finite diff approximation to 2nd derivative c***c c*** do i = 3,nn-2 c*** d2c(i) = (con(i+1) - 2.*con(i) + con(i-1))/6. c*** enddo c***c c***c-----No discontinuity detection near the boundary: cells 2, 3, NN-2, NN-1 c***c c*** do i = 4,nn-3 c***c c***c-----Compute etabars c***c c*** if ((-d2c(i+1)*d2c(i-1).gt.0.) .and. c*** & (abs(con(i+1) - con(i-1)) - c*** & EPS*min(abs(con(i+1)),abs(con(i-1))).gt.0.)) then c*** etabar(i) = -zeta*(d2c(i+1) - d2c(i-1))/ c*** & (con(i+1) - con(i-1)) c*** else c*** etabar(i) = 0. c*** endif c*** eta(i) = max(0.,min(ETA1*(etabar(i) - ETA2),1.)) c*** crd(i) = con(i+1) - 0.5*dc(i+1) c*** cld(i) = con(i-1) + 0.5*dc(i-1) c*** enddo c***c c*** do i = 2,nn-1 c*** cr(i) = cm(i+1) + eta(i)*(crd(i) - cm(i+1)) c*** cl(i) = cm(i) + eta(i)*(cld(i) - cm(i)) c*** enddo c*** endif c c-----Generate piecewise parabolic distributions c do i = 2,nn-1 c c-----Monotonicity c if ((cr(i) - con(i))*(con(i) - cl(i)).gt.0.) then dc(i) = cr(i) - cl(i) c6(i) = 6.*(con(i) - 0.5*(cl(i) + cr(i))) c c-----Overshoot cases c if (dc(i)*c6(i) .gt. dc(i)*dc(i)) then cl(i) = 3.*con(i) - 2.*cr(i) elseif (-dc(i)*dc(i) .gt. dc(i)*c6(i)) then cr(i) = 3.*con(i) - 2.*cl(i) endif else cl(i) = con(i) cr(i) = con(i) endif dc(i) = cr(i) - cl(i) c6(i) = 6.*(con(i) - 0.5*(cl(i) + cr(i))) enddo c c-----Compute fluxes from the parabolic distribution c do i = 2,nn-1 x = max(0., -vel(i-1)*dt/dx) fm(i) = x*(cl(i) + 0.5*x*(dc(i) + c6(i)*(1. - TWO3RDS*x))) x = max(0., vel(i)*dt/dx) fp(i) = x*(cr(i) - 0.5*x*(dc(i) - c6(i)*(1. - TWO3RDS*x))) enddo c c-----Compute fluxes from boundary cells assuming uniform distribution c if (vel(1).gt.0.) then x = vel(1)*dt/dx fp(1) = x*con(1) endif c if (vel(nn-1).lt.0.) then x = -vel(nn-1)*dt/dx fm(nn) = x*con(nn) endif c c-----Update concentrations c flxarr(1) = (fp(1) - fm(2)) do i = 2,nn-1 flxarr(i) = (fp(i) - fm(i+1)) con(i) = con(i) - & area(i)*(areav(i)*flxarr(i) - areav(i-1)*flxarr(i-1)) enddo c return end