#define MODAL_AERO #define WRF_PORT ! module_cam_mam_wateruptake.F ! adapted from cam3 modal_aero_wateruptake.F90 by r.c.easter, june 2010 !Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov ! 2010-06-24 rce notes !-------------------------------------------------------------- ! modal_aero_wateruptake.F90 !---------------------------------------------------------------------- !BOP ! ! !MODULE: modal_aero_wateruptake --- modal aerosol mode merging (renaming) ! ! !INTERFACE: module modal_aero_wateruptake ! !USES: use shr_kind_mod, only : r8 => shr_kind_r8 #ifndef WRF_PORT use cam_logfile, only: iulog #else use module_cam_support, only: iulog #endif implicit none private save ! !PUBLIC MEMBER FUNCTIONS: public modal_aero_wateruptake_sub ! !PUBLIC DATA MEMBERS: ! currently none ! !DESCRIPTION: This module implements ... ! ! !REVISION HISTORY: ! ! RCE 07.04.13: Adapted from MIRAGE2 code ! !EOP !---------------------------------------------------------------------- !BOC ! list private module data here !EOC !---------------------------------------------------------------------- contains !---------------------------------------------------------------------- subroutine modal_aero_wateruptake_sub( & lchnk, ncol, nstep, & iwaterup_flag, loffset, & aero_mmr_flag, h2o_mmr_flag, & deltat, h2ommr, t, pmid, pdel, cldn, & raer, raertend, dotend, qaerwat, & dgncur_a, dgncur_awet, wetdens & #ifdef OBSRH , obsrh & #endif ) ! following are local for now -- wait and see ! maer, naer, wetrad, density ) !----------------------------------------------------------------------- ! ! Purpose: Compute aerosol wet radius ! ! Method: Kohler theory ! ! Author: S. Ghan ! !----------------------------------------------------------------------- use modal_aero_data #ifndef WRF_PORT use mo_constants, only: pi use pmgrid, only: plat, plon use ppgrid, only: begchunk, endchunk, pcols, pver #endif use physconst, only: cpair, epsilo, gravit, mwdry, mwh2o, & rair, rga, rhoh2o, rh2o, latvap #ifndef WRF_PORT use constituents, only: pcnst, qmin #endif use wv_saturation, only: aqsat, qsat_water #ifndef WRF_PORT use phys_grid, only: get_lat_all_p, get_lon_all_p use abortutils, only : endrun use cam_history, only : outfld use spmd_utils, only : masterproc #else use module_cam_support, only: pcnst => pcnst_runtime, & pcols, pver, & endrun, masterproc, outfld use physconst, only: pi #endif implicit none integer, intent(in) :: lchnk ! chunk index integer, intent(in) :: ncol ! number of columns integer, intent(in) :: nstep ! time step integer, intent(in) :: iwaterup_flag ! identifies call location from typhysbc (1,2) integer, intent(in) :: loffset ! offset applied to modal aero "pointers" logical, intent(in) :: aero_mmr_flag ! if .true., aerosol q are kg/kg-air ! if .false., aerosol q are mol/mol-air logical, intent(in) :: h2o_mmr_flag ! if .true., h2ommr is kg/kg-air logical, intent(inout)::dotend(pcnst) ! identifies species for which tendencies are computed real(r8), intent(in) :: deltat ! time step (s) real(r8), intent(in) :: h2ommr(pcols,pver) ! layer specific humidity real(r8), intent(in) :: t(pcols,pver) ! layer temperatures (K) real(r8), intent(in) :: pmid(pcols,pver) ! layer pressure (Pa) real(r8), intent(in) :: pdel(pcols,pver) ! layer pressure thickness (Pa) real(r8), intent(in) :: cldn(pcols,pver) ! layer cloud fraction (0-1) real(r8), intent(in) :: raer(pcols,pver,pcnst) ! aerosol species MRs (kg/kg and #/kg) real(r8), intent(inout)::raertend(pcols,pver,pcnst) ! aerosol MR tendencies (kg/kg/s) ! only defined for aerosol water real(r8), intent(out) :: qaerwat(pcols,pver,ntot_amode) real(r8), intent(in) :: dgncur_a(pcols,pver,ntot_amode) real(r8), intent(out) :: dgncur_awet(pcols,pver,ntot_amode) real(r8), intent(out) :: wetdens(pcols,pver,ntot_amode) ! following are local for now -- wait and see ! real(r8), intent(out) :: maer(pcols,pver,ntot_amode) ! ! aerosol wet mass MR (including water) (kg/kg-air) ! real(r8), intent(out) :: naer(pcols,pver,ntot_amode) ! ! aerosol number MR (bounded!) (#/kg-air) ! real(r8), intent(out) :: wetrad(pcols,pver,ntot_amode) ! ! wet radius of aerosol (m) ! real(r8), intent(out) :: density(pcols,pver,ntot_amode) ! ! wet mean density of aerosol (kg/m3) ! local variables integer i,k,m integer icol_diag integer l ! species index integer lmass ! pointer for aerosol mass integer ltype ! pointer for aerosol type ! integer lwater ! pointer for water on aerosol integer lat(pcols), lon(pcols) ! lat,lon indices real(r8) density_water ! density of water (kg/m3) real(r8) drydens(ntot_amode) ! dry particle density (kg/m^3) real(r8) drymass(ntot_amode) ! single-particle-mean dry mass (kg) real(r8) dryrad(pcols,pver,ntot_amode) ! dry volume mean radius of aerosol (m) real(r8) dryvol(ntot_amode) ! single-particle-mean dry volume (m3) real(r8) dryvolmr(ntot_amode) ! volume MR for aerosol mode (m3/kg) real(r8) duma, dumb real(r8) es(pcols,pver) ! saturation vapor pressure (Pa) real(r8) hygro(ntot_amode) ! volume-weighted mean hygroscopicity (--) real(r8) hystfac(ntot_amode) ! working variable for hysteresis real(r8) pi43 real(r8) qs(pcols,pver) ! saturation specific humidity real(r8) qwater ! aerosol water MR real(r8) rh(pcols,pver) ! relative humidity (0-1) real(r8) third real(r8) v2ncur_a(pcols,pver,ntot_amode) real(r8) wtrvol(ntot_amode) ! single-particle-mean water volume in wet aerosol (m3) real(r8) wetvol(ntot_amode) ! single-particle-mean wet volume (m3) real(r8) :: maer(pcols,pver,ntot_amode) ! aerosol wet mass MR (including water) (kg/kg-air) real(r8) :: naer(pcols,pver,ntot_amode) ! aerosol number MR (bounded!) (#/kg-air) real(r8) :: wetrad(pcols,pver,ntot_amode) ! wet radius of aerosol (m) character(len=3) :: trnum ! used to hold mode number (as characters) !----------------------------------------------------------------------- ! set the dotend's ! currently, this is one of several routines called from aerosol_wet_intr ! so DO NOT re-initialize dotend=.false. and raertend=0.0 ! dotend(:) = .false. do m=1,ntot_amode ! lwater = lwaterptr_amode(m) - loffset ! dotend(lwater) = .true. ! raertend(1:ncol,:,lwater) = 0.0 end do third=1./3. pi43 = pi*4.0/3.0 density_water = rhoh2o ! is (kg/m3) ! compute size-related factors ! when numptr_amode(m)=0, earlier code set dryrad = prescribed volume. ! this is no longer needed because ! dryrad(i,k,m) = (1.0/v2ncur_a(i,k,m)/pi43)**third ! works in all cases do m=1,ntot_amode hystfac(m) = 1.0 / max( 1.0e-5_r8, & (rhdeliques_amode(m)-rhcrystal_amode(m)) ) enddo ! main loops over i, k ! call aqsat( t, pmid, es, qs, pcols, ncol, pver, 1, pver ) do k=1,pver do i=1,ncol qs(i,k)=qsat_water(t(i,k),pmid(i,k)) if ( h2o_mmr_flag ) then rh(i,k) = h2ommr(i,k)/qs(i,k) else rh(i,k) = h2ommr(i,k)*mwh2o/(mwdry*qs(i,k)) end if rh(i,k) = max(rh(i,k),0.0_r8) rh(i,k) = min(rh(i,k),0.98_r8) if (cldn(i,k) .lt. 1.0_r8) then rh(i,k) = (rh(i,k) - cldn(i,k)) / (1.0_r8 - cldn(i,k)) ! clear portion end if rh(i,k) = max(rh(i,k),0.0_r8) ! compute dryvolmr, maer, naer for each mode do m=1,ntot_amode maer(i,k,m)=0. dryvolmr(m)=0. hygro(m)=0. do l = 1, nspec_amode(m) lmass = lmassptr_amode(l,m) - loffset ltype = lspectype_amode(l,m) if ( aero_mmr_flag ) then duma = raer(i,k,lmass) else duma = raer(i,k,lmass)*(specmw_amode(ltype)/mwdry) end if maer(i,k,m) = maer(i,k,m) + duma dumb = duma/specdens_amode(ltype) dryvolmr(m) = dryvolmr(m) + dumb hygro(m) = hygro(m) + dumb*spechygro(ltype) enddo if (dryvolmr(m) > 1.0e-30_r8) then hygro(m) = hygro(m)/dryvolmr(m) else hygro(m) = spechygro( lspectype_amode(1,m) ) end if ! naer = aerosol number (#/kg) ! the new v2ncur_a replaces old coding here v2ncur_a(i,k,m) = 1. / ( (pi/6.)* & (dgncur_a(i,k,m)**3.)*exp(4.5*alnsg_amode(m)**2.) ) naer(i,k,m) = dryvolmr(m)*v2ncur_a(i,k,m) enddo !m=1,ntot_amode ! compute mean (1 particle) dry volume and mass for each mode ! old coding is replaced because the new (1/v2ncur_a) is equal to ! the mean particle volume ! also moletomass forces maer >= 1.0e-30, so (maer/dryvolmr) ! should never cause problems (but check for maer < 1.0e-31 anyway) do m=1,ntot_amode if (maer(i,k,m) .gt. 1.0e-31) then drydens(m) = maer(i,k,m)/dryvolmr(m) else drydens(m) = 1.0 end if dryvol(m) = 1.0/v2ncur_a(i,k,m) drymass(m) = drydens(m)*dryvol(m) dryrad(i,k,m) = (dryvol(m)/pi43)**third enddo ! compute wet radius for each mode do m=1,ntot_amode call modal_aero_kohler( & dryrad(i,k,m), hygro(m), rh(i,k), & wetrad(i,k,m), 1, 1 ) wetrad(i,k,m)=max(wetrad(i,k,m),dryrad(i,k,m)) dgncur_awet(i,k,m) = dgncur_a(i,k,m)* & (wetrad(i,k,m)/dryrad(i,k,m)) wetvol(m) = pi43*wetrad(i,k,m)*wetrad(i,k,m)*wetrad(i,k,m) wetvol(m) = max(wetvol(m),dryvol(m)) wtrvol(m) = wetvol(m) - dryvol(m) wtrvol(m) = max( wtrvol(m), 0.0_r8 ) ! apply simple treatment of deliquesence/crystallization hysteresis ! for rhcrystal < rh < rhdeliques, aerosol water is a fraction of ! the "upper curve" value, and the fraction is a linear function of rh if (rh(i,k) < rhcrystal_amode(m)) then wetrad(i,k,m) = dryrad(i,k,m) wetvol(m) = dryvol(m) wtrvol(m) = 0.0_r8 else if (rh(i,k) < rhdeliques_amode(m)) then wtrvol(m) = wtrvol(m)*hystfac(m) & *(rh(i,k) - rhcrystal_amode(m)) wtrvol(m) = max( wtrvol(m), 0.0_r8 ) wetvol(m) = dryvol(m) + wtrvol(m) wetrad(i,k,m) = (wetvol(m)/pi43)**third end if ! compute aer. water tendency = (new_water - old_water)/deltat ! [ either (kg-h2o/kg-air/s) or (mol-h2o/mol-air/s) ] ! lwater = lwaterptr_amode(m) - loffset if ( aero_mmr_flag ) then duma = 1.0_r8 else #if (defined MIMIC_CAM3) duma = mwdry/18.0_r8 #else duma = mwdry/mwh2o #endif end if qwater = density_water*naer(i,k,m)*wtrvol(m)*duma ! old_water (after modal_aero_calcsize) is ! qwater_old = raer(i,k,lwater) + raertend(i,k,lwater)*deltat ! and water tendency is ! raertend(i,k,lwater) = raertend(i,k,lwater) & ! + (qwater - qwater_old)/deltat ! which is equivalent to ! raertend(i,k,lwater) = (qwater - raer(i,k,lwater))/deltat qaerwat(i,k,m) = qwater ! compute aerosol wet density (kg/m3) if (wetvol(m) > 1.0e-30_r8) then wetdens(i,k,m) = (drymass(m) + density_water*wtrvol(m))/wetvol(m) else wetdens(i,k,m) = specdens_amode( lspectype_amode(1,m) ) end if enddo end do ! "i=1,ncol" end do ! "k=1,pver" ! output to history do m = 1, ntot_amode write( trnum, '(i3.3)' ) m call outfld( 'wat_a'//trnum(3:3), qaerwat(:,:,m), pcols, lchnk) ! note - eventually we should change these from "dgnd_a0N" to "dgnd_aN" call outfld( 'dgnd_a'//trnum(2:3), dgncur_a(:,:,m), pcols, lchnk) call outfld( 'dgnw_a'//trnum(2:3), dgncur_awet(:,:,m), pcols, lchnk) end do return end subroutine modal_aero_wateruptake_sub !----------------------------------------------------------------------- subroutine modal_aero_kohler( & rdry_in, hygro, s, rwet_out, im, imx ) ! calculates equlibrium radius r of haze droplets as function of ! dry particle mass and relative humidity s using kohler solution ! given in pruppacher and klett (eqn 6-35) ! for multiple aerosol types, assumes an internal mixture of aerosols implicit none ! arguments integer :: im ! number of grid points to be processed integer :: imx ! dimensioned number of grid points real(r8) :: rdry_in(imx) ! aerosol dry radius (m) real(r8) :: hygro(imx) ! aerosol volume-mean hygroscopicity (--) real(r8) :: s(imx) ! relative humidity (1 = saturated) real(r8) :: rwet_out(imx) ! aerosol wet radius (m) ! local variables integer, parameter :: imax=200 integer :: i, n, nsol real(r8) :: a, b real(r8) :: p40(imax),p41(imax),p42(imax),p43(imax) ! coefficients of polynomial real(r8) :: p30(imax),p31(imax),p32(imax) ! coefficients of polynomial real(r8) :: p real(r8) :: r3, r4 real(r8) :: r(imx) ! wet radius (microns) real(r8) :: rdry(imax) ! radius of dry particle (microns) real(r8) :: ss ! relative humidity (1 = saturated) real(r8) :: slog(imax) ! log relative humidity real(r8) :: vol(imax) ! total volume of particle (microns**3) real(r8) :: xi, xr complex(r8) :: cx4(4,imax),cx3(3,imax) real(r8), parameter :: eps = 1.e-4 real(r8), parameter :: mw = 18. real(r8), parameter :: pi = 3.14159 real(r8), parameter :: rhow = 1. real(r8), parameter :: surften = 76. real(r8), parameter :: tair = 273. real(r8), parameter :: third = 1./3. real(r8), parameter :: ugascon = 8.3e7 ! effect of organics on surface tension is neglected a=2.e4*mw*surften/(ugascon*tair*rhow) do i=1,im rdry(i) = rdry_in(i)*1.0e6 ! convert (m) to (microns) vol(i) = rdry(i)**3 ! vol is r**3, not volume b = vol(i)*hygro(i) ! quartic ss=min(s(i),1.-eps) ss=max(ss,1.e-10_r8) slog(i)=log(ss) p43(i)=-a/slog(i) p42(i)=0. p41(i)=b/slog(i)-vol(i) p40(i)=a*vol(i)/slog(i) ! cubic for rh=1 p32(i)=0. p31(i)=-b/a p30(i)=-vol(i) end do do 100 i=1,im ! if(vol(i).le.1.e-20)then if(vol(i).le.1.e-12)then r(i)=rdry(i) go to 100 endif p=abs(p31(i))/(rdry(i)*rdry(i)) if(p.lt.eps)then ! approximate solution for small particles r(i)=rdry(i)*(1.+p*third/(1.-slog(i)*rdry(i)/a)) else call makoh_quartic(cx4(1,i),p43(i),p42(i),p41(i),p40(i),1) ! find smallest real(r8) solution r(i)=1000.*rdry(i) nsol=0 do n=1,4 xr=real(cx4(n,i)) xi=imag(cx4(n,i)) if(abs(xi).gt.abs(xr)*eps) cycle if(xr.gt.r(i)) cycle if(xr.lt.rdry(i)*(1.-eps)) cycle if(xr.ne.xr) cycle r(i)=xr nsol=n end do if(nsol.eq.0)then write(iulog,*) & 'ccm kohlerc - no real(r8) solution found (quartic)' #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*)'roots =', (cx4(n,i),n=1,4) #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*)'p0-p3 =', p40(i), p41(i), p42(i), p43(i) #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*)'rh=',s(i) #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*)'setting radius to dry radius=',rdry(i) #ifdef WRF_PORT call wrf_message(iulog) #endif r(i)=rdry(i) ! stop endif endif if(s(i).gt.1.-eps)then ! save quartic solution at s=1-eps r4=r(i) ! cubic for rh=1 p=abs(p31(i))/(rdry(i)*rdry(i)) if(p.lt.eps)then r(i)=rdry(i)*(1.+p*third) else call makoh_cubic(cx3,p32,p31,p30,im) ! find smallest real(r8) solution r(i)=1000.*rdry(i) nsol=0 do n=1,3 xr=real(cx3(n,i)) xi=imag(cx3(n,i)) if(abs(xi).gt.abs(xr)*eps) cycle if(xr.gt.r(i)) cycle if(xr.lt.rdry(i)*(1.-eps)) cycle if(xr.ne.xr) cycle r(i)=xr nsol=n end do if(nsol.eq.0)then write(iulog,*) & 'ccm kohlerc - no real(r8) solution found (cubic)' #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*)'roots =', (cx3(n,i),n=1,3) #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*)'p0-p2 =', p30(i), p31(i), p32(i) #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*)'rh=',s(i) #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*)'setting radius to dry radius=',rdry(i) #ifdef WRF_PORT call wrf_message(iulog) #endif r(i)=rdry(i) ! stop endif endif r3=r(i) ! now interpolate between quartic, cubic solutions r(i)=(r4*(1.-s(i))+r3*(s(i)-1.+eps))/eps endif 100 continue ! bound and convert from microns to m do i=1,im r(i) = min(r(i),30._r8) ! upper bound based on 1 day lifetime rwet_out(i) = r(i)*1.e-6 end do return end subroutine modal_aero_kohler !----------------------------------------------------------------------- subroutine makoh_cubic( cx, p2, p1, p0, im ) ! ! solves x**3 + p2 x**2 + p1 x + p0 = 0 ! where p0, p1, p2 are real ! integer, parameter :: imx=200 integer :: im real(r8) :: p0(imx), p1(imx), p2(imx) complex(r8) :: cx(3,imx) integer :: i real(r8) :: eps, q(imx), r(imx), sqrt3, third complex(r8) :: ci, cq, crad(imx), cw, cwsq, cy(imx), cz(imx) save eps data eps/1.e-20/ third=1./3. ci=dcmplx(0.,1.) sqrt3=sqrt(3.) cw=0.5*(-1+ci*sqrt3) cwsq=0.5*(-1-ci*sqrt3) do i=1,im if(p1(i).eq.0.)then ! completely insoluble particle cx(1,i)=(-p0(i))**third cx(2,i)=cx(1,i) cx(3,i)=cx(1,i) else q(i)=p1(i)/3. r(i)=p0(i)/2. crad(i)=r(i)*r(i)+q(i)*q(i)*q(i) crad(i)=sqrt(crad(i)) cy(i)=r(i)-crad(i) if (abs(cy(i)).gt.eps) cy(i)=cy(i)**third cq=q(i) cz(i)=-cq/cy(i) cx(1,i)=-cy(i)-cz(i) cx(2,i)=-cw*cy(i)-cwsq*cz(i) cx(3,i)=-cwsq*cy(i)-cw*cz(i) endif enddo return end subroutine makoh_cubic !----------------------------------------------------------------------- subroutine makoh_quartic( cx, p3, p2, p1, p0, im ) ! solves x**4 + p3 x**3 + p2 x**2 + p1 x + p0 = 0 ! where p0, p1, p2, p3 are real ! integer, parameter :: imx=200 integer :: im real(r8) :: p0(imx), p1(imx), p2(imx), p3(imx) complex(r8) :: cx(4,imx) integer :: i real(r8) :: third, q(imx), r(imx) complex(r8) :: cb(imx), cb0(imx), cb1(imx), & crad(imx), cy(imx), czero czero=cmplx(0.0,0.0) third=1./3. do 10 i=1,im q(i)=-p2(i)*p2(i)/36.+(p3(i)*p1(i)-4*p0(i))/12. r(i)=-(p2(i)/6)**3+p2(i)*(p3(i)*p1(i)-4*p0(i))/48. & +(4*p0(i)*p2(i)-p0(i)*p3(i)*p3(i)-p1(i)*p1(i))/16 crad(i)=r(i)*r(i)+q(i)*q(i)*q(i) crad(i)=sqrt(crad(i)) cb(i)=r(i)-crad(i) if(cb(i).eq.czero)then ! insoluble particle cx(1,i)=(-p1(i))**third cx(2,i)=cx(1,i) cx(3,i)=cx(1,i) cx(4,i)=cx(1,i) else cb(i)=cb(i)**third cy(i)=-cb(i)+q(i)/cb(i)+p2(i)/6 cb0(i)=sqrt(cy(i)*cy(i)-p0(i)) cb1(i)=(p3(i)*cy(i)-p1(i))/(2*cb0(i)) cb(i)=p3(i)/2+cb1(i) crad(i)=cb(i)*cb(i)-4*(cy(i)+cb0(i)) crad(i)=sqrt(crad(i)) cx(1,i)=(-cb(i)+crad(i))/2. cx(2,i)=(-cb(i)-crad(i))/2. cb(i)=p3(i)/2-cb1(i) crad(i)=cb(i)*cb(i)-4*(cy(i)-cb0(i)) crad(i)=sqrt(crad(i)) cx(3,i)=(-cb(i)+crad(i))/2. cx(4,i)=(-cb(i)-crad(i))/2. endif 10 continue return end subroutine makoh_quartic !---------------------------------------------------------------------- end module modal_aero_wateruptake