chem/module_cam_mam_wateruptake.F

   1 #define MODAL_AERO
   2 #define WRF_PORT
   3 ! module_cam_mam_wateruptake.F
   4 ! adapted from cam3 modal_aero_wateruptake.F90 by r.c.easter, june 2010
   5 !Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
   6 ! 2010-06-24 rce notes
   7 !--------------------------------------------------------------
   8 ! modal_aero_wateruptake.F90
   9
  10
  11 !----------------------------------------------------------------------
  12 !BOP
  13 !
  14 ! !MODULE: modal_aero_wateruptake --- modal aerosol mode merging (renaming)
  15 !
  16 ! !INTERFACE:
  17    module modal_aero_wateruptake
  18
  19 ! !USES:
  20    use shr_kind_mod, only : r8 => shr_kind_r8
  21 #ifndef WRF_PORT
  22    use cam_logfile,  only: iulog
  23 #else
  24    use module_cam_support, only: iulog
  25 #endif
  26
  27    implicit none
  28    private
  29    save
  30
  31 !  !PUBLIC MEMBER FUNCTIONS:
  32    public modal_aero_wateruptake_sub
  33
  34 ! !PUBLIC DATA MEMBERS:
  35 !  currently none
  36
  37
  38 ! !DESCRIPTION: This module implements ...
  39 !
  40 ! !REVISION HISTORY:
  41 !
  42 !   RCE 07.04.13:  Adapted from MIRAGE2 code
  43 !
  44 !EOP
  45 !----------------------------------------------------------------------
  46 !BOC
  47 ! list private module data here
  48 !EOC
  49 !----------------------------------------------------------------------
  50
  51
  52    contains
  53
  54
  55 !----------------------------------------------------------------------
  56       subroutine modal_aero_wateruptake_sub(                &
  57                  lchnk, ncol, nstep,                        &
  58                  iwaterup_flag, loffset,                    &
  59                  aero_mmr_flag, h2o_mmr_flag,               &
  60                  deltat, h2ommr, t, pmid, pdel, cldn,       &
  61                  raer, raertend, dotend, qaerwat,           &
  62                  dgncur_a, dgncur_awet, wetdens             &
  63 #ifdef OBSRH
  64                  , obsrh                                    &
  65 #endif
  66                                                             )
  67 ! following are local for now -- wait and see
  68 !                maer, naer, wetrad, density                )
  69
  70 !-----------------------------------------------------------------------
  71 !
  72 ! Purpose: Compute aerosol wet radius
  73 !
  74 ! Method:  Kohler theory
  75 !
  76 ! Author:  S. Ghan
  77 !
  78 !-----------------------------------------------------------------------
  79       use modal_aero_data
  80 #ifndef WRF_PORT
  81       use mo_constants,  only: pi
  82       use pmgrid,        only: plat, plon
  83       use ppgrid,        only: begchunk, endchunk, pcols, pver
  84 #endif
  85       use physconst,     only: cpair, epsilo, gravit, mwdry, mwh2o,   &
  86                                rair, rga, rhoh2o, rh2o, latvap
  87 #ifndef WRF_PORT
  88       use constituents,  only: pcnst, qmin
  89 #endif
  90       use wv_saturation, only: aqsat, qsat_water
  91 #ifndef WRF_PORT
  92       use phys_grid,     only: get_lat_all_p, get_lon_all_p
  93
  94       use abortutils,    only : endrun
  95       use cam_history,   only : outfld
  96       use spmd_utils,    only : masterproc
  97 #else
  98       use module_cam_support, only: pcnst => pcnst_runtime, &
  99                                     pcols, pver, &
 100                                     endrun, masterproc, outfld
 101       use physconst,     only: pi
 102 #endif
 103
 104       implicit none
 105
 106       integer,  intent(in)  :: lchnk              ! chunk index
 107       integer,  intent(in)  :: ncol               ! number of columns
 108       integer,  intent(in)  :: nstep              ! time step
 109       integer,  intent(in)  :: iwaterup_flag
 110                                ! identifies call location from typhysbc (1,2)
 111       integer,  intent(in)  :: loffset            ! offset applied to modal aero "pointers"
 112
 113       logical,  intent(in)  :: aero_mmr_flag      ! if .true.,  aerosol q are kg/kg-air
 114                                                   ! if .false., aerosol q are mol/mol-air
 115       logical,  intent(in)  :: h2o_mmr_flag       ! if .true.,  h2ommr is kg/kg-air
 116       logical,  intent(inout)::dotend(pcnst)
 117                                ! identifies species for which tendencies are computed
 118
 119       real(r8), intent(in)  :: deltat             ! time step (s)
 120       real(r8), intent(in)  :: h2ommr(pcols,pver) ! layer specific humidity
 121       real(r8), intent(in)  :: t(pcols,pver)      ! layer temperatures (K)
 122       real(r8), intent(in)  :: pmid(pcols,pver)   ! layer pressure (Pa)
 123       real(r8), intent(in)  :: pdel(pcols,pver)   ! layer pressure thickness (Pa)
 124       real(r8), intent(in)  :: cldn(pcols,pver)   ! layer cloud fraction (0-1)
 125       real(r8), intent(in)  :: raer(pcols,pver,pcnst)
 126                                ! aerosol species MRs (kg/kg and #/kg)
 127       real(r8), intent(inout)::raertend(pcols,pver,pcnst)
 128                                ! aerosol MR tendencies (kg/kg/s)
 129                                ! only defined for aerosol water
 130       real(r8), intent(out)   :: qaerwat(pcols,pver,ntot_amode)
 131       real(r8), intent(in)    :: dgncur_a(pcols,pver,ntot_amode)
 132       real(r8), intent(out)   :: dgncur_awet(pcols,pver,ntot_amode)
 133       real(r8), intent(out)   :: wetdens(pcols,pver,ntot_amode)
 134
 135 ! following are local for now -- wait and see
 136 !     real(r8), intent(out) :: maer(pcols,pver,ntot_amode)
 137 !                              ! aerosol wet mass MR (including water) (kg/kg-air)
 138 !     real(r8), intent(out) :: naer(pcols,pver,ntot_amode)
 139 !                              ! aerosol number MR (bounded!) (#/kg-air)
 140 !     real(r8), intent(out) :: wetrad(pcols,pver,ntot_amode)
 141 !                              ! wet radius of aerosol (m)
 142 !     real(r8), intent(out) :: density(pcols,pver,ntot_amode)
 143 !                              ! wet mean density of aerosol (kg/m3)
 144
 145 !     local variables
 146
 147       integer i,k,m
 148       integer icol_diag
 149       integer l ! species index
 150       integer lmass ! pointer for aerosol mass
 151       integer ltype ! pointer for aerosol type
 152 !     integer lwater ! pointer for water on aerosol
 153       integer  lat(pcols), lon(pcols)      ! lat,lon indices
 154
 155       real(r8) density_water                   ! density of water (kg/m3)
 156       real(r8) drydens(ntot_amode)   ! dry particle density  (kg/m^3)
 157       real(r8) drymass(ntot_amode)   ! single-particle-mean dry mass  (kg)
 158       real(r8) dryrad(pcols,pver,ntot_amode) ! dry volume mean radius of aerosol (m)
 159       real(r8) dryvol(ntot_amode)    ! single-particle-mean dry volume (m3)
 160       real(r8) dryvolmr(ntot_amode)  ! volume MR for aerosol mode (m3/kg)
 161       real(r8) duma, dumb
 162       real(r8) es(pcols,pver)        ! saturation vapor pressure (Pa)
 163       real(r8) hygro(ntot_amode)     ! volume-weighted mean hygroscopicity (--)
 164       real(r8) hystfac(ntot_amode)   ! working variable for hysteresis
 165       real(r8) pi43
 166       real(r8) qs(pcols,pver)        ! saturation specific humidity
 167       real(r8) qwater                ! aerosol water MR
 168       real(r8) rh(pcols,pver)        ! relative humidity (0-1)
 169       real(r8) third
 170       real(r8) v2ncur_a(pcols,pver,ntot_amode)
 171       real(r8) wtrvol(ntot_amode)    ! single-particle-mean water volume in wet aerosol (m3)
 172       real(r8) wetvol(ntot_amode)    ! single-particle-mean wet volume (m3)
 173
 174       real(r8) :: maer(pcols,pver,ntot_amode)
 175                               ! aerosol wet mass MR (including water) (kg/kg-air)
 176       real(r8) :: naer(pcols,pver,ntot_amode)
 177                               ! aerosol number MR (bounded!) (#/kg-air)
 178       real(r8) :: wetrad(pcols,pver,ntot_amode)
 179                               ! wet radius of aerosol (m)
 180
 181       character(len=3) :: trnum       ! used to hold mode number (as characters)
 182
 183 !-----------------------------------------------------------------------
 184
 185
 186
 187 ! set the dotend's
 188 ! currently, this is one of several routines called from aerosol_wet_intr
 189 ! so DO NOT re-initialize dotend=.false. and raertend=0.0
 190 !     dotend(:) = .false.
 191       do m=1,ntot_amode
 192 !        lwater = lwaterptr_amode(m) - loffset
 193 !        dotend(lwater) = .true.
 194 !        raertend(1:ncol,:,lwater) = 0.0
 195       end do
 196
 197
 198       third=1./3.
 199       pi43 = pi*4.0/3.0
 200       density_water = rhoh2o   ! is (kg/m3)
 201
 202 ! compute size-related factors
 203 !    when numptr_amode(m)=0, earlier code set dryrad = prescribed volume.
 204 !    this is no longer needed because
 205 !           dryrad(i,k,m) = (1.0/v2ncur_a(i,k,m)/pi43)**third
 206 !    works in all cases
 207       do m=1,ntot_amode
 208            hystfac(m) = 1.0 / max( 1.0e-5_r8,   &
 209                               (rhdeliques_amode(m)-rhcrystal_amode(m)) )
 210       enddo
 211
 212 ! main loops over i, k
 213
 214 !     call aqsat( t, pmid, es, qs, pcols, ncol, pver, 1, pver )
 215       do k=1,pver
 216       do i=1,ncol
 217
 218          qs(i,k)=qsat_water(t(i,k),pmid(i,k))
 219
 220          if ( h2o_mmr_flag ) then
 221             rh(i,k) = h2ommr(i,k)/qs(i,k)
 222          else
 223             rh(i,k) = h2ommr(i,k)*mwh2o/(mwdry*qs(i,k))
 224          end if
 225          rh(i,k) = max(rh(i,k),0.0_r8)
 226          rh(i,k) = min(rh(i,k),0.98_r8)
 227          if (cldn(i,k) .lt. 1.0_r8) then
 228            rh(i,k) = (rh(i,k) - cldn(i,k)) / (1.0_r8 - cldn(i,k))  ! clear portion
 229          end if
 230          rh(i,k) = max(rh(i,k),0.0_r8)
 231
 232
 233 !     compute dryvolmr, maer, naer for each mode
 234          do m=1,ntot_amode
 235
 236             maer(i,k,m)=0.
 237             dryvolmr(m)=0.
 238             hygro(m)=0.
 239             do l = 1, nspec_amode(m)
 240                lmass = lmassptr_amode(l,m) - loffset
 241                ltype = lspectype_amode(l,m)
 242                if ( aero_mmr_flag ) then
 243                   duma = raer(i,k,lmass)
 244                else
 245                   duma = raer(i,k,lmass)*(specmw_amode(ltype)/mwdry)
 246                end if
 247                maer(i,k,m) = maer(i,k,m) + duma
 248                dumb = duma/specdens_amode(ltype)
 249                dryvolmr(m) = dryvolmr(m) + dumb
 250                hygro(m) = hygro(m) + dumb*spechygro(ltype)
 251             enddo
 252             if (dryvolmr(m) > 1.0e-30_r8) then
 253                hygro(m) = hygro(m)/dryvolmr(m)
 254             else
 255                hygro(m) = spechygro( lspectype_amode(1,m) )
 256             end if
 257
 258 !     naer = aerosol number (#/kg)
 259 !     the new v2ncur_a replaces old coding here
 260             v2ncur_a(i,k,m) = 1. / ( (pi/6.)*                            &
 261                 (dgncur_a(i,k,m)**3.)*exp(4.5*alnsg_amode(m)**2.) )
 262             naer(i,k,m) = dryvolmr(m)*v2ncur_a(i,k,m)
 263          enddo   !m=1,ntot_amode
 264
 265 !     compute mean (1 particle) dry volume and mass for each mode
 266 !     old coding is replaced because the new (1/v2ncur_a) is equal to
 267 !        the mean particle volume
 268 !     also moletomass forces maer >= 1.0e-30, so (maer/dryvolmr)
 269 !        should never cause problems (but check for maer < 1.0e-31 anyway)
 270          do m=1,ntot_amode
 271             if (maer(i,k,m) .gt. 1.0e-31) then
 272                drydens(m) = maer(i,k,m)/dryvolmr(m)
 273             else
 274                drydens(m) = 1.0
 275             end if
 276             dryvol(m) = 1.0/v2ncur_a(i,k,m)
 277             drymass(m) = drydens(m)*dryvol(m)
 278             dryrad(i,k,m) = (dryvol(m)/pi43)**third
 279          enddo
 280
 281 !     compute wet radius for each mode
 282          do m=1,ntot_amode
 283             call modal_aero_kohler(                          &
 284                     dryrad(i,k,m), hygro(m), rh(i,k),        &
 285                     wetrad(i,k,m), 1, 1                      )
 286
 287             wetrad(i,k,m)=max(wetrad(i,k,m),dryrad(i,k,m))
 288             dgncur_awet(i,k,m) = dgncur_a(i,k,m)*   &
 289                                      (wetrad(i,k,m)/dryrad(i,k,m))
 290             wetvol(m) = pi43*wetrad(i,k,m)*wetrad(i,k,m)*wetrad(i,k,m)
 291             wetvol(m) = max(wetvol(m),dryvol(m))
 292             wtrvol(m) = wetvol(m) - dryvol(m)
 293             wtrvol(m) = max( wtrvol(m), 0.0_r8 )
 294
 295 !     apply simple treatment of deliquesence/crystallization hysteresis
 296 !     for rhcrystal < rh < rhdeliques, aerosol water is a fraction of
 297 !     the "upper curve" value, and the fraction is a linear function of rh
 298             if (rh(i,k) < rhcrystal_amode(m)) then
 299                wetrad(i,k,m) = dryrad(i,k,m)
 300                wetvol(m) = dryvol(m)
 301                wtrvol(m) = 0.0_r8
 302             else if (rh(i,k) < rhdeliques_amode(m)) then
 303                wtrvol(m) = wtrvol(m)*hystfac(m)   &
 304                                     *(rh(i,k) - rhcrystal_amode(m))
 305                wtrvol(m) = max( wtrvol(m), 0.0_r8 )
 306                wetvol(m) = dryvol(m) + wtrvol(m)
 307                wetrad(i,k,m) = (wetvol(m)/pi43)**third
 308             end if
 309
 310 !     compute aer. water tendency = (new_water - old_water)/deltat
 311 !     [ either (kg-h2o/kg-air/s) or (mol-h2o/mol-air/s) ]
 312 !           lwater = lwaterptr_amode(m) - loffset
 313             if ( aero_mmr_flag ) then
 314                duma = 1.0_r8
 315             else
 316 #if (defined MIMIC_CAM3)
 317                duma = mwdry/18.0_r8
 318 #else
 319                duma = mwdry/mwh2o
 320 #endif
 321             end if
 322             qwater = density_water*naer(i,k,m)*wtrvol(m)*duma
 323
 324 !     old_water (after modal_aero_calcsize) is
 325 !           qwater_old = raer(i,k,lwater) + raertend(i,k,lwater)*deltat
 326 !     and water tendency is
 327 !           raertend(i,k,lwater) = raertend(i,k,lwater)   &
 328 !                                + (qwater - qwater_old)/deltat
 329 !     which is equivalent to
 330 !           raertend(i,k,lwater) = (qwater - raer(i,k,lwater))/deltat
 331             qaerwat(i,k,m) = qwater
 332
 333 !     compute aerosol wet density (kg/m3)
 334             if (wetvol(m) > 1.0e-30_r8) then
 335                wetdens(i,k,m) = (drymass(m) + density_water*wtrvol(m))/wetvol(m)
 336             else
 337                wetdens(i,k,m) = specdens_amode( lspectype_amode(1,m) )
 338             end if
 339
 340          enddo
 341
 342
 343       end do   ! "i=1,ncol"
 344       end do   ! "k=1,pver"
 345
 346
 347 ! output to history
 348       do m = 1, ntot_amode
 349          write( trnum, '(i3.3)' ) m
 350          call outfld( 'wat_a'//trnum(3:3),  qaerwat(:,:,m),     pcols, lchnk)
 351 ! note - eventually we should change these from "dgnd_a0N" to "dgnd_aN"
 352          call outfld( 'dgnd_a'//trnum(2:3), dgncur_a(:,:,m),    pcols, lchnk)
 353          call outfld( 'dgnw_a'//trnum(2:3), dgncur_awet(:,:,m), pcols, lchnk)
 354       end do
 355
 356
 357       return
 358       end subroutine modal_aero_wateruptake_sub
 359
 360
 361 !-----------------------------------------------------------------------
 362       subroutine modal_aero_kohler(   &
 363           rdry_in, hygro, s, rwet_out, im, imx )
 364
 365 ! calculates equlibrium radius r of haze droplets as function of
 366 ! dry particle mass and relative humidity s using kohler solution
 367 ! given in pruppacher and klett (eqn 6-35)
 368
 369 ! for multiple aerosol types, assumes an internal mixture of aerosols
 370
 371       implicit none
 372
 373 ! arguments
 374       integer :: im         ! number of grid points to be processed
 375       integer :: imx        ! dimensioned number of grid points
 376       real(r8) :: rdry_in(imx)    ! aerosol dry radius (m)
 377       real(r8) :: hygro(imx)      ! aerosol volume-mean hygroscopicity (--)
 378       real(r8) :: s(imx)          ! relative humidity (1 = saturated)
 379       real(r8) :: rwet_out(imx)   ! aerosol wet radius (m)
 380
 381 ! local variables
 382       integer, parameter :: imax=200
 383       integer :: i, n, nsol
 384
 385       real(r8) :: a, b
 386       real(r8) :: p40(imax),p41(imax),p42(imax),p43(imax) ! coefficients of polynomial
 387       real(r8) :: p30(imax),p31(imax),p32(imax) ! coefficients of polynomial
 388       real(r8) :: p
 389       real(r8) :: r3, r4
 390       real(r8) :: r(imx)        ! wet radius (microns)
 391       real(r8) :: rdry(imax)    ! radius of dry particle (microns)
 392       real(r8) :: ss            ! relative humidity (1 = saturated)
 393       real(r8) :: slog(imax)    ! log relative humidity
 394       real(r8) :: vol(imax)     ! total volume of particle (microns**3)
 395       real(r8) :: xi, xr
 396
 397       complex(r8) :: cx4(4,imax),cx3(3,imax)
 398
 399       real(r8), parameter :: eps = 1.e-4
 400       real(r8), parameter :: mw = 18.
 401       real(r8), parameter :: pi = 3.14159
 402       real(r8), parameter :: rhow = 1.
 403       real(r8), parameter :: surften = 76.
 404       real(r8), parameter :: tair = 273.
 405       real(r8), parameter :: third = 1./3.
 406       real(r8), parameter :: ugascon = 8.3e7
 407
 408
 409 !     effect of organics on surface tension is neglected
 410       a=2.e4*mw*surften/(ugascon*tair*rhow)
 411
 412       do i=1,im
 413            rdry(i) = rdry_in(i)*1.0e6   ! convert (m) to (microns)
 414            vol(i) = rdry(i)**3          ! vol is r**3, not volume
 415            b = vol(i)*hygro(i)
 416
 417 !          quartic
 418            ss=min(s(i),1.-eps)
 419            ss=max(ss,1.e-10_r8)
 420            slog(i)=log(ss)
 421            p43(i)=-a/slog(i)
 422            p42(i)=0.
 423            p41(i)=b/slog(i)-vol(i)
 424            p40(i)=a*vol(i)/slog(i)
 425 !          cubic for rh=1
 426            p32(i)=0.
 427            p31(i)=-b/a
 428            p30(i)=-vol(i)
 429       end do
 430
 431
 432        do 100 i=1,im
 433
 434 !       if(vol(i).le.1.e-20)then
 435         if(vol(i).le.1.e-12)then
 436            r(i)=rdry(i)
 437            go to 100
 438         endif
 439
 440         p=abs(p31(i))/(rdry(i)*rdry(i))
 441         if(p.lt.eps)then
 442 !          approximate solution for small particles
 443            r(i)=rdry(i)*(1.+p*third/(1.-slog(i)*rdry(i)/a))
 444         else
 445            call makoh_quartic(cx4(1,i),p43(i),p42(i),p41(i),p40(i),1)
 446 !          find smallest real(r8) solution
 447            r(i)=1000.*rdry(i)
 448            nsol=0
 449            do n=1,4
 450               xr=real(cx4(n,i))
 451               xi=imag(cx4(n,i))
 452               if(abs(xi).gt.abs(xr)*eps) cycle
 453               if(xr.gt.r(i)) cycle
 454               if(xr.lt.rdry(i)*(1.-eps)) cycle
 455               if(xr.ne.xr) cycle
 456               r(i)=xr
 457               nsol=n
 458            end do
 459            if(nsol.eq.0)then
 460               write(iulog,*)   &
 461                'ccm kohlerc - no real(r8) solution found (quartic)'
 462 #ifdef WRF_PORT
 463               call wrf_message(iulog)
 464 #endif
 465               write(iulog,*)'roots =', (cx4(n,i),n=1,4)
 466 #ifdef WRF_PORT
 467               call wrf_message(iulog)
 468 #endif
 469               write(iulog,*)'p0-p3 =', p40(i), p41(i), p42(i), p43(i)
 470 #ifdef WRF_PORT
 471               call wrf_message(iulog)
 472 #endif
 473               write(iulog,*)'rh=',s(i)
 474 #ifdef WRF_PORT
 475               call wrf_message(iulog)
 476 #endif
 477               write(iulog,*)'setting radius to dry radius=',rdry(i)
 478 #ifdef WRF_PORT
 479               call wrf_message(iulog)
 480 #endif
 481               r(i)=rdry(i)
 482 !             stop
 483            endif
 484         endif
 485
 486         if(s(i).gt.1.-eps)then
 487 !          save quartic solution at s=1-eps
 488            r4=r(i)
 489 !          cubic for rh=1
 490            p=abs(p31(i))/(rdry(i)*rdry(i))
 491            if(p.lt.eps)then
 492               r(i)=rdry(i)*(1.+p*third)
 493            else
 494               call makoh_cubic(cx3,p32,p31,p30,im)
 495 !             find smallest real(r8) solution
 496               r(i)=1000.*rdry(i)
 497               nsol=0
 498               do n=1,3
 499                  xr=real(cx3(n,i))
 500                  xi=imag(cx3(n,i))
 501                  if(abs(xi).gt.abs(xr)*eps) cycle
 502                  if(xr.gt.r(i)) cycle
 503                  if(xr.lt.rdry(i)*(1.-eps)) cycle
 504                  if(xr.ne.xr) cycle
 505                  r(i)=xr
 506                  nsol=n
 507               end do
 508               if(nsol.eq.0)then
 509                  write(iulog,*)   &
 510                   'ccm kohlerc - no real(r8) solution found (cubic)'
 511 #ifdef WRF_PORT
 512                  call wrf_message(iulog)
 513 #endif
 514                  write(iulog,*)'roots =', (cx3(n,i),n=1,3)
 515 #ifdef WRF_PORT
 516                  call wrf_message(iulog)
 517 #endif
 518                  write(iulog,*)'p0-p2 =', p30(i), p31(i), p32(i)
 519 #ifdef WRF_PORT
 520                  call wrf_message(iulog)
 521 #endif
 522                  write(iulog,*)'rh=',s(i)
 523 #ifdef WRF_PORT
 524                  call wrf_message(iulog)
 525 #endif
 526                  write(iulog,*)'setting radius to dry radius=',rdry(i)
 527 #ifdef WRF_PORT
 528                  call wrf_message(iulog)
 529 #endif
 530                  r(i)=rdry(i)
 531 !                stop
 532               endif
 533            endif
 534            r3=r(i)
 535 !          now interpolate between quartic, cubic solutions
 536            r(i)=(r4*(1.-s(i))+r3*(s(i)-1.+eps))/eps
 537         endif
 538
 539   100 continue
 540
 541 ! bound and convert from microns to m
 542       do i=1,im
 543          r(i) = min(r(i),30._r8) ! upper bound based on 1 day lifetime
 544          rwet_out(i) = r(i)*1.e-6
 545       end do
 546
 547       return
 548       end subroutine modal_aero_kohler
 549
 550
 551 !-----------------------------------------------------------------------
 552       subroutine makoh_cubic( cx, p2, p1, p0, im )
 553 !
 554 !     solves  x**3 + p2 x**2 + p1 x + p0 = 0
 555 !     where p0, p1, p2 are real
 556 !
 557       integer, parameter :: imx=200
 558       integer :: im
 559       real(r8) :: p0(imx), p1(imx), p2(imx)
 560       complex(r8) :: cx(3,imx)
 561
 562       integer :: i
 563       real(r8) :: eps, q(imx), r(imx), sqrt3, third
 564       complex(r8) :: ci, cq, crad(imx), cw, cwsq, cy(imx), cz(imx)
 565
 566       save eps
 567       data eps/1.e-20/
 568
 569       third=1./3.
 570       ci=dcmplx(0.,1.)
 571       sqrt3=sqrt(3.)
 572       cw=0.5*(-1+ci*sqrt3)
 573       cwsq=0.5*(-1-ci*sqrt3)
 574
 575       do i=1,im
 576       if(p1(i).eq.0.)then
 577 !        completely insoluble particle
 578          cx(1,i)=(-p0(i))**third
 579          cx(2,i)=cx(1,i)
 580          cx(3,i)=cx(1,i)
 581       else
 582          q(i)=p1(i)/3.
 583          r(i)=p0(i)/2.
 584          crad(i)=r(i)*r(i)+q(i)*q(i)*q(i)
 585          crad(i)=sqrt(crad(i))
 586
 587          cy(i)=r(i)-crad(i)
 588          if (abs(cy(i)).gt.eps) cy(i)=cy(i)**third
 589          cq=q(i)
 590          cz(i)=-cq/cy(i)
 591
 592          cx(1,i)=-cy(i)-cz(i)
 593          cx(2,i)=-cw*cy(i)-cwsq*cz(i)
 594          cx(3,i)=-cwsq*cy(i)-cw*cz(i)
 595       endif
 596       enddo
 597
 598       return
 599       end subroutine makoh_cubic
 600
 601
 602 !-----------------------------------------------------------------------
 603       subroutine makoh_quartic( cx, p3, p2, p1, p0, im )
 604
 605 !     solves x**4 + p3 x**3 + p2 x**2 + p1 x + p0 = 0
 606 !     where p0, p1, p2, p3 are real
 607 !
 608       integer, parameter :: imx=200
 609       integer :: im
 610       real(r8) :: p0(imx), p1(imx), p2(imx), p3(imx)
 611       complex(r8) :: cx(4,imx)
 612
 613       integer :: i
 614       real(r8) :: third, q(imx), r(imx)
 615       complex(r8) :: cb(imx), cb0(imx), cb1(imx),   &
 616                      crad(imx), cy(imx), czero
 617
 618
 619       czero=cmplx(0.0,0.0)
 620       third=1./3.
 621
 622       do 10 i=1,im
 623
 624       q(i)=-p2(i)*p2(i)/36.+(p3(i)*p1(i)-4*p0(i))/12.
 625       r(i)=-(p2(i)/6)**3+p2(i)*(p3(i)*p1(i)-4*p0(i))/48.   &
 626        +(4*p0(i)*p2(i)-p0(i)*p3(i)*p3(i)-p1(i)*p1(i))/16
 627
 628       crad(i)=r(i)*r(i)+q(i)*q(i)*q(i)
 629       crad(i)=sqrt(crad(i))
 630
 631       cb(i)=r(i)-crad(i)
 632       if(cb(i).eq.czero)then
 633 !        insoluble particle
 634          cx(1,i)=(-p1(i))**third
 635          cx(2,i)=cx(1,i)
 636          cx(3,i)=cx(1,i)
 637          cx(4,i)=cx(1,i)
 638       else
 639          cb(i)=cb(i)**third
 640
 641          cy(i)=-cb(i)+q(i)/cb(i)+p2(i)/6
 642
 643          cb0(i)=sqrt(cy(i)*cy(i)-p0(i))
 644          cb1(i)=(p3(i)*cy(i)-p1(i))/(2*cb0(i))
 645
 646          cb(i)=p3(i)/2+cb1(i)
 647          crad(i)=cb(i)*cb(i)-4*(cy(i)+cb0(i))
 648          crad(i)=sqrt(crad(i))
 649          cx(1,i)=(-cb(i)+crad(i))/2.
 650          cx(2,i)=(-cb(i)-crad(i))/2.
 651
 652          cb(i)=p3(i)/2-cb1(i)
 653          crad(i)=cb(i)*cb(i)-4*(cy(i)-cb0(i))
 654          crad(i)=sqrt(crad(i))
 655          cx(3,i)=(-cb(i)+crad(i))/2.
 656          cx(4,i)=(-cb(i)-crad(i))/2.
 657       endif
 658    10 continue
 659
 660       return
 661       end subroutine makoh_quartic
 662
 663 !----------------------------------------------------------------------
 664
 665    end module modal_aero_wateruptake
 666
 667