subroutine zadvec(m1,m2,m3,i0,j0,ia,iz,ja,jz,ibcon, & losat,ltopcon,igrd,ncol,nrow,nlay,nrads,nspc, & nsen,delt,dx, & dy,idfin,depth,entrn,dilut,tempk,press,species, & conc,rhoscl,ctop,lvupsolv,ptrtop,fluxes,fluxtmp, & isaptr,ipa_xy,ipa_lay,iproc_id) use filunit use chmstry use bndary use procan use rtracchm use tracer c c----CAMx v7Beta6 190902 c c ZADVEC drives vertical transport of concentrations. The operation c is performed on the regular model species (losat = False) or source c apportionment concentrations (losat = True). c This version also performs vertical transport of sensitivities c if DDM is enabled. c c Copyright 1996 - 2018 c Ramboll c c Modifications: c 12/3/99 Top boundary condition is converted from ppm to umol/m3 c using extrapolated density in an overlying layer c that is assumed to be the same thickness as the top c layer of the model c 12/07/01 added instructions for OMP c 01/30/02 Added code for RTRAC probing tool c 3/06/06 Revised top BC approach; removed top dummy layer c 05/01/07 Fixed bug that was causing PM PSAT species to c be converted to micro-moles c 07/16/07 -bkoo- Added check for HDDM c 07/16/08 -bkoo- Added DDM turn-off flag c 04/09/09 Improved definition of atmospheric density at top of model c 11/04/09 Revised vertical advection solver technique, now c employs zero-gradient top boundary condition c 11/12/09 Added code to save fluxes and apply to the tracers c 4/07/14 Added top con file c 11/28/16 -bkoo- Updated DDM for new top con c c Input arguments: c ltopcon top con flag c igrd grid index c ncol number of columns c nrow number of rows c nlay number of layers c nrads number of radicals c nspc number of species c nsen number of species X number of DDM parameters c or number of tracer species c or 1, whichever is larger c delt time step (s) c dx cell size in x-direction (m) c dy cell size in y-direction (m) c idfin map of nested grids in this grid c depth layer depth (m) c entrn entrainment rate (m/s) c dilut dilution rate (m/s) c tempk temperature (K) c press pressure (mb) c species species names c conc species concentrations (umol/m3) c rhoscl density scaling above top of model (unitless) c ctop initial concentrations in top layer (umol/m3) c lvupsolv 3-D flag to use upstream donor solver c ptrtop initial sensitivities in top layer c fluxtmp temporary array for fluxes c ipa_xy 2-D gridded array to identify if cell is c in a IPRM sub-domain c ipa_lay 3-D gridded array to identify IPRM sub-domain c each layer is in c iproc_id ID for this processor in MPI rank c c Output arguments: c conc species concentrations (umol/m3) c fluxes fluxes across the boundaries (umol) c c Routines Called: c VRTSLV c c Called by: c EMISTRNS c implicit none include "camx.prm" include "camx.inc" c integer :: m1,m2,m3,i0,j0,ia,iz,ja,jz,ibcon integer :: m_zi1,m_zi2,m_zj1,m_zj2 c integer ncol integer nrow integer nlay integer nrads integer nspc integer igrd integer nsen character*10 species(nspc) real conc(m1,m2,m3,nspc) real fluxtmp(m1,m2,m3,nspc) real ctop(m1,m2,nspc) real rhoscl(m1,m2) logical lvupsolv(m1,m2,m3) logical ltopcon real ptrtop(m1,m2,nsen) integer*8 isaptr real, dimension(m1,m2,m3) :: entrn,dilut,tempk,press real, dimension(m1,m2,m3) :: depth real dx(nrow) integer idfin(m1,m2) logical losat real*8 fluxes(nspc,11) real fluxtop real delt real dy integer iproc_id integer ispc,j,i,k,numspcs real convfac real c1d(MXLAYER+1) real d1d(MXLAYER) real ent1d(MXLAYER) real dil1d(MXLAYER) logical lvus(MXLAYER) c c===========================DDM Begin================================= c integer ls,isen,ioff integer*8 idx real sen1d((MXLAYER+1)*MXFDDM) c c===========================DDM End=================================== c c=================== Process Analysis Begin ========================== c integer ipa_xy(ncol,nrow) integer ipa_lay(ncol,nrow,nlay) integer ipa_idx logical ldoipts real fc1(MXLAYER) real fc2(MXLAYER) real fc3(MXLAYER) c c==================== Process Analysis End =========================== c c=================== Source Apportionment Begin ===================== c real fluxlay(MXLAYER) c c=================== Source Apportionment End ======================= c c-----Entry point c ! Set grid point limits m_zi1= ia-1 m_zi2= iz+1 m_zj1= ja-1 m_zj2= jz+1 if(btest(ibcon,0)) m_zi1 = ia if(btest(ibcon,1)) m_zi2 = iz if(btest(ibcon,2)) m_zj1 = ja if(btest(ibcon,3)) m_zj2 = jz c c ---- Initialize the temp flux array to zero --- c if( .NOT. losat ) call zeros(fluxtmp,m1*m2*m3*nspc) c c-----Vertical advection, entrainment, and dilution c numspcs = nspc if( losat ) numspcs = nsaspc if( iproc_id .LE. 1 ) then if( .NOT. losat ) then write(*,'(a20,$)') 'z advection ......' else write(*,'(a20,$)') 'SA z advection ...' endif endif if( .NOT. losat ) then write(iout,'(a20,$)') 'z advection ......' else write(iout,'(a20,$)') 'SA z advection ...' endif c c$omp parallel default(shared) c$omp& private(ispc,i,j,k,ent1d,dil1d,d1d,c1d,lvus,fluxlay, c$omp& convfac,ls,isen,ioff,idx,sen1d,fluxtop, c$omp& ldoipts,ipa_idx,fc1, c$omp& fc2,fc3) c c$omp do schedule(guided) c do 40 ispc = nrads+1,numspcs do 60 j = m_zj1, m_zj2 !2,nrow-1 do 50 i = m_zi1, m_zi2 !2,ncol-1 c c-----Skip cells occupied by child grids c if (idfin(i,j).gt.igrd) goto 50 do k = 1,nlay ent1d(k) = entrn(i,j,k) dil1d(k) = dilut(i,j,k) d1d(k) = depth(i,j,k) c1d(k) = conc(i,j,k,ispc) lvus(k) = lvupsolv(i,j,k) enddo c1d(nlay+1) = rhoscl(i,j)*ctop(i,j,ispc) if (ltopcon) then convfac = 1. if (ispc .LE. ngas) & convfac = densfac*(273./tempk(i,j,nlay))* & (press(i,j,nlay)/1013.) c1d(nlay+1) = convfac*c1d(nlay+1) endif c c================================DDM Begin============================ c if( (lddm .OR. lhddm) .AND. lddmcalc(igrd) ) then ls = 0 do isen = 1,nddmsp ioff = iptddm(ispc)+isen-1 do k = 1,nlay ls = ls + 1 idx = DBLE(isaptr-1) + DBLE(i) + & DBLE(m1)*DBLE(j-1) + DBLE(m1)*DBLE(m2)*DBLE(k-1) + & DBLE(m1)*DBLE(m2)*DBLE(m3)*DBLE(ioff-1) sen1d(ls) = ptconc(idx) enddo ls = ls + 1 sen1d(ls) = rhoscl(i,j)*ptrtop(i,j,ioff) if (ltopcon) then convfac = 1. if (ispc .LE. ngas) & convfac = densfac*(273./tempk(i,j,nlay))* & (press(i,j,nlay)/1013.) sen1d(ls) = convfac*sen1d(ls) endif enddo endif c c================================DDM End============================== c c c======================== Process Analysis Begin ===================== c ldoipts = .FALSE. if ( .NOT. ltrace .AND. lipr ) then if( i .GE. ia .AND. i .LE. iz .AND. j .GE. ja & .AND. j .LE. jz ) then if( ipa_xy(i+i0,j+j0) .GT. 0 ) ldoipts = .TRUE. endif endif c c c========================= Process Analysis End ====================== c c-----Solve vertical mass adjustments using Crank-Nicholson solver c call vrtslv(nlay,i,j,igrd,delt,ent1d,dil1d,d1d,c1d, & lvus,fluxtop,fluxlay,sen1d, & species(ispc),fc1,fc2,fc3,ldoipts) c c======================== Process Analysis Begin ===================== c if( ldoipts ) then do k=1,nlay if( ipa_lay(i+i0,j+j0,k) .GT. 0 ) then ipa_idx = ipa_lay(i+i0,j+j0,k) c c-----Concentration change from vertical advection c cipr(IPR_BADV, ipa_idx, ispc) = & cipr(IPR_BADV, ipa_idx, ispc) + fc1(k) c cipr(IPR_TADV, ipa_idx, ispc) = & cipr(IPR_TADV, ipa_idx, ispc) + fc2(k) c cipr(IPR_DADV, ipa_idx, ispc) = 0. c endif enddo endif c c========================= Process Analysis End ====================== c do k = 1,nlay conc(i,j,k,ispc) = c1d(k) enddo if( .NOT. losat ) then if(i .GE. ia .AND. i .LE. iz .AND. j .GE. ja & .AND. j .LE. jz ) then if (fluxtop.lt.0) then fluxtmp(i,j,1,ispc) = -fluxtop*dx(j+j0)*dy*delt else fluxtmp(i,j,2,ispc) = -fluxtop*dx(j+j0)*dy*delt endif endif endif c c================================DDM Begin============================ c if( (lddm .OR. lhddm) .AND. lddmcalc(igrd) ) then ls = -1 do isen = 1, nddmsp ls = ls + 1 ioff = iptddm(ispc)+isen-1 do k = 1, nlay ls = ls + 1 idx = DBLE(isaptr-1) + DBLE(i) + & DBLE(m1)*DBLE(j-1) + DBLE(m1)*DBLE(m2)*DBLE(k-1) + & DBLE(m1)*DBLE(m2)*DBLE(m3)*DBLE(ioff-1) ptconc(idx) = sen1d(ls) enddo enddo endif c c================================DDM End============================== c 50 continue 60 continue 40 continue c c$omp end parallel c c-----Put fluxes in global array c if( .NOT. losat ) then do j = 2,m2-1 do i = 2,m1-1 do ispc = 1,nspc fluxes(ispc,9) = fluxes(ispc,9) + DBLE(fluxtmp(i,j,1,ispc)) fluxes(ispc,10) = fluxes(ispc,10) + DBLE(fluxtmp(i,j,2,ispc)) enddo enddo enddo endif c if( iproc_id .LE. 1 ) then write(*,'(a)') ' Done' call flush(6) endif write(iout,'(a)') ' Done' call flush(iout) c return end