subroutine drvtuv(igrd,ncol,nrow,nlay,nspc,time,date,itzon, & idfin,cellat,cellon,height,press,tsurf, & tempk,water,conc,cod,alb,cldtrns) use o3colmap use chmstry c c----CAMx v7Beta6 190902 c c DRVTUV is the driver for the cloud/aerosol adjustment to c clear-sky photolysis rates. A streamlined version of the TUV c radiative transfer model (based on v4.8) is called to calculate c actinic fluxes through the grid column for 2 cases: clear with c default haze, and cloudy with actual haze. The ratio c of the two actinic flux profiles is used to apply to c the clear-sky lookup photolysis rates in KPHOTO.F. c c Copyright 1996 - 2018 c Ramboll c c Modifications: c 03/29/11 Support in-line TUV with aerosol optical depth c 04/02/12 Removed RADM cloud adjustment option, cloud/aerosol c adjustments now always done with in-line TUV c 04/20/13 Moved call to EMISTRNS, removed averaging across c input times c 08/11/14 Albedo now includes effects of snow c c Input arguments: c igrd grid index c ncol number of columns c nrow number of rows c nlay number of layers c nspc number of species c time current time (HHMM) c date current date (YYJJJ) c itzon time zone c idfin map of nested grids in this grid c cellat cell centroid latitude (deg) c cellon cell centroid longitude (deg) c height current layer interface height field (m) c press current layer pressure field (mb) c tsurf current surface temperature field (K) c tempk current layer temperature field (K) c water current layer water vapor field (mb) c conc concentration field (umol/m3, ug/m3) c cod layer cloud optical depth field (unitless) c alb surface UV albedo c c Output arguments: c cldtrns Cloud adjustment factor field (unitless) c c Routines called: c GETZNTH c TUV c c Called by: c EMISTRNS c implicit none include 'camx.prm' include 'camx_aero.inc' include 'flags.inc' c c-----Arguments c integer igrd, ncol, nrow, nlay, nspc, date, itzon integer idfin(ncol,nrow) real time real cellat(ncol,nrow), cellon(ncol,nrow), alb(ncol,nrow) real height(ncol,nrow,nlay) real press(ncol,nrow,nlay) real tsurf(ncol,nrow) real tempk(ncol,nrow,nlay) real water(ncol,nrow,nlay) real conc(ncol,nrow,nlay,nspc) real cod(ncol,nrow,nlay) real cldtrns(ncol,nrow,nlay) c c-----Local variables c integer NELT parameter (NELT = 51) integer i, j, k, l integer iz, izz, nz, nn integer rh1d(MXLAYER) real tuvalb real prslev,tmplev real zen1, zen2, zen, coszen real dtdz, tavg real wvapor,qwatr,ev,es,rhfac real totext,totssa,totcon,sumcod real factr,odint real temp1d(MXLAYER),pres1d(MXLAYER) real rafcld(MXLAYER+1) real z(MXLAYER+1),midht1d(MXLAYER) real airlev(MXLAYER+1) real odcld(MXLAYER+1), omcld(MXLAYER+1), gcld(MXLAYER+1) real odaer1(MXLAYER+1), omaer1(MXLAYER+1), gaer(MXLAYER+1) real odaer2(MXLAYER+1), omaer2(MXLAYER+1) real frh(100) real eltod(NELT) logical ldark real RGAS parameter(RGAS = 1.3806503d-19) !cm3 mbar K-1 molec-1 c c-----Aerosol optical depth profile from Elterman (1968). c These are vertical optical depths per km, in 1 km c intervals from 0 km to 50 km, at 340 nm. c Total integrated OD is 0.38 (~0.23 at 550 nm). c data eltod/ & 2.40E-01,1.06E-01,4.56E-02,1.91E-02,1.01E-02,7.63E-03, & 5.38E-03,5.00E-03,5.15E-03,4.94E-03,4.82E-03,4.51E-03, & 4.74E-03,4.37E-03,4.28E-03,4.03E-03,3.83E-03,3.78E-03, & 3.88E-03,3.08E-03,2.26E-03,1.64E-03,1.23E-03,9.45E-04, & 7.49E-04,6.30E-04,5.50E-04,4.21E-04,3.22E-04,2.48E-04, & 1.90E-04,1.45E-04,1.11E-04,8.51E-05,6.52E-05,5.00E-05, & 3.83E-05,2.93E-05,2.25E-05,1.72E-05,1.32E-05,1.01E-05, & 7.72E-06,5.91E-06,4.53E-06,3.46E-06,2.66E-06,2.04E-06, & 1.56E-06,1.19E-06,9.14E-07/ c c-----Aerosol humidity adjustment parameter as f(RH); FLAG (2000) c data frh/1.,1.,1.,1.,1., & 1.,1.,1.,1.,1., & 1.,1.,1.,1.0001,1.0001, & 1.0004,1.0006,1.0024,1.0056,1.0089, & 1.0097,1.0105,1.0111,1.0115,1.0118, & 1.0122,1.0126,1.0130,1.0135,1.0139, & 1.0173,1.0206,1.0254,1.0315,1.0377, & 1.0486,1.0596,1.0751,1.0951,1.1151, & 1.1247,1.1343,1.1436,1.1525,1.1615, & 1.1724,1.1833,1.1955,1.2090,1.2224, & 1.2368,1.2512,1.2671,1.2844,1.3018, & 1.3234,1.3450,1.3695,1.3969,1.4243, & 1.4628,1.5014,1.5468,1.5992,1.6516, & 1.6991,1.7466,1.7985,1.8549,1.9113, & 1.9596,2.0080,2.0596,2.1146,2.1695, & 2.2630,2.3565,2.4692,2.6011,2.7330, & 2.8461,2.9592,3.0853,3.2245,3.3637, & 3.5743,3.7849,4.0466,4.3594,4.6721, & 5.3067,5.9412,6.9627,8.3710,9.7793, & 12.4288,15.0773,18.0590,21.3709,22. / c c-----Entry point c nz = nlay + 1 c c-----Loop over all vertical grid columns; skip columns containing nested grids c c$omp parallel default(shared) c$omp& private(i,j,k,rafcld,zen,ldark,coszen, c$omp& sumcod,tuvalb,z,temp1d,pres1d,wvapor, c$omp& qwatr,ev,es,rh1d,midht1d,airlev,prslev, c$omp& tmplev,dtdz,tavg,iz,odcld,omcld,gcld, c$omp& izz,factr,odint,odaer1,omaer1,gaer, c$omp& odaer2,omaer2,totext,totssa,totcon,l,rhfac) c$omp do schedule(static) collapse(2) do j = 2,nrow-1 do i = 2,ncol-1 do k = 1,nz rafcld(k) = 1. enddo if (idfin(i,j).gt.igrd) goto 100 c c-----Get solar zenith angle c call getznth(cellat(i,j),cellon(i,j),time,date,itzon,zen, & ldark) if (ldark) goto 100 zen = min(zen,75.) coszen = COS(zen*3.1415/180.) c c-----Check if this is a cloud-free column; adjustmnet ratio = 1 c sumcod = 0. do k = 1,nlay sumcod = sumcod + cod(i,j,k) enddo if (sumcod.eq.0. .AND. naero.eq.0) goto 100 c c-----Set albedo c tuvalb = alb(i,j) c c-----Prepare 1-D met variables for TUV c z(1) = 0. do k = 1,nlay z(k+1) = height(i,j,k)*1.e-3 temp1d(k) = tempk(i,j,k) pres1d(k) = press(i,j,k) wvapor = water(i,j,k) qwatr = 1.e-6*wvapor*18./28.8 ev = qwatr*pres1d(k)/(qwatr + eps) es = e0*exp((lv/rv)*(1./273. - 1./temp1d(k))) rh1d(k) = nint(100.*max(0.01,min(0.95,ev/es))) enddo do k = 1,nlay midht1d(k) = (z(k) + z(k+1))*0.5 enddo c c-----Interpolate temperature and pressure to layer interfaces, then c calculate air density c prslev = pres1d(1) - (midht1d(1)-z(1))* & (pres1d(2)-pres1d(1))/(midht1d(2)-midht1d(1)) airlev(1) = prslev/(Rgas*tsurf(i,j)) do iz = 2,nz-1 prslev = pres1d(iz-1) + (z(iz)-midht1d(iz-1))* & (pres1d(iz)-pres1d(iz-1))/(midht1d(iz)-midht1d(iz-1)) tmplev = temp1d(iz-1) + (z(iz)-midht1d(iz-1))* & (temp1d(iz)-temp1d(iz-1))/(midht1d(iz)-midht1d(iz-1)) airlev(iz) = prslev/(Rgas*tmplev) enddo dtdz = (temp1d(nz-1) - temp1d(nz-2))/ & (midht1d(nz-1) - midht1d(nz-2)) tmplev = temp1d(nz-1) + dtdz*(z(nz) - midht1d(nz-1)) tavg = (tmplev + temp1d(nz-1))/2. prslev = pres1d(nz-1)* & exp(-9.8*(z(nz) - midht1d(nz-1))/(2.*287.*tavg)) airlev(nz) = prslev/(Rgas*tmplev) c c-----Set cloud params c do iz = 1,nz-1 odcld(iz) = cod(i,j,iz) omcld(iz) = 0.9999 gcld(iz) = 0.85 enddo c c-----Set aerosol params: all are assumed to be applicable at 340 nm. c For the default case, use the Elterman profile. c If aerosols are run, use PM concentration fields and optical c params provided by the chemparam file. c izz = 1 do 30 k = 1,NELT-1 do iz = izz,nz-1 if (midht1d(iz) .lt. float(k)) then factr = midht1d(iz) - float(k-1) odint = factr*eltod(k+1) + (1. - factr)*eltod(k) odaer1(iz) = odint*(z(iz+1) - z(iz)) omaer1(iz) = 0.90 gaer(iz) = 0.61 else izz = iz goto 30 endif enddo goto 31 30 continue c 31 do iz = 1,nz-1 if (naero.eq.0) then odaer2(iz) = odaer1(iz) omaer2(iz) = omaer1(iz) else totext = 0. totssa = 0. totcon = 0. do l = ngas+1,nspc rhfac = 1. if (rhadj(l).eq.1) rhfac = frh(rh1d(iz)) totext = totext + bext(l)*rhfac*conc(i,j,iz,l) totssa = totssa + ssa(l)*conc(i,j,iz,l) totcon = totcon + conc(i,j,iz,l) enddo odaer2(iz) = totext*(z(iz+1)-z(iz))*1000. omaer2(iz) = max(0.20,min(0.99,totssa/totcon)) endif enddo c c-----Call the tuv-cloud routine c call tuv(nz,z,airlev,tuvalb,coszen,odcld,omcld,gcld, & odaer1,odaer2,omaer1,omaer2,gaer,rafcld) c c-----Load 3-D CLDTRNS array c 100 do k = 1,nlay cldtrns(i,j,k) = 0.5 * (rafcld(k) + rafcld(k+1)) enddo enddo enddo c$omp end parallel c c-----Set boundary to be clear sky c do j = 1,nrow do k = 1,nlay cldtrns(1,j,k) = 1. cldtrns(ncol,j,k) = 1. enddo enddo do i = 1,ncol do k = 1,nlay cldtrns(i,1,k) = 1. cldtrns(i,nrow,k) = 1. enddo enddo return end