186 & wparc_ls, Aer_Props, npre_in, dpre_in, ccn_diagr8, Ndropr8, &
187 & cdncr8, smaxliqr8, incr8, smaxicer8, nheticer8, INimmr8, &
188 & dINimmr8, Ncdepr8, Ncdhfr8, sc_icer8, fdust_immr8, fdust_depr8, &
189 & fdust_dhfr8, nlimr8, use_average_v, CCN_param, IN_param, fd_dust,&
190 & fd_soot, pfrz_inc_r8, sigma_nuc, rhi_cell,nccn)
196 type(
aerprops),
intent(in) :: aer_props
198 logical :: use_average_v
201 real(r8),
intent(in) :: tparc_in, pparc_in, sigwparc_in, &
202 & wparc_ls, npre_in, dpre_in, Ndropr8, fd_soot, fd_dust, &
203 & sigma_nuc, rhi_cell
204 integer,
intent(in) :: ccn_param, in_param, nccn
206 real(r8),
dimension(:),
intent(inout) :: ccn_diagr8
208 real(r8),
intent(out) :: cdncr8, smaxliqr8, incr8, smaxicer8, &
209 & nheticer8, INimmr8, dINimmr8, Ncdepr8, Ncdhfr8, sc_icer8, &
210 & fdust_immr8, fdust_depr8, fdust_dhfr8, nlimr8, pfrz_inc_r8
215 integer :: k, n, i, j, naux
218 real*8 :: nact, wparc, tparc,pparc, accom,sigw, smax, antot, &
221 real*8,
dimension(nccn) :: smax_diag
224 real*8 :: nhet, nice, smaxice, nlim, air_den, frac, norg, nbc, &
225 & nhom, dorg, dbc, kappa, inimm, dinimm, aux
229 real*8,
dimension(nsmx_par) :: sg_par, tp_par, dpg_par
230 &, sig_par, vhf_par, ams_par
231 &, dens_par, deni_par, amfs_par
232 &, kappa_par, ndust_ice, sigdust_ice
236 real*8 :: temp_par, pres_par
237 real*8 :: akoh_par, alfa_par, bet2_par
238 real*8 aka_par, dv_par, psat_par, dair_par,surt_par,ddry_ice,
239 & np_ice,nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
240 & g1_ice, g2_ice,gdoin_ice, z_ice, norg_ice, sigorg_ice,
241 & dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
242 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
243 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
244 & si0bc_ice,nbc_ice,sc_ice,nhet_dep,fdust_imm,fdust_dep,fdust_dhf,
245 & waux_ice,fdrop_dust,fdrop_bc,d_preex, n_preex,
246 & one_over_tao_preex,p_ice, t_ice,miuv_ice,nhet_dhf,denice_ice,
247 & vpresw_ice,vpresi_ice,denair_ice
249 integer :: nmodes,act_param,nbindust_ice
255 sigwparc = sigwparc_in
257 air_den = pparc*28.8d-3/rgas_par/tparc
258 n_preex = max(npre_in*air_den, zero_par)
259 d_preex = max(dpre_in, 1.0e-9)
260 use_av_v = use_average_v
271 smaxice = max(2.349d0-(tparc/259d0) -1.0 , 0.0)
274 sc_ice = max(2.349d0-(tparc/259d0), 1.0)
275 If (tparc > thom) sc_ice =1.0
281 fdust_immr8 = zero_par
282 fdust_dhfr8 = zero_par
283 fdust_depr8 = zero_par
287 pfrz_inc_r8 = zero_par
295 if (sum(aer_props%num) <= 1.0e2)
then
299 nmodes = max(aer_props%nmods, 1)
318 dpg_par(n) = zero_par
319 vhf_par(n) = zero_par
320 ams_par(n) = zero_par
321 dens_par(n) = zero_par
324 amfs_par(n) = zero_par
325 deni_par(n) = zero_par
326 kappa_par(n) = zero_par
345 tp_par(n) = aer_props%num(n) * air_den
346 dpg_par(n) = max(aer_props%dpg(n), 1.0e-10)
347 sig_par(n) = aer_props%sig(n)
348 kappa_par(n) = max(aer_props%kap(n), 0.001)
349 dens_par(n) = aer_props%den(n)
351 if (kappa_par(n) > 0.01)
then
352 ams_par(n) = 18.0e-3*1.7*3.0/kappa_par(n)
354 ams_par(n) = 900.0e-3
358 deni_par(n) = dens_par(n)
359 antot = antot + tp_par(n)
368 temp_par = max(tparc, 245.0)
369 pres_par = max(pparc, 34000.0)
374 wparc = max(max(0.8d0*sigwparc, 0.01)+ wparc_ls, 0.01)
377 act_param = ccn_param
384 if (tparc > 245.0)
then
385 if (antot > 1.0)
then
387 call ccnspec (tparc,pparc,nmodes
388 &, amfs_par,dens_par,deni_par,vhf_par,ams_par
389 &, sg_par,tp_par,akoh_par,surt_par,temp_par,pres_par
390 &, dv_par,act_param,aka_par, psat_par,dair_par
393 if (wparc >= 0.005)
then
394 if (act_param > 1)
then
396 call arg_activ (wparc,0.d0,nact,smax,nmodes,tp_par
397 &, dpg_par,kappa_par,sig_par,temp_par, pres_par)
401 call pdfactiv (wparc,0.d0,nact,smax,nmodes
402 &, alfa_par,bet2_par,akoh_par,sg_par,tp_par
403 &, temp_par, pres_par,aka_par, dv_par, psat_par
404 &, dair_par,ntot,sig_par)
409 cdncr8 = max(nact/air_den, zero_par)
410 smaxliqr8 = max(smax, zero_par)
438 sigorg_ice = zero_par
455 if (kappa_par(k) > 0.1)
then
456 np_ice = np_ice + tp_par(k)
457 ddry_ice = ddry_ice + dpg_par(k)
461 ddry_ice = ddry_ice/max(naux , 1)
468 nbindust_ice = max(aeraux%nmods, 1)
476 sigdust_ice = zero_par
478 ddust_ice(n) = dble(aeraux%dpg(n))
479 ndust_ice(n) = dble(aeraux%num(n))*air_den
480 sigdust_ice(n) = dble(aeraux%sig(n))
487 naux = max(aeraux%nmods, 1)
488 dbc_ice = dble(sum(aeraux%dpg(1:naux)))/naux
489 nbc_ice = dble(sum(aeraux%num(1:naux)))*air_den
490 sigbc_ice = dble(sum(aeraux%sig(1:naux)))/naux
495 naux = max(aeraux%nmods, 1)
496 dorg_ice = dble(sum(aeraux%dpg(1:naux)))/naux
497 norg_ice = dble(sum(aeraux%num(1:naux)))*air_den
498 sigorg_ice = dble(sum(aeraux%sig(1:naux)))/naux
508 antot=sum(ndust_ice)+ norg_ice+ nbc_ice+ np_ice
510 sigwparc=max(0.01, sigwparc)
511 sigwparc=min(5.0, sigwparc)
512 waux_ice=max(wparc_ls + sigwparc*0.8, 0.01)
521 if (antot > 1.0e2)
then
522 if (tparc < to_ice)
then
524 CALL prop_ice(tparc, pparc, denice_ice,ddry_ice
525 &, nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice
526 &, g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice
527 &, kdust_ice, kbc_ice, shdust_ice, shbc_ice
528 &, effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice
529 &, si0bc_ice,nbc_ice,d_preex, n_preex
530 &, one_over_tao_preex,p_ice, t_ice,act_param
531 &, ndust_ice,vpresw_ice,vpresi_ice,denair_ice)
534 if (tparc > thom)
then
539 if (sum(ndust_ice)+ norg_ice+ nbc_ice .gt. 1.e3)
then
541 call inimmersion(inimm, dinimm, waux_ice,dbc_ice,sigbc_ice
542 &, nbc_ice,fdust_imm,fdrop_dust,fdrop_bc,ndust_ice
543 &, sigdust_ice,ddust_ice,nbindust_ice
544 &, vpresw_ice,vpresi_ice,t_ice)
547 ndust_ice = max(ndust_ice*(1.0-fdrop_dust), 0.0)
548 nbc_ice = max(nbc_ice*(1.0-fdrop_bc), 0.0)
551 call iceparam (sigwparc, denice_ice,ddry_ice,np_ice
552 &, nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice
553 &, g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,sc_ice
554 &, norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice
555 &, kdust_ice, kbc_ice, shdust_ice, shbc_ice
556 &, effdust_ice, effbc_ice, del1dust_ice, si0dust_ice
557 &, del1bc_ice, si0bc_ice,nbc_ice,one_over_tao_preex
558 &, nhet, nice, smaxice, nlim,nhet_dep,nhet_dhf,fdust_dep
559 &, p_ice,t_ice,ndust_ice,sigdust_ice,ddust_ice,nbindust_ice
560 &, use_av_v,miuv_ice,vpresw_ice,vpresi_ice,denair_ice)
567 call iceparam (sigwparc, denice_ice,ddry_ice,np_ice
568 &, nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice
569 &, g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,sc_ice
570 &, norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice
571 &, kdust_ice, kbc_ice, shdust_ice, shbc_ice
572 &, effdust_ice, effbc_ice, del1dust_ice, si0dust_ice
573 &, del1bc_ice, si0bc_ice,nbc_ice,one_over_tao_preex
574 &, nhet, nice, smaxice, nlim,nhet_dep,nhet_dhf,fdust_dep
575 &, p_ice, t_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice
576 &, use_av_v,miuv_ice,vpresw_ice,vpresi_ice,denair_ice)
580 aux = (sc_ice - rhi_cell)/(sq2_par*sigma_nuc)
582 pfrz_inc_r8 = 1.0d0- 0.5d0*(1.0d0+erf(aux))
583 pfrz_inc_r8 = min(max(pfrz_inc_r8, 0.0), 0.999)
590 If (tparc > thom) sc_ice =1.0
594 smaxicer8 = min(max(smaxice, zero_par), 2.0)
595 nheticer8 = min(max(nhet, zero_par), 1e10)
596 incr8 = min(max(nice/air_den, zero_par), 1e10)
597 nlimr8 = min(max(nlim, zero_par), 1e10)
598 sc_icer8 = min(max(sc_ice, 1.0), 2.0)
599 inimmr8 = min(max(inimm, zero_par), 1e10)
600 dinimmr8 = min(max(dinimm/air_den, zero_par), 1e10)
601 ncdepr8 = min(max(nhet_dep, zero_par), 1e10)
602 ncdhfr8 = min(max(nhet_dhf, zero_par), 1e10)
603 fdust_immr8 = min(max(fdust_imm, zero_par), 1e10)
604 fdust_depr8 = min(max(fdust_dep, zero_par), 1e10)
605 fdust_dhfr8 = min(max(fdust_dhf, zero_par), 1e10)
1493 subroutine activate (wparc,ndroplet,smax,nmodes,
1494 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,temp_par, pres_par,
1495 & aka_par, dv_par, psat_par,dair_par,ntot,sig_par)
1498 integer :: i,niter,nmodes
1500 real*8 :: ndrpl, wparc, beta,cf1, cf2,x1,sinteg1,sinteg2,y1, x2,
1501 &y2,x3,y3, sign,smax, ent_par, ndroplet, nrdpl,bet1_par,
1502 & alfa_par,bet2_par,akoh_par,sg_par(nmodes),tp_par(nmodes),
1503 & temp_par, pres_par,aka_par, dv_par, psat_par,dair_par,ntot,
1504 & wparcel,sig_par(nmodes)
1516 alfa_par = grav_par*amw_par*dhv_par/cpair_par/rgas_par/temp_par /
1517 &temp_par -grav_par*ama_par/rgas_par/temp_par
1519 bet1_par = pres_par*ama_par/psat_par/amw_par + amw_par*dhv_par*
1520 & dhv_par/cpair_par/rgas_par/temp_par/temp_par
1521 bet2_par = rgas_par*temp_par*denw_par/psat_par/dv_par/amw_par/
1522 & 4d0 +dhv_par*denw_par/4d0/aka_par/temp_par*(dhv_par* amw_par/
1523 &rgas_par/temp_par - 1d0)
1524 beta = 0.5d0*pi_par*bet1_par*denw_par/bet2_par/alfa_par/ wparc/
1528 cf1 = 0.5d0*(((1.d0/bet2_par)/(alfa_par*wparc))**0.5d0)
1538 if (ntot .gt. zero_par)
then
1539 call sintegral (x1,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
1540 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
1543 y1 = (sinteg1*cf1+sinteg2*cf2)*beta*x1 - 1d0
1546 if (ntot .gt. zero_par)
then
1547 call sintegral (x2,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
1548 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
1551 y2 = (sinteg1*cf1+sinteg2*cf2)*beta*x2 - 1d0
155520
do 30 i=1,maxit_par
1558 if (ntot .gt. zero_par)
then
1559 call sintegral (x3,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
1560 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
1563 y3 = (sinteg1*cf1+sinteg2*cf2)*beta*x3 - 1d0
1565 if (sign(1.d0,y1)*sign(1.d0,y3) .le. zero_par)
then
1574 if (abs(x2-x1) .le. eps_par*x1)
goto 40
1584 if (ntot .gt. zero_par)
then
1585 call sintegral (x2,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
1586 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
1590 y3 = (sinteg1*cf1+sinteg2*cf2)*beta*x3 - 1d0
2074 SUBROUTINE iceparam (sigma_w, denice_ice,ddry_ice,np_ice,
2075 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
2076 & g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,sc_ice,
2077 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,
2078 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2079 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2080 & si0bc_ice,nbc_ice,one_over_tao_preex,
2081 & nhet, nice, smax, nlim,Nhet_dep,Nhet_dhf,fdust_dep
2082 & ,P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2083 & use_av_v,miuv_ice,vpresw_ice,vpresi_ice,denair_ice)
2085 real*8,
intent(in) :: sigma_w,denice_ice,ddry_ice,np_ice,
2086 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
2087 & g1_ice, g2_ice,gdoin_ice, z_ice, lambda_ice,
2088 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,
2089 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2090 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2091 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,p_ice, t_ice,
2092 & vpresw_ice,vpresi_ice,denair_ice
2093 real*8,
dimension(:) :: ndust_ice, sigdust_ice,ddust_ice
2096 real*8,
intent(out) :: nice, nhet, smax, nlim,nhet_dep,nhet_dhf
2099 real*8 :: wpar_ice, sigmav_ice,vmax_ice,miuv_ice
2100 integer,
intent(in)::nbindust_ice
2109 vmax_ice= max(min(miuv_ice+(4d0*sigmav_ice), vmax_ice),
2112 if ((sigmav_ice .lt. 0.05) .or. (t_ice .gt. thom))
then
2116 if (vmax_ice .gt. vmin_ice + 0.01)
then
2118 wpar_ice=min(max(miuv_ice + sigma_w*0.8, vmin_ice), vmax_ice)
2119 CALL nice_param(wpar_ice, denice_ice,ddry_ice,np_ice,nin_ice,
2120 & alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2121 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2122 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2123 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2124 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2125 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
2126 & nice, smax, nhet, nlim,nhet_dep,nhet_dhf,fdust_dep,
2127 & p_ice, t_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2128 & vpresw_ice,vpresi_ice,denair_ice)
2130 call nice_vdist(denice_ice,ddry_ice,np_ice,
2131 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2132 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2133 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2134 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2135 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2136 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
2137 & nice, smax, nhet, nlim,nhet_dep,nhet_dhf,fdust_dep,
2138 & p_ice, t_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice
2139 & ,vpresw_ice,vpresi_ice,denair_ice)
2164 subroutine nice_vdist(denice_ice,ddry_ice,np_ice,
2165 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2166 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2167 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2168 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2169 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2170 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
2171 & nice, smax, nhet, nlim,Nhet_dep,Nhet_dhf,fdust_dep,
2172 & P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2173 & vpresw_ice,vpresi_ice,denair_ice)
2175 real*8 :: quadx(6), wpar, sum1, quadw(6), dp, sum2, sum3, sum4,
2177 real*8,
intent(out) :: nice, smax, nhet, nlim,nhet_dep,nhet_dhf
2179 real*8 wpar_icex,denice_ice,ddry_ice,np_ice,
2180 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2181 & normv_ice,sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2182 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2183 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,sc_ice,
2184 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2185 & si0bc_ice,nbc_ice,one_over_tao_preex,p_ice, t_ice,
2186 & vpresw_ice,vpresi_ice,denair_ice
2187 real*8,
dimension(:) :: ndust_ice, sigdust_ice,ddust_ice
2189 INTEGER :: INDEX,nbindust_ice
2191 DATA quadx/0.23861918d0, -0.23861918d0, 0.66120939d0, -
2192 &0.66120939d0, 0.93246951d0, -0.93246951d0/
2194 DATA quadw/0.46791393d0, 0.46791393d0, 0.36076157d0,
2195 & 0.36076157d0, 0.17132449d0, 0.17132449d0/
2199 x1=(vmin_ice-miuv_ice)/(sq2_par*sigmav_ice)
2200 x2=(vmax_ice-miuv_ice)/(sq2_par*sigmav_ice)
2202 x2=max(x1 +0.01, x2)
2203 normv_ice=(erfapp(x2)-erfapp(x1))*0.5d0
2212 wpar=max(0.5d0*(((vmax_ice-vmin_ice)*quadx(index)) +(vmax_ice+
2216 CALL nice_param(wpar, denice_ice,ddry_ice,np_ice,nin_ice,
2217 & alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2218 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2219 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2220 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2221 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2222 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
2223 & nice, smax, nhet, nlim,nhet_dep,nhet_dhf,fdust_dep,
2224 & p_ice, t_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2225 & vpresw_ice,vpresi_ice,denair_ice)
2227 CALL gausspdf(wpar, dp, sigmav_ice,miuv_ice,normv_ice)
2228 sum1=sum1+(nice*dp*quadw(index))
2229 sum2=sum2+(smax*dp*quadw(index))
2230 sum3=sum3+(nhet*dp*quadw(index))
2231 sum4=sum4+(nlim*dp*quadw(index))
2232 sum5=sum5+(sc_ice*dp*quadw(index))
2237 nice=sum1*(vmax_ice-vmin_ice)*0.5d0
2238 smax=sum2*(vmax_ice-vmin_ice)*0.5d0
2239 nhet=sum3*(vmax_ice-vmin_ice)*0.5d0
2240 nlim=sum4*(vmax_ice-vmin_ice)*0.5d0
2241 sc_ice=sum5*(vmax_ice-vmin_ice)*0.5d0
2301 subroutine nice_param(wpar_icex,denice_ice,ddry_ice,np_ice,
2302 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2303 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2304 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2305 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2306 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2307 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
2308 & nice, smax, nhet, nlim_,Nhet_dep,Nhet_dhf,fdust_dep,
2309 & P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2310 & vpresw_ice,vpresi_ice,denair_ice)
2312 real*8,
intent(in) :: wpar_icex,denice_ice,np_ice,nin_ice,
2313 & alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2314 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2315 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,
2316 & kdust_ice, kbc_ice, shdust_ice, shbc_ice, ddry_ice,
2317 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2318 & si0bc_ice,nbc_ice,one_over_tao_preex,p_ice, t_ice,
2319 & vpresw_ice,vpresi_ice,denair_ice
2320 real*8,
intent(out) :: nice, smax, nhet, nlim_,nhet_dep,nhet_dhf
2323 real*8 :: aux1, aux2, g, dpmax, miu, monvol, fds, nlim, dlim,
2324 & dstar, ds, nstar, nhom, fc, phido, auxnc, sizecorr, dsh, nhet_,
2325 & f1, f2, saux, sup, slow, shalf, fhalf, dpmin, gam, wpar_ice,
2326 & preex_effect, swsat,
sstep,sc_ice,lambda_ice
2327 real*8,
dimension(:) :: ndust_ice, sigdust_ice,ddust_ice
2329 integer :: INDEX, maxiter_s,nbindust_ice
2335 preex_effect=1.0-(one_over_tao_preex*(shom_ice)/alfa_ice/
2336 & ((shom_ice+1.0))/wpar_icex)
2337 swsat = vpresw_ice/vpresi_ice -1d0
2340 if (preex_effect .le. 0.0)
then
2347 sc_ice = shom_ice + 1.d0
2353 wpar_ice = wpar_icex*preex_effect
2359 if (np_ice .gt. 1.0)
then
2361 monvol=np_ice*1.0d-6*ddry_ice*ddry_ice*ddry_ice
2362 aux1=1.6397d-14*t_ice-3.1769d-12
2363 dpmax=aux1*(monvol**(-0.373d0))*(wpar_ice**(-0.05))
2364 IF (dpmax.GT.1.0d-4)
THEN
2374 dpmin=dliq_ice+(0.02/sqrt(alfa_ice*wpar_ice*g1_ice))
2375 dpmax=max(dpmax,dpmin)
2378 aux2=dlog((g2_ice+(g1_ice*dpmax))/(g2_ice+(g1_ice*dliq_ice)))
2379 g=1.3346d0*((g1_ice*aux1)-(g2_ice*aux2))/(aux1*g1_ice*g1_ice)
2380 lambda_ice=lambda_ice/sqrt(wpar_ice)
2381 aux1 = g1_ice*alfa_ice*wpar_ice
2382 nstar=(aux1*sqrt(aux1))/beta_ice/z_ice/sq2_par
2394 if (typeofspec_ice .ge. 0d0)
then
2396 call inspec_ice(shom_ice,nhet_,dsh,np_ice,norg_ice,sigorg_ice,
2397 & dorg_ice, dbc_ice,sigbc_ice,
2398 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2399 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2400 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2401 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice
2402 & ,vpresw_ice,vpresi_ice)
2403 sizecorr=exp(-2d0/lambda_ice/shom_ice)
2404 dstar=((4d0*dsh*dsh/3d0)+(2d0*dsh*(shom_ice-dsh))) /(shom_ice-
2407 if (dstar .gt. 0.0)
then
2408 nlim=min(nstar*(shom_ice+1d0)/shom_ice/sqrt(dstar)/sizecorr,
2416 dstar=((4d0*dsh*dsh/3d0)+(2d0*dsh*(shom_ice-dsh))) /(shom_ice-
2418 dlim=-gam+sqrt((gam*gam)+(2d0*dstar/g1_ice/alfa_ice/wpar_ice))
2419 nlim=alfa_ice*wpar_ice*(shom_ice+1d0)/z_ice/beta_ice/shom_ice
2420 nlim=nlim*((g1_ice*dlim)+g2_ice)/dlim/dlim
2425 nlim_=min(nlim, 1d10)
2426 nlim_=max(nlim, 1d-6)
2428 if (nhet_ .gt. 0.0)
then
2430 fds=1d0-(aux1*sqrt(aux1))
2437 if (purehom_ice)
then
2441 if ((purehet_ice) .or. (t_ice .GE. thom))
then
2446 IF (fds .GE. 1.e-10)
THEN
2448 miu=fds*alfa_ice*(shom_ice+1d0)*wpar_ice*koft_ice/shom_ice
2450 phido=sqrt(pi_ice*g/miu/2d0)*(g/miu)
2451 auxnc=2d0*denair_ice/beta_ice/koft_ice/denice_ice/pi_ice/np_ice
2456 if (np_ice .gt.0d0)
then
2458 IF (fc .le. 0.6d0)
then
2459 nhom=np_ice*exp(-fc)*(1.0d0-exp(-fc))
2461 nhom=np_ice/(1d0+exp((9d0-2d0*fc)/7d0))
2470 if (purehom_ice) nhet_ = 0.d0
2483 if (typeofspec_ice .lt. 0d0) saux=sh_ice+0.00000000001d0
2485 f1= findsmax(saux,dsh,
2486 & nhet_,np_ice,norg_ice,sigorg_ice,
2487 & dorg_ice, dbc_ice,sigbc_ice,
2488 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2489 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2490 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2491 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
2492 & gdoin_ice,alfa_ice,wpar_ice,gam,g1_ice,g2_ice,z_ice,
2493 & beta_ice,nin_ice,vpresw_ice,vpresi_ice,nstar)
2496 maxiter_s = int(1d0/
sstep) + 1
2498 do index =1, maxiter_s
2500 if (saux .ge. swsat)
then
2508 f2=findsmax(saux,dsh,
2509 & nhet_,np_ice,norg_ice,sigorg_ice,
2510 & dorg_ice, dbc_ice,sigbc_ice,
2511 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2512 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2513 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2514 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
2515 & gdoin_ice,alfa_ice,wpar_ice,gam,g1_ice,g2_ice,z_ice,
2516 & beta_ice,nin_ice,vpresw_ice,vpresi_ice,nstar)
2517 IF (f2*f1 .lt. 0d0)
exit
2521 if (f2*f1 .gt.0d0)
then
2526 if (saux .lt. swsat)
then
2530 slow=saux-(
sstep + 0.001d0)
2533 shalf=0.5d0*(sup+slow)
2534 fhalf=findsmax(shalf,dsh,
2535 & nhet_,np_ice,norg_ice,sigorg_ice,
2536 & dorg_ice, dbc_ice,sigbc_ice,
2537 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2538 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2539 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2540 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
2541 & gdoin_ice,alfa_ice,wpar_ice,gam,g1_ice,g2_ice,z_ice,
2542 & beta_ice,nin_ice,vpresw_ice,vpresi_ice,nstar)
2544 IF (sign(1.d0,f1)*sign(1.d0,fhalf) .LE. 0d0)
THEN
2553 IF (abs(slow-sup) .LE. 5d-3)
exit
2563 if ((smax .gt. shom_ice) .and.(t_ice .LE. thom)) smax=shom_ice
2565 if (typeofspec_ice .ge. 0d0)
then
2566 call inspec_ice(smax, nhet,dsh,np_ice,norg_ice,sigorg_ice,
2567 & dorg_ice, dbc_ice,sigbc_ice,
2568 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2569 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2570 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2571 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2572 & vpresw_ice,vpresi_ice)
2581 sc_ice=max(smax+1.0-dsh, 1.0)
2586 if (fds .gt. 0.0)
then
2587 sc_ice = (shom_ice+1.0)*fds + sc_ice*(1.0-fds)
2591 sc_ice=min(shom_ice+1.0, sc_ice)
2598 real*8 function findsmax(SX,DSH,
2599 & NHET_,np_ice,norg_ice,sigorg_ice,
2600 & dorg_ice, dbc_ice,sigbc_ice,
2601 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2602 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2603 & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
2604 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
2605 & gdoin_ice,alfa_ice,wpar_ice,GAM,g1_ice,g2_ice,z_ice,
2606 & beta_ice,nin_ice,vpresw_ice,vpresi_ice,NSTAR)
2607 real*8,
intent(in) :: sx
2609 real*8 :: dstar, sizecorr, dsh, nstar,dlim
2610 integer nbindust_ice
2611 real*8,
dimension(:) :: ndust_ice, sigdust_ice,ddust_ice
2612 real*8 shom_ice,nhet_,np_ice,norg_ice,sigorg_ice,
2613 & dorg_ice, dbc_ice,sigbc_ice,
2614 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2615 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2616 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2618 & gdoin_ice,alfa_ice,wpar_ice,gam,g1_ice,g2_ice,z_ice,
2619 & beta_ice,nin_ice,vpresw_ice,vpresi_ice
2621 if (typeofspec_ice .ge. 0d0)
then
2623 call inspec_ice(sx,nhet_,dsh,np_ice,norg_ice,sigorg_ice,
2624 & dorg_ice, dbc_ice,sigbc_ice,
2625 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2626 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2627 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2628 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2629 & vpresw_ice,vpresi_ice)
2630 sizecorr=exp(-2d0/lambda_ice/sx)
2631 dstar=((4d0*dsh*dsh/3d0)+(2d0*dsh*(sx-dsh)))/(sx-dsh+1d0)
2632 dstar=dstar+(gdoin_ice*alfa_ice*wpar_ice)
2635 tao=nhet_*sizecorr*sx*sqrt(dstar)/(sx+1d0)/nstar
2640 dstar=((4d0*dsh*dsh/3d0)+(2d0*dsh*(sx-dsh))) /(sx-dsh+1d0)
2641 dlim=-gam+sqrt((gam*gam)+(2d0*dstar/g1_ice/alfa_ice/wpar_ice))
2642 tao=alfa_ice*wpar_ice*(sx+1d0)/z_ice/beta_ice/sx
2643 tao=tao*((g1_ice*dlim)+g2_ice)/dlim/dlim/nin_ice
2781 SUBROUTINE prop_ice(T, P, denice_ice,ddry_ice,
2782 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
2783 & g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,
2784 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2785 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2786 & si0bc_ice,nbc_ice,D_preex, N_preex,one_over_tao_preex,
2787 & P_ice, T_ice,act_param,ndust_ice,vpresw_ice,vpresi_ice,
2790 real*8,
intent(in) :: t, p,ddry_ice,nbc_ice,d_preex, n_preex
2791 real*8,
intent(out) :: alfa_ice,beta_ice,shom_ice,
2792 & koft_ice, dliq_ice,g1_ice, g2_ice,gdoin_ice, z_ice,
2793 & lambda_ice,kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2794 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2795 & si0bc_ice,one_over_tao_preex,p_ice, t_ice
2797 real*8 :: aux, aux1, aux2, sw, fice, mice, tc, hdust, hbc, b0,
2798 & b1, b2, b3, x, t0bc, t0dust, gam, gamma,nin_ice,
2799 & tr, vw, den_m, eact, toact, dact, acc, n1, siw, rgo,
2800 & ngo, mw_molec,dhs_ice,denwat_ice,denice_ice,vpresw_ice,
2801 & vpresi_ice,denair_ice,diff_ice,aircond_ice
2803 real*8,
dimension(:) :: ndust_ice
2804 integer,
intent(in) :: act_param
2807 t_ice = min(max(t, tmin_ice), to_ice)
2808 p_ice = max(p, pmin_ice)
2813 dhs_ice = dhsub_ice(t_ice)
2814 vpresw_ice = vpreswater_ice(t_ice)
2815 vpresi_ice = vpresice(t_ice)
2816 denice_ice = densityice(t_ice)
2817 denwat_ice = watdensity_ice(t_ice)
2818 denair_ice = p_ice*amw_ice/rgas_ice/t_ice
2821 diff_ice=(0.211d0*101325d0/p_ice)*((t_ice/273d0)**1.94d0)*1.0d-4
2822 aux1=1.0e-3*(4.39d0+0.071d0*t_ice)
2825 aux2=(2d0*aux1/(thaccom_ice*1.0d-6*denair_ice*cpa_ice)) *((58.0d-
2826 &3*pi_ice/(rgas_ice*t_ice))**0.5d0)
2828 aircond_ice=aux1/(1.d0+aux2)
2833 aux1=grav_ice*dhs_ice/rv_ice/t_ice/t_ice/cpa_ice
2834 aux2=grav_ice*amw_ice/rgas_ice/t_ice
2836 beta_ice=amw_ice*p_ice/wmw_ice/vpresi_ice
2837 gamma=1.5d0*dhs_ice*dhs_ice/rv_ice/t_ice/t_ice/cpa_ice
2839 beta_ice=beta_ice+gamma
2843 shom_ice=2.349d0-(t_ice/259d0)
2844 sw=shom_ice*vpresi_ice/vpresw_ice
2845 shom_ice=shom_ice-1d0
2846 koft_ice=(0.0240d0*t_ice*t_ice)-(8.035d0*t_ice)+934.0d0
2852 if (sw .lt. 0.99)
then
2853 aux1=(1d0/(1d0-sw))-1.1764d0
2855 aux1=(1d0/0.01)-1.1764d0
2857 dliq_ice=ddry_ice*0.9344d0*(aux1**0.333)
2861 aux1=denice_ice*rv_ice*t_ice/vpresi_ice/diff_ice
2862 aux2=dhs_ice*denice_ice/aircond_ice/t_ice
2863 aux2=aux2*((dhs_ice/rv_ice/t_ice)-1.0d0)
2864 g1_ice=(aux1+aux2)/4.0d0
2867 one_over_tao_preex = beta_ice*denice_ice*pi_ice*0.5* d_preex*
2868 &n_preex/g1_ice/denair_ice
2872 g2_ice=denice_ice*rv_ice*t_ice/2.0d0/vpresi_ice/depcoef_ice
2873 g2_ice=g2_ice*((2.0d0*pi_ice/rv_ice/t_ice)**0.5d0)
2876 gdoin_ice=(g1_ice*0.5d0*doin_ice*doin_ice)+(g2_ice*doin_ice)
2877 z_ice=denice_ice*pi_ice/2.0d0/denair_ice
2880 lambda_ice=1d0/sqrt(alfa_ice*g1_ice*gam*gam)
2886 if (typeofspec_ice .lt. 0)
then
2888 nin_ice=(sum(ndust_ice)+nbc_ice)*0.05d0
2889 elseif (typeofspec_ice .eq. 3)
then
2896 fice=0.25d0*((mice*mice*mice)-(3d0*mice)+2d0)
2897 kdust_ice=koft_ice*fice
2899 fice=0.25d0*((mice*mice*mice)-(3d0*mice)+2d0)
2900 kbc_ice=koft_ice*fice
2901 elseif (typeofspec_ice .eq. 4)
then
2906 b0=-1.0261d0; b1=3.1656d-3; b2=5.3938d-4; b3=8.2584d-6
2907 x=b0+(b1*tc)+(b2*tc*tc)+(b3*tc*tc*tc)
2908 si0dust_ice=1d0+(10d0**x)
2909 del1dust_ice=cubicint_ice(tc, t0dust, t0dust+5d0, 1d0, hdust)
2912 b0=0.5652d0; b1=1.085d-2; b2=-3.118d-5
2913 si0bc_ice=b0+(b1*t_ice)+(b2*t_ice*t_ice)-0.1d0
2914 del1bc_ice=cubicint_ice(tc, t0bc, t0bc+5d0, 1d0, hbc)
3057 subroutine inspec_ice(six, N, Dsh,np_ice,norg_ice, sigorg_ice,
3058 & dorg_ice, dbc_ice,sigbc_ice,
3059 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
3060 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
3061 & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
3062 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
3063 & vpresw_ice,vpresi_ice)
3065 real*8,
intent(in) :: six,np_ice,norg_ice, sigorg_ice,
3066 & dorg_ice, dbc_ice,sigbc_ice,
3067 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
3068 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
3069 & si0bc_ice,nbc_ice,p_ice, t_ice
3070 real*8,
intent(out) :: n, dsh,nhet_dep,nhet_dhf,fdust_dep
3071 real*8 :: nd, nbc, aux, si_, sw, del0, ddel0, fc, delw0, ddelw0,
3072 & sw0, hdust, hbc, nbase, dnd, dnbc, dnbase, dh, dfc, ndaux,
3073 & dndaux, dnorg, norg, ndustaux, frac, aux2, dx2, fdep, ndep, ndhf,
3074 & dndep, dndhf, si, dfrac, ncdep_, ncdhf_ , fglassy, nglassy,
3075 & dnglassy, siw, d_grid_bio, n_grid_bio,vpresw_ice,vpresi_ice
3077 real*8,
dimension(3) :: sig_array, the_array, frac_array
3078 real*8,
dimension(1:nbindust_ice) :: ndust_ice, sigdust_ice,
3081 real :: n_iw, DSh_s , nbc_s, dbc_s, Asolo
3082 real,
dimension (nbindust_ice) :: ndust_s, ddust_s
3083 integer :: index, kindex, mode,nbindust_ice
3087 sw=si_*vpresi_ice/vpresw_ice
3098 if ((six .lt. 0.02) .or. (t_ice .gt. 270.0))
then
3104 if (sw .ge. 1.0)
then
3106 si_=vpresw_ice/vpresi_ice
3109 siw=vpresw_ice/vpresi_ice
3116 select case (typeofspec_ice)
3119 n=1000d0*exp(-0.639d0)*(exp(12.96d0*si)-1d0)
3123 n=pdg07_ice(si, t_ice)
3124 if (t_ice .le. 243d0)
then
3133 DO index=1, nbindust_ice
3135 ndustaux=ndustaux+ndust_ice(index)
3139 if (si .le.shdust_ice)
then
3140 nd=(si/shdust_ice)*ndustaux*effdust_ice* exp(-kdust_ice*
3142 dnd=nd*((1d0/si)+kdust_ice)
3145 nd=ndustaux*effdust_ice
3150 if (si .le.shbc_ice)
then
3151 nbc=(si/shbc_ice)*nbc_ice*effbc_ice* exp(-kbc_ice*(shbc_ice-si))
3152 dnbc=nbc*((1d0/si)+kbc_ice)
3154 nbc=nbc_ice*effbc_ice
3159 if (((dnd+dnbc) .gt. 0d0) .and. (n .gt. 0.0))
then
3170 delw0=cubicint_ice(sw, sw0, 1d0, 0d0, 1d0)
3171 ddelw0=dcubicint_ice(sw, sw0, 1d0, 0d0, 1d0)
3173 nbase=pdg07_ice(si, t_ice)
3176 if (t_ice .le. 243d0)
then
3179 dnbase=12.96d0*nbase
3183 del0=cubicint_ice(si_, si0dust_ice, si0dust_ice+0.1d0, 0d0, 1d0)
3184 ddel0=dcubicint_ice(si_, si0dust_ice, si0dust_ice+0.1d0, 0d0,
3187 fc=0.5d0*del1dust_ice*del0
3188 dfc=0.5d0*del1dust_ice*ddel0
3190 hdust=fc+((1d0-fc)*delw0)
3191 dh=(dfc*(1d0-delw0))+(ddelw0*(1d0-fc))
3193 if (hdust .gt. 1d0)
then
3199 aux=(2d0/3d0)*hdust*(nbase/0.76d0)*pi_ice/5.0d-7/4d0
3200 aux2=(2d0/3d0)*pi_ice/0.76d0/5.0d-7/4d0
3205 DO index =1, nbindust_ice
3208 dx2=ddust_ice(index)*ddust_ice(index)*ddust_ice(index)*0.52*
3211 frac=0.5d0*(1d0-erfapp(-log(ddust_ice(index)/0.1e-6) /
3212 &sigdust_ice(index)/sq2_par))
3214 ndaux=frac*ndust_ice(index)*(1d0-exp(-aux*dx2))
3217 ndaux=(frac*ndust_ice(index)-ndaux)
3218 dndaux=ndaux*((dh*nbase)+(hdust*dnbase))*aux2*dx2
3227 del0=cubicint_ice(si_, si0bc_ice, si0bc_ice+0.1d0, 0d0, 1d0)
3228 ddel0=dcubicint_ice(si_, si0bc_ice, si0bc_ice+0.1d0, 0d0, 1d0)
3230 fc=0.5d0*del1bc_ice*del0
3231 hbc=fc+((1d0-fc)*delw0)
3232 dfc=0.5d0*del1bc_ice*ddel0
3233 dh=(dfc*(1d0-delw0))+(ddelw0*(1d0-fc))
3236 if (hbc .gt. 1d0)
then
3241 frac=0.5d0*(1d0 -erfapp(-log(dbc_ice/0.1e-6) /sigbc_ice/
3245 dx2=dbc_ice*dbc_ice*dbc_ice*0.52*acorr_bc
3247 aux=((1d0/3d0)-0.06d0)*hbc*(nbase/0.76d0)*pi_ice/2.7d-7
3248 aux2=((1d0/3d0)-0.06d0)*pi_ice/0.76d0/2.7d-7
3250 nbc=nbc_ice*frac*(1d0-exp(-aux*dx2))
3251 dnbc=(nbc_ice*frac-nbc)*((dh*nbase)+(hbc*dnbase))*aux2*dx2
3256 frac=0.5d0*(1d0-erfapp(-log(dorg_ice/0.1e-6) /sigorg_ice/
3259 dx2=dorg_ice*dorg_ice
3261 aux=0.06d0*hbc*(nbase/0.76d0)*pi_ice/9.1d-7
3262 aux2=0.06d0*pi_ice/0.76d0/9.1d-7
3265 norg=norg_ice*frac*(1d0-exp(-aux*dx2))
3266 dnorg=(norg_ice*frac-norg)*((dh*nbase)+(hbc*dnbase))*aux2*dx2
3279 if (t_ice .lt. 210.0)
then
3280 nglassy= min(0.01 +0.0045*(210.0 -t_ice), 0.1)
3281 fglassy= 7.7211*1e-3*si_ - 9.2688*1e-3
3282 fglassy=min(fglassy, 3.3587e-3)
3283 fglassy=max(fglassy, 0.0)
3288 nglassy = np_ice*nglassy*fglassy
3289 dnglassy =nglassy*7.7211*1e-3
3296 if ((dnd+dnbc+dnorg+dnglassy) .gt. 0d0)
then
3297 dsh=n/(dnd+dnbc+dnorg+dnglassy)
3307 d_grid_bio =dorg_ice
3314 ndust_s = sngl(frac_dust*ndust_ice)
3315 nbc_s = sngl(frac_bc*nbc_ice)
3316 asolo = sngl(0.25d0*frac_org*norg_ice*pi_par*dorg_ice*dorg_ice)
3320 DO index =1, nbindust_ice
3321 frac=0.5d0*(1d0-erfapp(log(0.1e-6/ddust_ice(index)) /
3322 &sigdust_ice(index)/sq2_par))
3323 ndust_s(index) = sngl(frac*ndust_ice(index))
3326 frac=0.5d0*(1d0-erfapp(log(0.1e-6/dbc_ice) /sigbc_ice/sq2_par))
3328 nbc_s = sngl(frac*nbc_ice)
3330 frac=0.5d0*(1d0-erfapp(log(0.1e-6/dorg_ice) /sigorg_ice/
3333 asolo = sngl(frac*norg_ice*pi_par*dorg_ice*dorg_ice)
3336 ddust_s=sngl(ddust_ice)
3337 dbc_s = sngl(dbc_ice*1.0)
3339 call empirical_param_phillips(sngl(si_), sngl(siw), sngl(sw), (/
3340 &ddust_s/), (/ndust_s/), 5, (/dbc_s/), (/nbc_s/), 1, (/
3341 &sngl(d_grid_bio)/), (/sngl(n_grid_bio)/), 1, asolo, n_iw, dsh_s
3342 & ,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice)
3352 if (dsh .ge. si)
then
3356 if (t_ice .gt. 255.0)
then
3475 SUBROUTINE empirical_param_phillips(SI, SIW, SW, D_grid_dust,
3476 & n_grid_dust, ijstop_dust, D_grid_soot, n_grid_soot, ijstop_soot,
3477 & D_grid_bio, n_grid_bio, ijstop_bio, A_solo, n_iw, DSH,
3478 & Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice)
3480 real,
intent(IN):: SI, SIW, SW, A_solo
3481 real*8,
intent(IN):: p_ice, t_ice
3482 real,
dimension(:),
intent(IN):: D_grid_dust, n_grid_dust,
3483 & d_grid_soot, n_grid_soot, d_grid_bio, n_grid_bio
3484 integer,
intent(IN):: ijstop_dust, ijstop_soot, ijstop_bio
3487 real :: nin_a_nuc_dust, nin_a_nuc_soot, nin_a_nuc_bio,
3488 & nin_a_nuc_solo, num_ic_dust_imm, num_ic_soot_imm, num_ic_bio_imm,
3491 real,
intent (inout) :: DSH, n_iw
3492 real*8,
intent (out) :: nhet_dep,nhet_dhf,fdust_dep
3494 real :: dn_in_dust, dn_in_soot, dn_in_bio, dn_in_solo, dNall,
3495 & dnaux, naux, ss_w, dh_frac_dust, dh_frac_soot, dh_frac_solo, aux,
3496 & dfdep, temperature_k, p_sat, ahet_aux
3500 REAL :: RHO_CFDC, BASE_DUST_OMEGA, BASE_SOOT_PHILIC_OMEGA,
3501 & base_bio_omega, alpha_dust, alpha_soot, alpha_bio,
3502 & fraction_depnucl_warm_dust, pie, base_solo_omega,
3503 & temp_max_dust_degc, temp_max_soot_degc, temp_max_bio_degc,
3505 parameter( base_dust_omega = 2.0e-6, base_soot_philic_omega = 1.e-
3506 &7, base_bio_omega = 0.89e-6, base_solo_omega = 5.6e-5,
3507 & glass_frac = 0.5, alpha_dust = 2./3., alpha_soot = 1./3. - 0.03,
3508 & alpha_bio = 0.03, rho_cfdc = 50000./(287.*228.15),
3509 & fraction_depnucl_warm_dust = 0.15, pie = 3.1415926,
3510 & temp_max_dust_degc = -10., temp_max_soot_degc = -15.,
3511 & temp_max_bio_degc = -2.)
3513 real :: FAC_CORRECT_RH = 2., rho_aida
3514 real:: H_frac_dust, n_in, n_in_dust, n_in_ultra, n_in_dust_ultra,
3515 & cihenc_dust, esw, esi, ss_i, n_in_soot_ultra, h_frac_soot,
3516 & h_frac_bio, n_in_soot, n_in_bio, n_in_bio_ultra, cihenc_soot,
3517 & cihenc_bio, delta_si, delta_t, delta_sw, n_in_max, ss_iw, rho
3519 real :: H_frac_solO, RHI, n_in_solO, n_in_solO_star, CIHENC_solO,
3522 real :: mu, S_i_0, RH_crit, S_i_w_warm, S_i_w_cold, S_i_w,
3524 real :: S_w_0, dep_frac, n_in_hat, n_in_tilde
3533 if(ijstop_dust .ne.
SIZE(n_grid_dust)) stop 6366
3534 if(ijstop_soot .ne.
SIZE(n_grid_soot)) stop 6366
3535 if(ijstop_bio .ne.
SIZE(n_grid_bio)) stop 6366
3547 nin_a_nuc_dust=0.0; nin_a_nuc_soot=0.0; nin_a_nuc_bio=
3548 &0.0; nin_a_nuc_solo=0.0
3549 num_ic_dust_imm=0.0; num_ic_soot_imm=0.0; num_ic_bio_imm=
3550 &0.0; num_ic_solo_imm=0.0
3552 n_in_dust=0.0; dn_in_dust=0.0;n_in_soot=0.0
3553 dn_in_soot=0.0; dn_in_bio=0.0; n_in_bio=0.0;
3554 dn_in_solo=0.0;n_in_solo=0.0;n_in_dust_ultra=0.0
3555 n_in_soot_ultra=0.0;n_in_bio_ultra=0.0;
3564 temperature_k=sngl(t_ice)
3575 rho_aida = 90000./(287.*205.)
3577 rho = p_sat/(287.*temperature_k)
3579 psi_solo = a_solo/base_solo_omega
3580 ss_i = min(max(si-1.0, 0.0), 1.0)
3581 ss_w = min(max(sw-1.0, -1.0), 1.0)
3582 ss_iw = min(max(siw - 1.0, 0.0), 1.0)
3586 if(temperature_k < 273.15 .and. temperature_k > 273.15 -
3599 delta_si = h_1_smooth(ss_i + 1, 1.1, 1.2, 0.0, 1.,dh1smooth);
3600 delta_t = h_1_smooth(-(temperature_k-273.15), 35., 40.,
3601 & fraction_depnucl_warm_dust, 1.,dh1smooth);
3602 delta_sw = h_1_smooth(ss_w + 1.0, 0.97, 1., 0., 1.,dh1smooth);
3604 tc_hm_degc = temperature_k - 273.15
3607 s_i_0 = 1. + 10.**(8.2584e-6*tc_hm_degc*tc_hm_degc*tc_hm_degc +
3608 & 5.3938e-4*tc_hm_degc*tc_hm_degc + 3.1656e-3*tc_hm_degc - 1.0261)
3613 aux =h_1_smooth(-(temperature_k-273.15), 35., 40.,
3614 & fraction_depnucl_warm_dust, 1.,dh1smooth)/fac_correct_rh
3616 dep_frac = h_1_smooth(ss_i + 1, s_i_0, s_i_0 + 0.1, 0.,1.,
3620 aux= h_1_smooth(ss_w + 1.0, s_w_0, 1., 0.,1.,dh1smooth)
3622 h_frac_dust = dep_frac + (1. - dep_frac)*aux
3624 dh_frac_dust = dfdep + (siw*(1. - dep_frac)*dh1smooth)- aux*
3627 if(h_frac_dust > 1.) h_frac_dust = 1.
3629 if ((h_frac_dust .gt. 1.0e-6) .and. (h_frac_dust .lt. 1.))
then
3630 dh_frac_dust = dh_frac_dust/h_frac_dust
3638 aux =h_1_smooth(-(temperature_k-273.15), 65., 75., 0.,1.,
3640 dep_frac = h_1_smooth(ss_i + 1, s_i_0, s_i_0+0.1, 0.,1.,
3642 h_frac_solo = dep_frac
3644 if ((h_frac_solo .gt. 1.0e-6) .and. (h_frac_solo .lt. 1.))
then
3645 dh_frac_solo = dh1smooth/h_frac_solo
3647 if(h_frac_solo > 1.) h_frac_solo = 1.
3654 aux = h_1_smooth(-(temperature_k-273.15), 40., 50., 0.,1.,
3655 & dh1smooth) / fac_correct_rh
3656 dep_frac = h_1_smooth(ss_i + 1, s_i_0, s_i_0+0.1, 0.,1.,
3659 dfdep= dh1smooth*aux
3661 aux = h_1_smooth(ss_w + 1.0, s_w_0, 1., 0.,1.,dh1smooth)
3662 h_frac_soot = dep_frac + (1. - dep_frac)*aux
3663 if(h_frac_soot > 1.) h_frac_soot = 1.
3665 dh_frac_soot = dfdep + (siw*(1. - dep_frac)*dh1smooth)- aux*
3667 if ((h_frac_soot .gt. 1.0e-6) .and. (h_frac_soot .lt. 1.))
then
3668 dh_frac_soot = dh_frac_soot/h_frac_soot
3674 h_frac_bio = h_frac_soot
3676 if(temperature_k < 273.15 .and. temperature_k >= 273.15 -
3678 n_in = 1.e3* (exp(12.96*ss_i - 0.639)/rho_cfdc) *0.0587*
3680 if( temperature_k > 273.15 -5. .and. temperature_k < 273.15 -
3682 n_in = n_in*h_1_smooth(-(temperature_k-273.15), 2., 5., 0., 1.,
3685 if(temperature_k >= 273.15 - 2. ) n_in = 0.
3688 if(temperature_k < 273.15 -25. )
then
3689 n_in_tilde = 1000.*(exp(0.1296*(ss_i*100.-10.))**0.3)*
3690 &fac_correct_rh/rho_cfdc
3693 if(temperature_k >= 273.15 - 30.) n_in_max = 1.e3* (exp(12.96*
3694 &ss_iw - 0.639)/rho_cfdc) *0.0587*fac_correct_rh
3695 if(temperature_k < 273.15 - 30.) n_in_max = 1000.*(exp(0.1296*
3696 &(ss_iw*100.-10.))**0.3)*fac_correct_rh/rho_cfdc
3698 if(n_in_hat > n_in_max) n_in_hat = n_in_max
3699 if(n_in_tilde > n_in_max) n_in_tilde = n_in_max
3703 n_in = n_in_hat * ((n_in_tilde/n_in_hat)**(h_1_smooth(-
3704 &(temperature_k-273.15), 25., 35., 0., 1.,dh1smooth)))
3707 if(n_in > n_in_max) n_in = n_in_max
3714 if(temperature_k < 273.15 - 30.)
then
3722 do ij = 1, ijstop_dust
3724 mu = n_in*alpha_dust*h_frac_dust*ahet_dust(ij)/base_dust_omega
3726 ahet_aux = pie*d_grid_dust(ij)*d_grid_dust(ij)*d_grid_dust(ij)*
3727 &4.73*acorr_dust/6.0
3728 mu = n_in*alpha_dust*h_frac_dust*ahet_aux/base_dust_omega
3730 naux = (1. - exp(-mu))*n_grid_dust(ij)
3731 n_in_dust = n_in_dust + naux
3732 dn_in_dust = max(mu*(n_grid_dust(ij)-naux)*(dnaux +
3733 & dh_frac_dust), 0.0) + dn_in_dust
3737 if( temperature_k > 273.15 +temp_max_dust_degc -
3738 & 20. .and. temperature_k < 273.15 + temp_max_dust_degc)
then
3739 n_in_dust = n_in_dust*h_1_smooth(-(temperature_k-273.15),-
3740 &temp_max_dust_degc,-temp_max_dust_degc+20., 0., 1.,dh1smooth)
3742 if(temperature_k >= 273.15 + temp_max_dust_degc) n_in_dust = 0.
3747 do ij = 1, ijstop_soot
3750 mu = n_in*alpha_soot*h_frac_soot*ahet_bc/base_soot_philic_omega
3752 ahet_aux = pie*d_grid_soot(ij)*d_grid_soot(ij)*d_grid_soot(ij)*
3754 mu = n_in*alpha_soot*h_frac_soot*ahet_aux/base_soot_philic_omega
3756 naux = (1. - exp(-mu))*n_grid_soot(ij)
3757 n_in_soot = n_in_soot + naux
3758 dn_in_soot = max(mu*(n_grid_soot(ij)-naux)*(dnaux+dh_frac_soot),
3763 if( temperature_k > 273.15 + temp_max_soot_degc -
3764 & 10. .and. temperature_k < 273.15 + temp_max_soot_degc)
then
3765 n_in_soot = n_in_soot*h_1_smooth(-(temperature_k-273.15),-
3766 &temp_max_soot_degc,-temp_max_soot_degc+10., 0., 1.,dh1smooth)
3769 if(temperature_k >= 273.15 + temp_max_soot_degc) n_in_soot = 0.
3775 do ij = 1, ijstop_bio
3776 mu = n_in*alpha_bio*h_frac_bio*pie*(d_grid_bio(ij)**2.) /
3779 mu = n_in*alpha_bio*h_frac_bio
3780 naux = (1. - exp(-mu))*n_grid_bio(ij)
3785 n_in_bio = n_in_bio + naux
3786 dn_in_bio = max(mu*(n_grid_bio(ij)-naux)*dnaux, 0.0) + dn_in_bio
3792 if( temperature_k > 273.15 + temp_max_bio_degc -
3793 & 3. .and. temperature_k < 273.15 + temp_max_bio_degc)
then
3794 n_in_bio = n_in_bio*h_1_smooth(-(temperature_k-273.15),-
3795 &temp_max_bio_degc,-temp_max_bio_degc+3., 0., 1.,dh1smooth)
3798 if(temperature_k >= 273.15 + temp_max_bio_degc ) n_in_bio = 0.
3803 n_in = 0.; n_in_ultra = 0.; n_in_dust = 0.; n_in_soot =
3804 & 0.; n_in_bio = 0.;
3807 if(temperature_k < 273.15 - 35.)
then
3808 n_in_ultra = 1000.*(exp(0.1296*(ss_i*100.-10.))**0.3)*
3809 &fac_correct_rh/rho_cfdc
3814 rhi = (ss_i+1.)*100.
3815 if(rhi < 0.) rhi = 0.
3816 n_in_solo_star = 1000.e6*(7.7211e-5 * rhi - 9.2688e-3)/rho_aida
3819 n_in_dust_ultra = 0.;
3821 do ij = 1, ijstop_dust
3824 mu = n_in_ultra*alpha_dust*h_frac_dust*ahet_dust(ij)/
3827 ahet_aux = pie*d_grid_dust(ij)*d_grid_dust(ij)*d_grid_dust(ij)*
3828 &4.73*acorr_dust/6.0
3829 mu = n_in_ultra*alpha_dust*h_frac_dust*ahet_aux/base_dust_omega
3834 naux = (1. - exp(-mu))*n_grid_dust(ij)
3835 n_in_dust_ultra = n_in_dust_ultra + naux
3836 dn_in_dust = max(mu*(n_grid_dust(ij)-naux)*(dnaux +dh_frac_dust),
3845 n_in_soot_ultra = 0.0
3847 do ij = 1, ijstop_soot
3850 mu = n_in_ultra*alpha_soot*h_frac_soot*ahet_bc/
3851 &base_soot_philic_omega
3853 ahet_aux = pie*d_grid_soot(ij)*d_grid_soot(ij)*d_grid_soot(ij)*
3855 mu = n_in_ultra*alpha_soot*h_frac_soot*ahet_aux/
3856 &base_soot_philic_omega
3861 naux = (1. - exp(-mu))*n_grid_soot(ij)
3862 n_in_soot_ultra = n_in_soot_ultra + naux
3863 dn_in_soot = max(mu*(n_grid_soot(ij)-naux)*(dnaux +dh_frac_soot),
3881 n_in_solo = psi_solo*glass_frac*h_frac_solo*n_in_solo_star
3882 dn_in_solo =max(psi_solo*glass_frac* (h_frac_solo*77211.0*100.0/
3883 &rho_aida + n_in_solo_star*dh_frac_solo), 0.0)
3887 n_in_ultra = 0.; n_in_dust_ultra = 0.; n_in_soot_ultra =
3888 & 0.; n_in_solo = 0.; n_in_bio_ultra = 0.;
3893 n_in_dust = n_in_dust + n_in_dust_ultra;
3894 n_in_soot = n_in_soot + n_in_soot_ultra;
3895 n_in_bio = n_in_bio + n_in_bio_ultra;
3903 if(n_in_dust + n_in_bio + n_in_soot + n_in_solo > 0.)
then
3905 cihenc_dust = n_in_dust - nin_a_nuc_dust
3906 if(cihenc_dust < 0.) cihenc_dust = 0.
3908 cihenc_soot = n_in_soot - nin_a_nuc_soot
3909 if(cihenc_soot < 0.) cihenc_soot = 0.
3911 cihenc_bio = n_in_bio - nin_a_nuc_bio
3912 if(cihenc_bio < 0.) cihenc_bio = 0.
3914 cihenc_solo = n_in_solo - nin_a_nuc_solo
3915 if(cihenc_solo < 0.) cihenc_solo = 0.
3918 n_iw = n_iw + cihenc_dust
3919 nin_a_nuc_dust = nin_a_nuc_dust + cihenc_dust
3920 num_ic_dust_imm = num_ic_dust_imm + cihenc_dust
3922 n_iw = n_iw + cihenc_soot
3923 nin_a_nuc_soot = nin_a_nuc_soot + cihenc_soot
3924 num_ic_soot_imm = num_ic_soot_imm + cihenc_soot
3926 n_iw = n_iw + cihenc_bio
3927 nin_a_nuc_bio = nin_a_nuc_bio + cihenc_bio
3928 num_ic_bio_imm = num_ic_bio_imm + cihenc_bio
3930 n_iw = n_iw + cihenc_solo
3931 nin_a_nuc_solo = nin_a_nuc_solo + cihenc_solo
3932 num_ic_solo_imm = num_ic_solo_imm + cihenc_solo
3940 n_iw = n_in_dust + n_in_bio + n_in_soot + n_in_solo
3941 dnall = dn_in_dust + dn_in_bio + dn_in_soot + dn_in_solo
3943 if (n_iw .gt. zero_par)
then
3944 fdust_dep = dble(n_in_dust/n_iw)
3946 nhet_dep = dble(n_in_dust)
3947 nhet_dhf = dble(n_in_solo)
3952 if (( dnall > 0.) .and. (n_iw .gt. 0.0))
then
3953 dsh=max(min(n_iw/dnall, ss_i), 0.005)