diff --git a/src/activation.F90 b/src/activation.F90 index 3c2ec85..4884da8 100644 --- a/src/activation.F90 +++ b/src/activation.F90 @@ -261,14 +261,12 @@ subroutine activate(dt, cloud_mass, cloud_number, w, rho, dnumber, dmac, T, p, activated_arg = sum(dnccn_all(:)) + dactive = zero_real_wp + dmass_d = zero_real_wp + if (Smax_act > smax_act_min .and. .not. l_warm) then - if (iopt_inuc < 4) then - ! For lower-order ice nucleation options, need to initialise - ! dactive and dmass_d - dactive = zero_real_wp - dmass_d = zero_real_wp - else + if (iopt_inuc > 3) then ! For higher-order ice nucleation schemes, dactive and dmass_d ! can be based on dustphys dactive = 0.01*dustphys%N(1) diff --git a/src/aerosol_routines.F90 b/src/aerosol_routines.F90 index ce2d505..c74bad0 100644 --- a/src/aerosol_routines.F90 +++ b/src/aerosol_routines.F90 @@ -14,6 +14,9 @@ module aerosol_routines use lognormal_funcs, only: MNtoRm ! DPG - added this for MNtoRm since was ! causing circular conflicts as wanted ! to use it in which_mode_to_use.F90 + USE umprintmgr, ONLY: umPrint, umMessage + + implicit none @@ -209,8 +212,8 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, & character(len=*), parameter :: RoutineName='EXAMINE_AEROSOL' - real, parameter :: aero_mact_mean_max=2.3e-15 !kg 3xsigma+mean(0.3micron) density 1777 - real, parameter :: dust_mact_mean_max=3.1e-13 !kg 3xsigma+mean(1micron) density 1777 + real, parameter :: aero_mact_mean_max=1e30!3e-13 !kg -make bigger for coarse soluble aerosol + real, parameter :: dust_mact_mean_max=1e30!3e-13 !kg 3xsigma+mean(1micron) density 1777 ! When l_bypass_whichmode is True (or no larger modes exist to move aerosol to) ! aerosol mean sizes in hydrometeors can become so large that they lead to instability ! during hydrometeor sedimentation. Limiting the maximum mean mass stops these edge @@ -315,7 +318,18 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, & aeroact(k)%mact=mac aeroact(k)%rcrit=0.0 aeroact(k)%mact_mean=aeroact(k)%mact /(aeroact(k)%nact + epsilon(ntot)) - aeroact(k)%mact_mean=min(aero_mact_mean_max, aeroact(k)%mact_mean) + if (aeroact(k)%mact_mean .gt. aero_mact_mean_max)then +! WRITE(umMessage, *) 'PRFx1 ',aeroact(k)%mact_mean,mac,ntot,ntot*aero_mact_mean_max/mac +! CALL umPrint( umMessage) + aeroact(k)%mact_mean=min(aero_mact_mean_max, aeroact(k)%mact_mean) +! aeroact(k)%mact=aero_mact_mean_max*aeroact(k)%nact +! mac=aeroact(k)%mact + + aeroact(k)%nact=aeroact(k)%mact/aero_mact_mean_max + ntot=aeroact(k)%nact + + + endif ! Get mean radius of distribution rm_arc=MNtoRm(aeroact(k)%mact,aeroact(k)%nact,density,sigma_arc) @@ -340,7 +354,7 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, & aeroact(k)%rcrit2=rcrit2 aeroact(k)%mact2=mact2 aeroact(k)%mact2_mean=aeroact(k)%mact2/(aeroact(k)%nact2 + epsilon(ntot)) - aeroact(k)%mact2_mean=min(aero_mact_mean_max, aeroact(k)%mact2_mean) +! aeroact(k)%mact2_mean=min(aero_mact_mean_max, aeroact(k)%mact2_mean) end if @@ -348,7 +362,7 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, & aeroact(k)%mact1=max(0.0_wp, aeroact(k)%mact-aeroact(k)%mact2) aeroact(k)%rcrit1=0.0 aeroact(k)%mact1_mean=aeroact(k)%mact1/(aeroact(k)%nact1+epsilon(mar)) - aeroact(k)%mact1_mean=min(aero_mact_mean_max, aeroact(k)%mact1_mean) +! aeroact(k)%mact1_mean=min(aero_mact_mean_max, aeroact(k)%mact1_mean) if (cloud_number > epsilon(1.0_wp)) then aeroact(k)%nratio1=max(0.0,min(1.0,aeroact(k)%nact1/cloud_number)) @@ -457,14 +471,27 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, & dustact(k)%mact=mad dustact(k)%rcrit=0.0 dustact(k)%mact_mean=dustact(k)%mact/(dustact(k)%nact+epsilon(mad)) - dustact(k)%mact_mean=min(dust_mact_mean_max, dustact(k)%mact_mean) +! dustact(k)%mact_mean=min(dust_mact_mean_max, dustact(k)%mact_mean) + if (dustact(k)%mact_mean .gt. dust_mact_mean_max)then +! WRITE(umMessage, *) 'PRFx2 ',dustact(k)%mact_mean,mad,nitot,nitot*dust_mact_mean_max/mad +! CALL umPrint( umMessage) + dustact(k)%mact_mean=min(dust_mact_mean_max, dustact(k)%mact_mean) +! dustact(k)%mact=dust_mact_mean_max*dustact(k)%nact +! mad=dustact(k)%mact + + dustact(k)%nact=dustact(k)%mact/dust_mact_mean_max + nitot=dustact(k)%nact + + + endif if (snow_number > epsilon(1.0_wp)) then dustact(k)%nact2=nitot*ratio_s dustact(k)%rcrit2=0.0 dustact(k)%mact2=mad *ratio_s dustact(k)%mact2_mean=dustact(k)%mact2/(dustact(k)%nact2+epsilon(nhtot)) - dustact(k)%mact2_mean=min(dust_mact_mean_max, dustact(k)%mact2_mean) +! dustact(k)%mact2_mean=min(dust_mact_mean_max, dustact(k)%mact2_mean) + dustact(k)%nratio2=max(0.0,min(1.0,dustact(k)%nact2/(snow_number+epsilon(nhtot)) )) end if @@ -473,7 +500,7 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, & dustact(k)%rcrit3=0.0 dustact(k)%mact3=mad *ratio_g dustact(k)%mact3_mean=dustact(k)%mact3/(dustact(k)%nact3+epsilon(nhtot)) - dustact(k)%mact3_mean=min(dust_mact_mean_max, dustact(k)%mact3_mean) +! dustact(k)%mact3_mean=min(dust_mact_mean_max, dustact(k)%mact3_mean) dustact(k)%nratio3=max(0.0,min(1.0,dustact(k)%nact3/(graupel_number+epsilon(nhtot)) )) end if @@ -482,7 +509,8 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, & dustact(k)%mact1=max(0.0_wp, dustact(k)%mact-dustact(k)%mact2-dustact(k)%mact3) dustact(k)%rcrit1=0.0 dustact(k)%mact1_mean=dustact(k)%mact1/(dustact(k)%nact1+epsilon(mar)) - dustact(k)%mact1_mean=min(dust_mact_mean_max, dustact(k)%mact1_mean) +! dustact(k)%mact1_mean=min(dust_mact_mean_max, dustact(k)%mact1_mean) + dustact(k)%nratio1=max(0.0,min(1.0,dustact(k)%nact1/(ice_number+epsilon(nhtot)) )) end if @@ -559,14 +587,25 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, & aeroice(k)%mact=maai aeroice(k)%rcrit=0.0 aeroice(k)%mact_mean=aeroice(k)%mact/(aeroice(k)%nact+epsilon(nhtot)) - aeroice(k)%mact_mean=min(aero_mact_mean_max, aeroice(k)%mact_mean) +! aeroice(k)%mact_mean=min(aero_mact_mean_max, aeroice(k)%mact_mean) + if (aeroice(k)%mact_mean .gt. aero_mact_mean_max)then +! WRITE(umMessage, *) 'PRFx3 ',aeroice(k)%mact_mean,maai,nitot,nitot*aero_mact_mean_max/maai +! CALL umPrint( umMessage) + aeroice(k)%mact_mean=min(aero_mact_mean_max, aeroice(k)%mact_mean) +! aeroice(k)%mact=aero_mact_mean_max*aeroice(k)%nact +! maai=aeroice(k)%mact + aeroice(k)%nact=aeroice(k)%mact/aero_mact_mean_max + nitot=aeroice(k)%nact + + endif + if (ratio_s > epsilon(1.0_wp)) then aeroice(k)%nact2=nitot*ratio_s aeroice(k)%rcrit2=0.0 aeroice(k)%mact2=maai*ratio_s aeroice(k)%mact2_mean=aeroice(k)%mact2/(aeroice(k)%nact2+epsilon(nhtot)) - aeroice(k)%mact2_mean=min(aero_mact_mean_max, aeroice(k)%mact2_mean) +! aeroice(k)%mact2_mean=min(aero_mact_mean_max, aeroice(k)%mact2_mean) aeroice(k)%nratio2=max(0.0, min(1.0,aeroice(k)%nact2/(snow_number+epsilon(nhtot)) )) end if @@ -647,6 +686,16 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, & dustliq(k)%rcrit=0.0 dustliq(k)%mact_mean=dustliq(k)%mact/(dustliq(k)%nact+epsilon(nhtot)) dustliq(k)%mact_mean=min(dust_mact_mean_max, dustliq(k)%mact_mean) + if (dustliq(k)%mact_mean .gt. dust_mact_mean_max)then +! WRITE(umMessage, *) 'PRFx4 ',dustliq(k)%mact_mean,madl,ntot,ntot*aero_mact_mean_max/madl +! CALL umPrint( umMessage) + dustliq(k)%mact_mean=min(dust_mact_mean_max, dustliq(k)%mact_mean) +! dustliq(k)%mact=dust_mact_mean_max*dustliq(k)%nact +! madl=dustliq(k)%mact + dustliq(k)%nact=dustliq(k)%mact/dust_mact_mean_max + ntot=dustliq(k)%nact + + endif if (ratio_r > epsilon(1.0_wp)) then dustliq(k)%nact2=ntot*ratio_r diff --git a/src/mphys_tidy.F90 b/src/mphys_tidy.F90 index fa0b2bb..04209cc 100644 --- a/src/mphys_tidy.F90 +++ b/src/mphys_tidy.F90 @@ -2,7 +2,8 @@ module mphys_tidy use variable_precision, only: wp use process_routines, only: process_rate, process_name use aerosol_routines, only: aerosol_active - use thresholds, only: thresh_tidy, thresh_atidy + use thresholds, only: thresh_tidy, thresh_atidy,& + aeromass_small, aeronumber_small use passive_fields, only: exner, pressure use mphys_switches, only: & i_qv, i_ql, i_nl, i_qr, i_nr, i_m3r, i_th, & @@ -14,6 +15,8 @@ module mphys_tidy i_an2, i_am2, i_am4, i_am5, l_warm, & i_an6, i_am6, i_am7, i_am8, i_am9, & i_an11, i_an12, & + i_an1,i_am1, i_an2, i_am2, i_an3, i_am3, & + i_am6, i_an6, i_am10, i_an10, & l_process, ntotalq, ntotala, & i_qstart, i_nstart, i_m3start, & l_separate_rain, l_tidy_conserve_E, l_tidy_conserve_q, & @@ -208,12 +211,16 @@ subroutine qtidy(ixy_inner, dt, nz, qfields, procs, aerofields, aeroact, dustact !-------------------------------------------------------------------------- IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) - l_qsig(0)=.false. - l_qsneg(0)=.false. + l_qsig(:)=.false. + l_qsneg(:)=.false. !should this be l_qsneg(:) do k = 1, nz nr_reset=.false. + qr_reset=.false. m3r_reset=.false. + nl_reset=.false. + ql_reset=.false. + ni_reset=.false. qi_reset=.false. qs_reset=.false. ns_reset=.false. @@ -865,12 +872,88 @@ subroutine tidy_ain(qfields, aerofields) end if if (i_an11 > 0) then if (((qfields(k, i_ql) + qfields(k,i_qr) + qfields(k,i_qi) + qfields(k,i_qs) + qfields(k,i_qg) & - <=0.0 .and. aerofields(k,i_an12)>0.0) & + <=0.0 .and. aerofields(k,i_an11)>0.0) & .or. aerofields(k,i_an11) < 0.0)) then aerofields(k,i_an11)=0.0 end if end if + + +! aerofields(:, i_am1)=a1(ks:ke,i,j) AitkenSolMass +! aerofields(:, i_an1)=a2(ks:ke,i,j) AitkenSolNumber +! aerofields(:, i_am2)=a3(ks:ke,i,j) AccumSolMass +! aerofields(:, i_an2)=a4(ks:ke,i,j) AccumSolNumber +! aerofields(:, i_am3)=a5(ks:ke,i,j) CoarseSolMass +! aerofields(:, i_an3)=a6(ks:ke,i,j) CoarseSolNumber +! aerofields(:, i_am4)=a7(ks:ke,i,j) ActSolLiq_casim +! aerofields(:, i_am5)=a8(ks:ke,i,j) ActSolRain_casim <-this is not currently available +! aerofields(:, i_am6)=a9(ks:ke,i,j) CoarseDustMass +! aerofields(:, i_an6)=a10(ks:ke,i,j) CoarseDustNumber + + IF (i_am1 > 0) then + if (aerofields(k, i_am1) < aeromass_small) then + aerofields(k, i_am1)=0.0 + aerofields(k, i_an1)=0.0 + end if + END IF + IF (i_am2 > 0) then + if (aerofields(k, i_am2) < aeromass_small) then + aerofields(k, i_am2)=0.0 + aerofields(k, i_an2)=0.0 + end if + END IF + IF (i_am3 > 0) then + if (aerofields(k, i_am3) < aeromass_small) then + aerofields(k, i_am3)=0.0 + aerofields(k, i_an3)=0.0 + end if + END IF + IF (i_am6 > 0) then + if (aerofields(k, i_am6) < aeromass_small) then + aerofields(k, i_am6)=0.0 + aerofields(k, i_an6)=0.0 + end if + END IF + IF (i_an10 > 0) then + if (aerofields(k, i_am10) < aeromass_small) then + aerofields(k, i_am10)=0.0 + aerofields(k, i_an10)=0.0 + end if + END IF + + IF (i_am1 > 0) then + if (aerofields(k, i_an1) < aeronumber_small) then + aerofields(k, i_am1)=0.0 + aerofields(k, i_an1)=0.0 + end if + END IF + IF (i_am2 > 0) then + if (aerofields(k, i_an2) < aeronumber_small) then + aerofields(k, i_am2)=0.0 + aerofields(k, i_an2)=0.0 + end if + END IF + IF (i_am3 > 0) then + if (aerofields(k, i_an3) < aeronumber_small) then + aerofields(k, i_am3)=0.0 + aerofields(k, i_an3)=0.0 + end if + END IF + IF (i_am6 > 0) then + if (aerofields(k, i_an6) < aeronumber_small) then + aerofields(k, i_am6)=0.0 + aerofields(k, i_an6)=0.0 + end if + END IF + IF (i_am10 > 0) then + if (aerofields(k, i_an10) < aeronumber_small) then + aerofields(k, i_am10)=0.0 + aerofields(k, i_an10)=0.0 + end if + END IF + + end do IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) diff --git a/src/sedimentation.F90 b/src/sedimentation.F90 index 700d372..70fb180 100644 --- a/src/sedimentation.F90 +++ b/src/sedimentation.F90 @@ -945,30 +945,38 @@ subroutine sedr_1M_2M(ixy_inner, step_length, qfields, aeroact, dustact, & !============================ if (l_ased .and. l_da_local) then if (params%id == cloud_params%id) then - dmac=(flux_n2(k+1)*aeroact(k+1)%nratio1*aeroact(k+1)%mact1_mean - & - flux_n2(k)*aeroact(k)%nratio1*aeroact(k)%mact1_mean)* rdz_on_rho(k,ixy_inner) + dmac=(min(flux_n2(k+1)*aeroact(k+1)%nratio1*aeroact(k+1)%mact1_mean, & + aeroact(k+1)%mact1/(step_length*rdz_on_rho(k+1,ixy_inner))) & + -min(flux_n2(k)*aeroact(k)%nratio1*aeroact(k)%mact1_mean, & + aeroact(k)%mact1/(step_length*rdz_on_rho(k,ixy_inner))))* rdz_on_rho(k,ixy_inner) if (l_passivenumbers) then dnumber_a=(flux_n2(k+1)*aeroact(k+1)%nratio1 - & flux_n2(k)*aeroact(k)%nratio1)* rdz_on_rho(k,ixy_inner) end if if (.not. l_warm) then - dmad=(flux_n2(k+1)*dustact(k+1)%nratio1*dustact(k+1)%mact1_mean- & - flux_n2(k)*dustact(k)%nratio1*dustact(k)%mact1_mean)*rdz_on_rho(k,ixy_inner) + dmad=(min(flux_n2(k+1)*dustact(k+1)%nratio1*dustact(k+1)%mact1_mean, & + dustact(k+1)%mact1/(step_length*rdz_on_rho(k+1,ixy_inner)))& + -min(flux_n2(k)*dustact(k)%nratio1*dustact(k)%mact1_mean, & + dustact(k)%mact1/(step_length*rdz_on_rho(k,ixy_inner))))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers_ice .and. dustact(k)%mact_mean > 0.0) then dnumber_d=(flux_n2(k+1)*dustact(k+1)%nratio1- & flux_n2(k)*dustact(k)%nratio1)*rdz_on_rho(k,ixy_inner) end if end if else if (params%id == rain_params%id) then - dmac=(flux_n2(k+1)*aeroact(k+1)%nratio2*aeroact(k+1)%mact2_mean- & - flux_n2(k)*aeroact(k)%nratio2*aeroact(k)%mact2_mean)*rdz_on_rho(k,ixy_inner) + dmac=(min(flux_n2(k+1)*aeroact(k+1)%nratio2*aeroact(k+1)%mact2_mean, & + aeroact(k+1)%mact2/(step_length*rdz_on_rho(k+1,ixy_inner))) & + -min(flux_n2(k)*aeroact(k)%nratio2*aeroact(k)%mact2_mean, & + aeroact(k)%mact2/(step_length*rdz_on_rho(k,ixy_inner))))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers .and. aeroact(k)%mact_mean > 0.0) then dnumber_a=(flux_n2(k+1)*aeroact(k+1)%nratio2- & flux_n2(k)*aeroact(k)%nratio2)*rdz_on_rho(k,ixy_inner) end if if (.not. l_warm) then - dmad=(flux_n2(k+1)*dustact(k+1)%nratio2*dustact(k+1)%mact2_mean- & - flux_n2(k)*dustact(k)%nratio2*dustact(k)%mact2_mean)*rdz_on_rho(k,ixy_inner) + dmad=(min(flux_n2(k+1)*dustact(k+1)%nratio2*dustact(k+1)%mact2_mean, & + dustact(k+1)%mact2/(step_length*rdz_on_rho(k+1,ixy_inner)))& + -min(flux_n2(k)*dustact(k)%nratio2*dustact(k)%mact2_mean, & + dustact(k)%mact2/(step_length*rdz_on_rho(k,ixy_inner))))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers_ice) then dnumber_d=(flux_n2(k+1)*dustact(k+1)%nratio2- & flux_n2(k)*dustact(k)%nratio2)*rdz_on_rho(k,ixy_inner) @@ -977,10 +985,14 @@ subroutine sedr_1M_2M(ixy_inner, step_length, qfields, aeroact, dustact, & end if if (params%id == ice_params%id) then - dmac = (flux_n2(k+1)*aeroact(k+1)%nratio1*aeroact(k+1)%mact1_mean- & - flux_n2(k)*aeroact(k)%nratio1*aeroact(k)%mact1_mean)*rdz_on_rho(k,ixy_inner) - dmad = (flux_n2(k+1)*dustact(k+1)%nratio1*dustact(k+1)%mact1_mean- & - flux_n2(k)*dustact(k)%nratio1*dustact(k)%mact1_mean)*rdz_on_rho(k,ixy_inner) + dmac = (min(flux_n2(k+1)*aeroact(k+1)%nratio1*aeroact(k+1)%mact1_mean, & + aeroact(k+1)%mact1/(step_length*rdz_on_rho(k+1,ixy_inner)))& + -min(flux_n2(k)*aeroact(k)%nratio1*aeroact(k)%mact1_mean, & + aeroact(k)%mact1/(step_length*rdz_on_rho(k,ixy_inner))))*rdz_on_rho(k,ixy_inner) + dmad = (min(flux_n2(k+1)*dustact(k+1)%nratio1*dustact(k+1)%mact1_mean, & + dustact(k+1)%mact1/(step_length*rdz_on_rho(k+1,ixy_inner)))& + -min(flux_n2(k)*dustact(k)%nratio1*dustact(k)%mact1_mean, & + dustact(k)%mact1/(step_length*rdz_on_rho(k,ixy_inner))))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers) then dnumber_a=(flux_n2(k+1)*aeroact(k+1)%nratio1- & flux_n2(k)*aeroact(k)%nratio1)*rdz_on_rho(k,ixy_inner) @@ -990,10 +1002,14 @@ subroutine sedr_1M_2M(ixy_inner, step_length, qfields, aeroact, dustact, & flux_n2(k)*dustact(k)%nratio1)*rdz_on_rho(k,ixy_inner) end if else if (params%id == snow_params%id) then - dmac=(flux_n2(k+1)*aeroact(k+1)%nratio2*aeroact(k+1)%mact2_mean- & - flux_n2(k)*aeroact(k)%nratio2*aeroact(k)%mact2_mean)*rdz_on_rho(k,ixy_inner) - dmad=(flux_n2(k+1)*dustact(k+1)%nratio2*dustact(k+1)%mact2_mean- & - flux_n2(k)*dustact(k)%nratio2*dustact(k)%mact2_mean)*rdz_on_rho(k,ixy_inner) + dmac=(min(flux_n2(k+1)*aeroact(k+1)%nratio2*aeroact(k+1)%mact2_mean, & + aeroact(k+1)%mact2/(step_length*rdz_on_rho(k+1,ixy_inner))) & + -min(flux_n2(k)*aeroact(k)%nratio2*aeroact(k)%mact2_mean, & + aeroact(k)%mact2/(step_length*rdz_on_rho(k,ixy_inner))))*rdz_on_rho(k,ixy_inner) + dmad=(min(flux_n2(k+1)*dustact(k+1)%nratio2*dustact(k+1)%mact2_mean, & + dustact(k+1)%mact2/(step_length*rdz_on_rho(k+1,ixy_inner))) & + -min(flux_n2(k)*dustact(k)%nratio2*dustact(k)%mact2_mean, & + dustact(k)%mact2/(step_length*rdz_on_rho(k,ixy_inner))))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers) then dnumber_a=(flux_n2(k+1)*aeroact(k+1)%nratio2- & flux_n2(k)*aeroact(k)%nratio2)*rdz_on_rho(k,ixy_inner) @@ -1003,10 +1019,14 @@ subroutine sedr_1M_2M(ixy_inner, step_length, qfields, aeroact, dustact, & flux_n2(k)*dustact(k)%nratio2)*rdz_on_rho(k,ixy_inner) end if else if (params%id == graupel_params%id) then - dmac=(flux_n2(k+1)*aeroact(k+1)%nratio3*aeroact(k+1)%mact3_mean- & - flux_n2(k)*aeroact(k)%nratio3*aeroact(k)%mact3_mean)*rdz_on_rho(k,ixy_inner) - dmad=(flux_n2(k+1)*dustact(k+1)%nratio3*dustact(k+1)%mact3_mean- & - flux_n2(k)*dustact(k)%nratio3*dustact(k)%mact3_mean)*rdz_on_rho(k,ixy_inner) + dmac=(min(flux_n2(k+1)*aeroact(k+1)%nratio3*aeroact(k+1)%mact3_mean, & + aeroact(k+1)%mact3/(step_length*rdz_on_rho(k+1,ixy_inner))) & + -min(flux_n2(k)*aeroact(k)%nratio3*aeroact(k)%mact3_mean, & + aeroact(k)%mact3/(step_length*rdz_on_rho(k,ixy_inner))))*rdz_on_rho(k,ixy_inner) + dmad=(min(flux_n2(k+1)*dustact(k+1)%nratio3*dustact(k+1)%mact3_mean, & + dustact(k+1)%mact3/(step_length*rdz_on_rho(k+1,ixy_inner))) & + -min(flux_n2(k)*dustact(k)%nratio3*dustact(k)%mact3_mean, & + dustact(k)%mact3/(step_length*rdz_on_rho(k,ixy_inner))))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers) then dnumber_a=(flux_n2(k+1)*aeroact(k+1)%nratio3- & @@ -1027,23 +1047,27 @@ subroutine sedr_1M_2M(ixy_inner, step_length, qfields, aeroact, dustact, & !============================ if (l_ased .and. l_da_local) then if (params%id == cloud_params%id) then - dmac=(flux_n2(k+1)*aeroact(k+1)%nratio1*aeroact(k+1)%mact1_mean)*rdz_on_rho(k,ixy_inner) + dmac=min(flux_n2(k+1)*aeroact(k+1)%nratio1*aeroact(k+1)%mact1_mean, & + aeroact(k+1)%mact1/(step_length*rdz_on_rho(k+1,ixy_inner)))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers) then dnumber_a=(flux_n2(k+1)*aeroact(k+1)%nratio1)*rdz_on_rho(k,ixy_inner) end if if (.not. l_warm) then - dmad=(flux_n2(k+1)*dustact(k+1)%nratio1*dustact(k+1)%mact1_mean)*rdz_on_rho(k,ixy_inner) + dmad=min(flux_n2(k+1)*dustact(k+1)%nratio1*dustact(k+1)%mact1_mean, & + dustact(k+1)%mact1/(step_length*rdz_on_rho(k+1,ixy_inner)))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers_ice) then dnumber_d=(flux_n2(k+1)*dustact(k+1)%nratio1)*rdz_on_rho(k,ixy_inner) end if end if else if (params%id == rain_params%id) then - dmac=(flux_n2(k+1)*aeroact(k+1)%nratio2*aeroact(k+1)%mact2_mean)*rdz_on_rho(k,ixy_inner) + dmac=min(flux_n2(k+1)*aeroact(k+1)%nratio2*aeroact(k+1)%mact2_mean, & + aeroact(k+1)%mact2/(step_length*rdz_on_rho(k+1,ixy_inner)))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers) then dnumber_a=(flux_n2(k+1)*aeroact(k+1)%nratio2)*rdz_on_rho(k,ixy_inner) end if if (.not. l_warm) then - dmad=(flux_n2(k+1)*dustact(k+1)%nratio2*dustact(k+1)%mact2_mean)*rdz_on_rho(k,ixy_inner) + dmad=min(flux_n2(k+1)*dustact(k+1)%nratio2*dustact(k+1)%mact2_mean, & + dustact(k+1)%mact2/(step_length*rdz_on_rho(k+1,ixy_inner)))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers_ice) then dnumber_d=flux_n2(k+1)*dustact(k+1)%nratio2*rdz_on_rho(k,ixy_inner) end if @@ -1051,8 +1075,10 @@ subroutine sedr_1M_2M(ixy_inner, step_length, qfields, aeroact, dustact, & end if if (params%id == ice_params%id) then - dmac=(flux_n2(k+1)*aeroact(k+1)%nratio1*aeroact(k+1)%mact1_mean)*rdz_on_rho(k,ixy_inner) - dmad=(flux_n2(k+1)*dustact(k+1)%nratio1*dustact(k+1)%mact1_mean)*rdz_on_rho(k,ixy_inner) + dmac=min(flux_n2(k+1)*aeroact(k+1)%nratio1*aeroact(k+1)%mact1_mean, & + aeroact(k+1)%mact1/(step_length*rdz_on_rho(k+1,ixy_inner)))*rdz_on_rho(k,ixy_inner) + dmad=min(flux_n2(k+1)*dustact(k+1)%nratio1*dustact(k+1)%mact1_mean, & + dustact(k+1)%mact1/(step_length*rdz_on_rho(k+1,ixy_inner)))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers) then dnumber_a=(flux_n2(k+1)*aeroact(k+1)%nratio1)*rdz_on_rho(k,ixy_inner) end if @@ -1060,8 +1086,10 @@ subroutine sedr_1M_2M(ixy_inner, step_length, qfields, aeroact, dustact, & dnumber_d=(flux_n2(k+1)*dustact(k+1)%nratio1)*rdz_on_rho(k,ixy_inner) end if else if (params%id == snow_params%id) then - dmac=(flux_n2(k+1)*aeroact(k+1)%nratio2*aeroact(k+1)%mact2_mean)*rdz_on_rho(k,ixy_inner) - dmad=(flux_n2(k+1)*dustact(k+1)%nratio2*dustact(k+1)%mact2_mean)*rdz_on_rho(k,ixy_inner) + dmac=min(flux_n2(k+1)*aeroact(k+1)%nratio2*aeroact(k+1)%mact2_mean, & + aeroact(k+1)%mact2/(step_length*rdz_on_rho(k+1,ixy_inner)))*rdz_on_rho(k,ixy_inner) + dmad=min(flux_n2(k+1)*dustact(k+1)%nratio2*dustact(k+1)%mact2_mean, & + dustact(k+1)%mact2/(step_length*rdz_on_rho(k+1,ixy_inner)))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers) then dnumber_a=(flux_n2(k+1)*aeroact(k+1)%nratio2)*rdz_on_rho(k,ixy_inner) end if @@ -1070,8 +1098,10 @@ subroutine sedr_1M_2M(ixy_inner, step_length, qfields, aeroact, dustact, & end if else if (params%id == graupel_params%id) then if (i_aerosed_method==1) then - dmac=(flux_n2(k+1)*aeroact(k+1)%nratio3*aeroact(k+1)%mact3_mean)*rdz_on_rho(k,ixy_inner) - dmad=(flux_n2(k+1)*dustact(k+1)%nratio3*dustact(k+1)%mact3_mean)*rdz_on_rho(k,ixy_inner) + dmac=min(flux_n2(k+1)*aeroact(k+1)%nratio3*aeroact(k+1)%mact3_mean, & + aeroact(k+1)%mact3/(step_length*rdz_on_rho(k+1,ixy_inner)))*rdz_on_rho(k,ixy_inner) + dmad=min(flux_n2(k+1)*dustact(k+1)%nratio3*dustact(k+1)%mact3_mean, & + dustact(k+1)%mact3/(step_length*rdz_on_rho(k+1,ixy_inner)))*rdz_on_rho(k,ixy_inner) if (l_passivenumbers) then dnumber_a=(flux_n2(k+1)*aeroact(k+1)%nratio3)*rdz_on_rho(k,ixy_inner) end if diff --git a/src/special.F90 b/src/special.F90 index 9286209..003cb56 100644 --- a/src/special.F90 +++ b/src/special.F90 @@ -240,7 +240,7 @@ function casim_erfc(x) !-------------------------------------------------------------------------- IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) - casim_erfc=erfg(x,c) + casim_erfc=erfc(x) !erfg(x,c) ! bug in erfg for c=1 and x<1e-10 IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle)