!*******************************************************************************
!$RCSfile$
!$Revision$
!$Date$
!*******************************************************************************
SUBROUTINE set_stack_rel_prime( frac )
!*******************************************************************************
!
! FUNCTION:  Set true release parameters for STACK release
!            Use PRIME
!
! PRECONDITIONS REQUIRED:
!
! SUBROUTINES AND FUNCTIONS CALLED:
!
! REVISION HISTORY:
!
!*******************************************************************************

! --- MODULES

USE sciprime_fi
USE scipuff_fi
USE met_fi
USE files_fi
USE error_fi
USE constants_fd
USE sampler_fi

IMPLICIT NONE

! --- ARGUMENTS

REAL, INTENT( OUT ) :: frac    ! - puff height (includes plume rise on output)
                             !   mass fraction in primary plume

! --- LOCALS

CHARACTER(20)             :: rdfrm
REAL,PARAMETER            :: gravi = 9.807

REAL                      :: xmlinv
REAL                      :: ddz, thv1, thv2
REAL                      :: xmap, ymap    ! - map factors

INTEGER                   :: ios, alloc_stat
INTEGER                   :: i,k,ifvsec,isrc
INTEGER                   :: nza,ntr
REAL,DIMENSION(nst)       :: dturb,dtrur,purb,prur
REAL,DIMENSION(mxnz)      :: uamb,ramb,tamb,zgpta   ! common data in ambient
REAL,DIMENSION(mxnzp1)    :: dedz,zfacea      ! common data in ambient
REAL,DIMENSION(mxntr)     :: xtr,ytr,ztr,rtr

REAL                      :: xr,yr,zr,ybadj,sz,sy,szc,syc
REAL                      :: dszdx,dsydx,dszcdx,dsycdx
REAL                      :: h,hx,hy,tamb0,ramb0,fqcav, hsav
!!!!REAL                      :: dsbh,dsbw,dsbl,xadj,yadj
!      DSBH - real    - Effective building height (m)
!      DSBW - real    - Effective building width (m) across flow
!      DSBL - real    - Effective building length (m) along flow
!      XADJ - real    - Distance (m) from source to upwind face of bldg along flow
!      YADJ - real    - Distance (m) from source to center of upwind face of bldg across flow

LOGICAL                   :: nopath, nokey
LOGICAL                   :: lfound

!***********************************************************************
!     Initialize Default Wind Profile Exponents and DTHETADZ
!***********************************************************************

!STAB.CLASS  A    B     C     D      E      F
!           ***  ***   ***   ***    ***    ***
DATA dturb  /0.,  0.,   0.,   0.,   0.02, 0.035/, &
     dtrur  /0.,  0.,   0.,   0.,   0.02, 0.035/, &
     purb  /0.15, 0.15, 0.20, 0.25, 0.30, 0.30/,  &
     prur  /0.07, 0.07, 0.10, 0.15, 0.35, 0.55/

IF( nsrc_prime == 0 )THEN

  !-- Initialise variables
  CALL init_stack_rel_prime()

  iline = 0
  WRITE(rdfrm,'("(A",I3.3,",T1,",I3.3,"A1)")') istrg, istrg
  OPEN(lun_tmp,FILE=TRIM(relDisplay),STATUS='OLD',ACTION='READ',IOSTAT=ios)
  IF( ios /= 0 )THEN
    nError   = IV_ERROR
    eRoutine ='set_stack_rel_prime'
    eMessage = 'Missing Prime input File: '//TRIM(relDisplay)
    GOTO 9999
  END IF
  DO
    ! Increment the Line Counter
    iline = iline + 1
    ! READ Record to Buffers, as A80 and 80A1 for ISTRG = 80.
    ! Length of ISTRG is Set in PARAMETER Statement in MAIN1.INC
    READ(lun_tmp,rdfrm,IOSTAT=ios) runst1, (runst(i), i = 1, istrg)
    IF( ios == 0 )THEN
      ! Convert Lower Case to Upper Case Letters           ---   CALL LWRUPR
      CALL lwrupr
      ! Define Fields on Card                              ---   CALL DEFINE
      CALL define
      ! Get the Contents of the Fields                     ---   CALL GETFLD
      CALL getfld
      !  If Blank Line, Then CYCLE to Next Card
      IF (bline) CYCLE
      ! Extract Pathway ID From Field 1                    ---   CALL EXPATH
      CALL expath(field(1),nopath)
      ! For Invalid Pathway and Comment Lines Skip to Next Record
      IF (nopath .or. path == '**') CYCLE
      ! Extract Keyword From Field 2                       ---   CALL EXKEY
      CALL exkey(field(2),nokey)
      ! When Keyword Is Wrong, Save Keyword and Skip To The Next Record
      IF (nokey) THEN
        pkeywd = keywrd
        CYCLE
      END IF
      ! Process Input Card Based on Pathway
      IF (path == 'SO') THEN
        ! Process SOurce Pathway Cards                    ---   CALL SOCARD
        CALL socard
      END IF
      ! Store the Current Keyword as the Previous Keyword
      pkeywd = keywrd
    ELSE
      EXIT
    END IF
  END DO
  CLOSE(lun_tmp,IOSTAT=ios)
  ! Check Source Array Limits for Too Few Values
  CALL srcqa
  IF (nError /= NO_ERROR) GO TO 9999
  ! Check for missing keywords
  DO i = 1,ikn
    IF (isstat(i) == 0) THEN
      nError = WN_ERROR
      eRoutine='set_stack_rel_prime'
      WRITE(eInform,*)'Using default value of 0.0'
      SELECT CASE(i)
        CASE(1)
          eMessage= 'Missing rural keyword'
          WRITE(eInform,*)'Using default/previously set value of ',rural
        CASE(2)
          eMessage= 'Missing BUILDHGT keyword'
        CASE(3)
          eMessage= 'Missing BUILDWID keyword'
        CASE(4)
          eMessage= 'Missing BUILDLEN keyword'
        CASE(5)
          eMessage= 'Missing XBADJ keyword'
        CASE(6)
          eMessage= 'Missing YBADJ keyword'
       END SELECT
       eAction  = CHAR(0)
       CALL WarningMessage(.true.)
       IF (nError /= NO_ERROR) GO TO 9999
    END IF
  END DO

  IF( nsmp > 0 )THEN
    ALLOCATE( concPrime(nsmp),STAT=alloc_stat )
    IF( alloc_stat /= 0 )THEN
      nError   = UK_ERROR
      eRoutine = 'set_stack_rel_prime'
      eMessage = 'Error allocating PRIME receptor array'
      GOTO 9999
    END IF
    concPrime = 0.
  END IF

END IF

CALL numpr1
CALL prime1

hstack    = zrel                          ! stack height (m)
!dstack    = size_rel                      ! stack diameter (m)
dstack    = 2.*sigx                       ! stack diameter (m) - DSH 2012/05/18 because size_rel=0
tstack    = buoy + 273.15                 ! exhaust temperature (deg K)
wstack    = wmom                          ! exhaust velocity (m/s)

! KST - integer - PG stability class
IF (xml /= 0.0) THEN
   xmlinv = 1./xml
   IF (xmlinv < -0.14) THEN
      kst = 1
   ELSE IF (xmlinv < -0.05) THEN
      kst = 2
   ELSE IF (xmlinv < -0.01) THEN
      kst = 3
   ELSE IF (xmlinv < 0.01) THEN
      kst = 4
   ELSE IF (xmlinv < 0.06) THEN
      kst = 5
   ELSE
      kst = 6
   END IF
ELSE
   kst = 4
END IF

CALL sindex(srcid,nsrc_prime,relName,isrc,lfound)
IF( .NOT.lfound )THEN
  eRoutine = 'set_stack_rel_prime'
  nError   = IV_ERROR
  eMessage = 'PRIME source not found'
  eInform  = 'Release name = '//TRIM(relName)
  GO TO 9999
END IF

IF (isrural(isrc)) THEN
    rural = .true.
    urban = .false.
ELSE
    urban = .true.
    rural = .false.
END IF

! Power law coefficient and temperature gradient
IF (rural) THEN
    pvel   = prur(kst)
    dtemdz = dtrur(kst)
ELSE
    pvel   = purb(kst)
    dtemdz = dturb(kst)
END IF

! ---- Set stack parameters
CALL get_met( xrel,yrel,zrel,0.,0.,1 )
CALL get_sector(ub,vb,uref,ifvsec,afv)
dsbh = adsbh(ifvsec,isrc)
dsbw = adsbw(ifvsec,isrc)
dsbl = adsbl(ifvsec,isrc)
xadj = adsxadj(ifvsec,isrc)
yadj = adsyadj(ifvsec,isrc)
CALL wakflg(hstack,dsbh,dsbw,wake)

IF( wake )THEN
  tamb0  = tb
   ! --- Set meteorological arrays for plume rise calculation
  IF (primeprof) THEN
    CALL nummet_prime(uref,zrel,pvel,tamb0,dtemdz,ldbhr)
    !ustack - Wind speed (m/s) at release height
    CALL wsadj(ustack,hstack,uref,zrel,pvel)
    !ubldg - Wind speed (m/s) at top of building
    CALL wsadj(ubldg,dsbh,ustack,zrel,pvel)
    CALL set_bl_prime(-1.,0.)
  ELSE
    CALL get_prime_grid(nza,mxnz,zgpta,zfacea,ramb0)
    CALL get_met(xrel,yrel,2.,0.,0.,1)
    thv1   = (1. + 0.608*hb)*thb
    CALL get_met(xrel,yrel,zgpta(1),0.,0.,1)
    uamb(1)  = SQRT(ub*ub + vb*vb)
    tamb(1)  = tb
    thv2     = (1. + 0.608*hb)*thb
    dedz(1)  = (thv2-thv1)/zgpta(1)
    CALL set_bl_prime(wts,zinv)
    DO k = 2, nza
      ddz = 1./(zgpta(k) - zgpta(k-1))
      thv1 = thv2
      CALL get_met(xrel,yrel,zgpta(k),0.,0.,1)
      uamb(k)  = SQRT(ub*ub + vb*vb)
      tamb(k)  = tb
      thv2     = (1. + 0.608*hb)*thb
      dedz(k) = ddz*(thv2-thv1)
    END DO
    dedz(nza+1) = dedz(nza)
    ! Create ambient density profile on grid center
    ! (Code section from nummet of PRIME.FOR
    DO i = 1,nza
      ! Compute average temperature in layer
      thv1 = 0.5*(tamb(i)+tamb0)
      ! Compute ambient air density at height, ZGPTA(i)
      ramb(i)=ramb0*(tamb0/tamb(i))*EXP(-gravi*zgpta(i)/(RGAS*thv1))
    ENDDO
    CALL set_ambient(nza,uamb,ramb,dedz,tamb,tamb0,ramb0)
    CALL get_met(xrel,yrel,hstack,0.,0.,1)
    ustack = MAX(uMin,SQRT(ub*ub+vb*vb)) ! Wind speed (m/s) at release height
    CALL get_met(xrel,yrel,dsbh,0.,0.,1)
    ubldg  = MAX(uMin,SQRT(ub*ub+vb*vb))  ! Wind speed (m/s) at top of building
  END IF
  CALL wake_ini(ldbhr,kst,rural,dsbh,dsbw,dsbl,xadj,yadj,ubldg,ustack)
  hseff = hstack
  reff  = 0.5*dstack
  ntr = mxntr
  CALL numrise(ldbhr,hseff,reff,tstack,wstack,ntr,xtr,ytr,ztr,rtr)
  CALL get_wakedat(xr,yr,zr,ybadj,sy,sz,dsydx,dszdx,syc,szc,dsycdx, &
                   dszcdx,ntr,xtr,ytr,ztr,hseff,fqcav,wmom, &
                   buoy,nError)
  IF (nError /= NO_ERROR) THEN
    eRoutine= 'get_wakedat'
    eMessage= 'xr Point lies outside range of tabulated values '
    WRITE(eInform,*)'PRIME Model not implemented'
    eAction  = 'Using Plume rise algorithm'
    CALL WarningMessage(.true.)
    wake = .false.
    GO TO 9999
  END IF
! Reset release location
  CALL mapfac( xrel , yrel , xmap , ymap )
  wdsin = SIN(afv)
  wdcos = COS(afv)
  xrel_prm(1) = xrel + xr*wdsin*xmap - yr*wdcos*xmap  ! afv = -pi to +pi from North.
  xrel_prm(2) = xrel + xr*wdsin*xmap - (yr+ybadj)*wdcos*xmap
  yrel_prm(1) = yrel + xr*wdcos*ymap + yr*wdsin*xmap
  yrel_prm(2) = yrel + xr*wdcos*ymap + (yr+ybadj)*wdsin*xmap

! set scichem release parameters for the two source
  if (lter) then
    call get_topog(xrel_prm(1),yrel_prm(1),h,hx,hy)
    if (nError /= NO_ERROR) go to 9999
  else
    h = 0.0
  end if
  hsav = h
  zrel_prm(1)  = zr  + h
  if (lter) then
    call get_topog(xrel_prm(2),yrel_prm(2),h,hx,hy)
    if (nError /= NO_ERROR) go to 9999
  endif
  zrel_prm(2)  = h
  sigy_prm(1)  = sy
  sigy_prm(2)  = syc
  sigz_prm(1)  = sz
  sigz_prm(2)  = szc
  frac_prm(1)  = 1.0 - fqcav
  frac_prm(2)  = fqcav
  cmass_prm    = cmass
  wmom_prm(1)  = wmom
  wmom_prm(2)  = 0.0
  buoy_prm(1)  = buoy
  buoy_prm(2)  = 0.0
  !   Reset bouyancy to zero
  wmom        = 0.0
  buoy        = 0.0
  size_rel    = 0.0
  frac        = frac_prm(1)
  CALL get_met(xrel_prm(2),yrel_prm(2),zrel_prm(2),0.,0.,0)
  ustack    = MAX(uMin,SQRT(ub*ub+vb*vb)) ! Wind speed (m/s) at secondary src release location
  ky_prm(2) = syc*dsycdx*ustack
  kz_prm(2) = szc*dszcdx*ustack

  CALL get_met(xrel_prm(1),yrel_prm(1),zrel_prm(1),0.,0.,0)
  ustack    = MAX(uMin,SQRT(ub*ub+vb*vb)) ! Wind speed (m/s) at primary src release height ** saved for receptors **
  ky_prm(1) = sy*dsydx*ustack
  kz_prm(1) = sz*dszdx*ustack

  zi_rel = zinv - hsav

  if (fqcav > 0.0) then
     write(lun_log,'("Prime calculation: Source 1(xrel,yrel,zrel,1-fqcav,cmass_prm): ")',ADVANCE='NO')
     write(lun_log,'(3(1pE10.3,1x),0pF4.2,1x,1pE10.3)')xrel_prm(1), yrel_prm(1), zrel_prm(1),frac_prm(1),cmass_prm
     write(lun_log,'("Prime calculation: Source 2(xrel,yrel,zrel,fqcav): ")',ADVANCE='NO')
     write(lun_log,'(3(1pE10.3,1x),0pF4.2)')xrel_prm(2),yrel_prm(2), zrel_prm(2),frac_prm(2)
  else
     write(lun_log,'("Prime calculation (xrel,yrel,zrel,cmass_prm): ")',ADVANCE='NO')
     write(lun_log,'(3(1pE10.3,1x),1pE10.3)')xrel_prm(1), yrel_prm(1), zrel_prm(1),cmass_prm
  end if


  DO i = 1,nsmp
    CALL AddPRIMEconc( i,xr,sy,xmap,ymap )
  END DO

END IF

umom = 0.0
vmom = 0.0

9999 RETURN

END

SUBROUTINE wakflg(hs,dsbh,dsbw,wake)
!***********************************************************************
!                 WAKFLG Module of ISC2 Short Term Model - ISCST2
! ----------------------------------------------------------------------
! ---    ISC-PRIME     Version 1.0    Level 970812              Modified
! ---        D. Strimaitis
! ---        Earth Tech, Inc.
!            Prepared for EPRI under contract WO3527-01
! ----------------------------------------------------------------------

!        PURPOSE: To Set Wake Flags for Building Downwash Algorithms

!        PROGRAMMER: Roger Brode, Jeff Wang

!        DATE:    March 2, 1992

!        INPUTS:  Building Dimensions
!                 Source Parameters
!                 Meteorological Variables for One Hour

!        OUTPUTS: Logical Flags for Wake Switche, WAKE
!***********************************************************************
IMPLICIT NONE
REAL            :: dsbh,dsbw,hs
LOGICAL         :: wake

!     Set Initial Wake Switches Based on Building Dimensions
IF (dsbh==0.0 .or. dsbw==0.0 .or. &
    hs >= (dsbh + 1.5*AMIN1(dsbh,dsbw))) THEN
    wake   = .false.
ELSE
    wake   = .true.
END IF

RETURN
END

SUBROUTINE wsadj(us,hs,uref,zref,p)
!***********************************************************************
!                 WSADJ Module of ISC2 Short Term Model - ISCST2
!
!        PURPOSE: Adjusts Wind Speed from Anemometer Height to Stack Height
!
!        PROGRAMMER: Roger Brode, Jeff Wang
!
!        INPUTS:  Arrays of Source Parameters
!                 Meteorological Variables for One Hour
!                 Wind Speed Profile Exponents (Default or User-defined)
!
!        OUTPUTS: Stack Top Wind Speed, US
!
!***********************************************************************
IMPLICIT NONE
REAL         :: us,hs,uref,zref,p

!     Adjust Wind Speed -- Assume Wind Speed Constant Below 10 meters
IF (hs >= 10.0) THEN
    us = uref * (hs/zref)**p
ELSE IF (zref > 10.0) THEN
    us = uref * (10.0/zref)**p
ELSE
    us = uref
END IF

!     Do Not Allow Stack Height Wind Speed < 1.0 m/s
IF (us < 1.0) THEN
    us = 1.0
END IF

RETURN
END

SUBROUTINE get_sector(u,v,s,ifvsec,afv)
  IMPLICIT NONE
  ! Input u,v; Output s,ifvsec,afv
  REAL    :: u,v,s,beta,deg2rad,afv
  INTEGER :: ifvsec

  deg2rad = ATAN(1.0)/45.0
  s = SQRT(u*u+v*v)
  IF( u == 0. .AND. v == 0. )THEN
    afv = 0.
  ELSE
    afv = ATAN2(u,v)    ! -pi to +pi from North.
  END IF
  beta = afv/deg2rad
  IF (beta < 0.0) beta = beta + 360.0
  ifvsec = INT(beta*0.1 + 0.4999)
  IF (ifvsec == 0) ifvsec = 36

RETURN
END
!==============================================================================

SUBROUTINE init_stack_rel_prime()

USE scipuff_fi
USE sciprime_fi

IMPLICIT NONE
! Initialise variables

! LDBHR - Debug output written when .TRUE.
ldbhr   = .FALSE.
iwrk2   = 0
isstat  = 0
wake    = .FALSE.
kyprm   = 0.0
kzprm   = 0.0
adsbh   = 0.0
adsbw   = 0.0
adsbl   = 0.0
adsxadj = 0.0
adsyadj = 0.0
isrural = .FALSE.
rural   = .FALSE.
urban   = .TRUE.
srcid   = ' '
! Call NUMMET to set up wind, temp profiles when .TRUE. else use GETMET
primeprof = .false.

keywd(1:ikn) = (/'RURAL   ','BUILDHGT','BUILDWID',&
                  'BUILDLEN','XBADJ   ','YBADJ   '/)

RETURN
END
!==============================================================================

SUBROUTINE set_stack_rel_prime_null()

USE scipuff_fi
USE sciprime_fi

IMPLICIT NONE

cmass_prm = 0.

frac_prm(1) = 1.
frac_prm(2) = 0.

sigy_prm(1:2) = 1.
sigy_prm(1:2) = 1.
sigz_prm(1:2) = 1.
xrel_prm(1:2) = 0.
yrel_prm(1:2) = 0.
zrel_prm(1:2) = 0.
ky_prm(1:2)   = 1.
kz_prm(1:2)   = 1.
wmom_prm(1:2) = 0.
buoy_prm(1:2) = 0.

umom = 0.0
vmom = 0.0

RETURN
END

!==============================================================================

SUBROUTINE load_prime_rel( irel,zbar )

USE scipuff_fi
USE sciprime_fi

IMPLICIT NONE

INTEGER :: irel
REAL    :: zbar

cmass   = cmass_prm*frac_prm(irel)

sigx    = sigy_prm(irel)
sigy    = sigy_prm(irel)
sigz    = sigz_prm(irel)
xrel    = xrel_prm(irel)
yrel    = yrel_prm(irel)
zbar    = zrel_prm(irel)
kyprm   = ky_prm(irel)
kzprm   = kz_prm(irel)
wmom    = wmom_prm(irel)
buoy    = buoy_prm(irel)

RETURN
END

!==============================================================================

SUBROUTINE AddPRIMEconc( ismp,xr,sy,xmap,ymap )

USE scipuff_fi
USE sampler_fi
USE sciprime_fi
USE met_fi

IMPLICIT NONE

INTEGER, INTENT( IN ) :: ismp
REAL,    INTENT( IN ) :: xr, sy
REAL,    INTENT( IN ) :: xmap, ymap

INTEGER ipositn
REAL    xs, ys, x, y, distr, xbrec, ybrec, zflag
REAL    szout,syout,szcav,sycav, conc, xfac

xs = (smp(ismp)%x-xrel)/xmap  !Distance from receptor to source in project coordinates
ys = (smp(ismp)%y-yrel)/ymap

!     Calculate Downwind (X) and Crosswind (Y) Distances
x =  xs*wdsin + ys*wdcos
y = -xs*wdcos + ys*wdsin

IF( x >= xr )GOTO 9999

!     Calculate Source-Receptor (Radial) Distance, DISTR
distr = SQRT( x*x + y*y )

xbrec = x-xadj
ybrec = y-yadj

!DO ii = 1,2

zflag = smp(ismp)%zh

call position( xbrec,ybrec,zflag,ipositn )

IF( ipositn.EQ.4 .AND. ABS(y).GT.1.191754*x )THEN
!              Receptor is at least 50 deg. off the plume centerline
! ---          and is not within a building wake
  GOTO 9999
ELSE IF( ipositn.NE.2 .AND. distr.LT.0.99) then
!              Receptor Too Close to Source for Calculation
! ---          and is not within a building near-wake (cavity)
  GOTO 9999
! -------------------------------------------------------------------

ELSE

  IF( wts < 0. .OR. (zrel_prm(1) <= zi_rel) .OR. (zrel_prm(1) <= zi_rel) )THEN
!                    Determine Simple Terrain Sigmas        ---   CALL PDIS
       !!!CALL PDIS(X,SY,SZ,XY,XZ,SBID)
       call WAKE_XSIG(x,zrel_prm(1)-zrel,.FALSE.,szout,syout,szcav,sycav)
!                 Apply Intermediate Terrain Logic
!                    Simple Terrain Model Only              ---   CALL PSIMPL
       CALL PSIMPL( x,y,zflag,zrel_prm(1),szout,syout,conc   )

       xfac = sy*( 1. + MIN( 2.*(uubl+vvbl)/(MAX(uu_calm,ub**2+vb**2)),1.) )
       xfac = 1. - EXP(-0.5*((x-xr)/(xfac))**2)
       concPrime(ismp) = conc * xfac

   END IF

END IF

9999 CONTINUE

RETURN
END

!==============================================================================

SUBROUTINE PSIMPL( xr,yr,zflagr,hes,szs,sys,hrval )

!USE MAIN1_INC
!USE MAIN2_INC
!USE MAIN3_INC
!***********************************************************************
!               PSIMPL Module of ISC2 Short Term Model - ISCST2
! ----------------------------------------------------------------------
! ---    ISC-PRIME     Version 1.0    Level 970812              Modified
! ---        D. Strimaitis
! ---        Earth Tech, Inc.
!            Prepared for EPRI under contract WO3527-01
! ----------------------------------------------------------------------
!
!        PURPOSE: Calculates Hourly Concentration or Deposition
!                 value for POINT Sources
!                 Using Gaussian Plume Equation for Simple Terrain
!
!                 (Replaces PCHI and PDEP)
!
!           NOTE: Particle settling is treated as a "tilted plume"
!                 until the centerline reaches the surface.  Thereafter
!                 the centroid height of the plume continues to be
!                 modified by gravity.  This process is simulated by
!                 altering the sigma-z for each particle-size.  Hence,
!                 sigma-z is now a function of particle-size.
!
!        PROGRAMMER: D. Strimaitis, SRC
!
!        DATE:    December 15, 1993
!
!        MODIFIED:   To call PDEP for call to SUB. DEPCOR; to use
!                    modified SUB. VERT.  R.W. Brode, PES, Inc. - 9/30/94
!
!        INPUTS:
!
!        OUTPUTS: HRVAL, Concentration or Deposition for Particular
!                 Source/Receptor Combination
!
!        CALLED FROM:   PCALC
!***********************************************************************

USE met_fi
USE sciprime_fi
USE constants_fd

IMPLICIT NONE

!     Variable Declarations

REAL, INTENT( IN  ) :: xr, yr, zflagr, hes, szs, sys
REAL, INTENT( OUT ) :: hrval

REAL, PARAMETER :: EXPLIM = -50.0

! --- Declare local PRIME arrays for "3-source" data
REAL q2(3),y2(3),sy2(3),z2(3),h2(3),sz2(3),qc2(3)

INTEGER is, n1, n2

REAL qtk, fqcav, yterm, a0, v, vterm, emifac
REAL y, sy, sz, he, us, zflag

! --- PRIME ---------------------------------------------------
! --- When there is a building wake, consider treatment of mass in
! --- cavity as additional sources, or as only source
! --- Place selected plume data into transfer arrays (first element)

qtk=cmass_prm

q2(1)=qtk
y2(1)=yr
sy2(1)=sys
z2(1)= zflagr
h2(1)=hes
sz2(1)=szs
n1=1
n2=1

! ---    Define cavity source                              ---   CALL CAV_SRC
CALL CAV_SRC(xr,yr,zflagr,fqcav,qc2,h2,y2,z2,sz2,sy2,n1,n2)
IF( fqcav>0.0 )THEN
! ---       Set source strengths
  q2(1)=qtk*(1.0-fqcav)
  q2(2)=qtk*fqcav*qc2(2)
  q2(3)=qtk*fqcav*qc2(3)
END IF

hrval = 0.0

! --- Loop over 3 possible sources
DO is=n1,n2

! --- Transfer data for current source
qtk=q2(is)
y=y2(is)
sy=sy2(is)
zflag=z2(is)
he=h2(is)
sz=sz2(is)
! ----------------------------------------------------------------------

emifac = 1. !or 3600?

us = ustack !Saved from above

YTERM = -0.5*(Y*Y)/(SY*SY)
IF (YTERM > EXPLIM) THEN

      V = 0.
!           Calculate the Vertical Term, V, for gases
!                 Calculate Concentration Form of V         ---   CALL VERT
      A0 = -0.5/(SZ*SZ)
      CALL VERT(HE,SZ,A0,ZFLAG,V)
!           Include SZ in the denomenator of V
      V = V/SZ


!           Complete VTERM (SZ already in denomenator of V)
      VTERM = V/(PI2*US*SY)

!           Check for Possible Underflow Condition
      IF (VTERM>0.0 .AND. (LOG(VTERM)+YTERM)>EXPLIM) THEN
         HRVAL = hrval + &
                 QTK * emifac * VTERM * EXP(YTERM)
      END IF

END IF

END DO

 999  RETURN
END

!==============================================================================

SUBROUTINE VERT(HEARG,SZARG,A0,ZARG,VOUT)

!USE MAIN1_INC
!USE MAIN2_INC
!USE MAIN3_INC
!***********************************************************************
!                 VERT Module of ISC2 Short Term Model - ISCST2
!
!        PURPOSE: Calculates Vertical Term for Use in Gaussian Plume Equation
!
!        PROGRAMMER: Roger Brode, Jeff Wang
!
!        MODIFIED BY R.W. Brode, PES, Inc. to use calling arguments - 9/30/94
!
!        MODIFIED BY D. Strimaitis, SRC (for Wet REMOVAL of Gases)
!
!        DATE:    November 8, 1993
!
!
!        INPUTS:  Plume Height
!                 Vertical Dispersion Parameter
!                 Stability Class
!                 Mixing Height
!                 Receptor Height Above Ground
!
!        OUTPUTS: Vertical Term, VOUT
!
!        CALLED FROM:   PSIMPL, PCOMPL, ASIMPL
!***********************************************************************

USE met_fi
USE sciprime_fi
USE constants_fd

IMPLICIT NONE

REAL, INTENT( IN  ) :: HEARG,SZARG,A0,ZARG
REAL, INTENT( OUT ) :: VOUT

REAL, PARAMETER :: EXPLIM = -50.0

!     Variable Declarations
INTEGER i
REAL    a1, a2, a3, a4, a5, a6, sum, twoizi, t, SRT2PI

!     Variable Initializations

VOUT = 0.0

SRT2PI = SQRT(PI2)

IF (ZARG == 0.0) THEN
!        Vertical Term for Case With No Flagpole Receptor
!   IF (STABLE .OR. ZI>=10000.) THEN
   IF (wts < 0. .OR. zi_rel >= 10000.) THEN
      A1 = A0 * HEARG * HEARG
      IF (A1 > EXPLIM)  VOUT = 2.*EXP(A1)
   ELSE IF ((SZARG/zi_rel) >= 1.6) THEN
      VOUT  = SRT2PI*(SZARG/zi_rel)
   ELSE
      A1 = A0 * HEARG * HEARG
      IF (A1 > EXPLIM)  VOUT = EXP(A1)
      SUM = 0.0
      DO 100 I = 1, 100
         T  = 0.0
         TWOIZI = 2.*FLOAT(I)*zi_rel
         A2 = A0 * (TWOIZI-HEARG) * (TWOIZI-HEARG)
         A3 = A0 * (TWOIZI+HEARG) * (TWOIZI+HEARG)
         IF (A2 > EXPLIM)  T = EXP(A2)
         IF (A3 > EXPLIM)  T = T + EXP(A3)
         SUM = SUM + T
         IF (ABS(T) <= 5.0E-9) THEN
!                 Exit Loop
            GOTO 200
         END IF
 100        CONTINUE
!           Calculate Total Vert. Term - (2.*) was Removed for Optimization
 200        VOUT  = 2.*(VOUT + SUM)
   END IF
ELSE
!        Vertical Term for Case of ZARG .NE. 0.0
   IF (wts < 0. .OR. zi_rel >= 10000.) THEN
      A1 = A0 * (ZARG-HEARG) * (ZARG-HEARG)
      A2 = A0 * (ZARG+HEARG) * (ZARG+HEARG)
      IF (A1 > EXPLIM)  VOUT = EXP(A1)
      IF (A2 > EXPLIM)  VOUT = VOUT + EXP(A2)
   ELSE IF (SZARG/zi_rel >= 1.6) THEN
      VOUT  = SRT2PI*(SZARG/zi_rel)
   ELSE
      A1 = A0 * (ZARG-HEARG) * (ZARG-HEARG)
      A2 = A0 * (ZARG+HEARG) * (ZARG+HEARG)
      IF (A1 > EXPLIM)  VOUT = EXP(A1)
      IF (A2 > EXPLIM)  VOUT = VOUT + EXP(A2)
      SUM = 0.0
      DO 300 I = 1, 100
         T  = 0.0
         TWOIZI = 2.*FLOAT(I)*zi_rel
         A3 = A0 * (ZARG-(TWOIZI-HEARG)) * (ZARG-(TWOIZI-HEARG))
         A4 = A0 * (ZARG+(TWOIZI-HEARG)) * (ZARG+(TWOIZI-HEARG))
         A5 = A0 * (ZARG-(TWOIZI+HEARG)) * (ZARG-(TWOIZI+HEARG))
         A6 = A0 * (ZARG+(TWOIZI+HEARG)) * (ZARG+(TWOIZI+HEARG))
         IF (A3 > EXPLIM)  T = T + EXP(A3)
         IF (A4 > EXPLIM)  T = T + EXP(A4)
         IF (A5 > EXPLIM)  T = T + EXP(A5)
         IF (A6 > EXPLIM)  T = T + EXP(A6)
         SUM = SUM + T
         IF (ABS(T) <= 1.0E-8) THEN
!                 Exit Loop
            GOTO 400
         END IF
 300        CONTINUE
 400        VOUT  = VOUT + SUM
   END IF
END IF

RETURN
END
