diff --git a/src/atmosphere/photochem_atmosphere_rhs.f90 b/src/atmosphere/photochem_atmosphere_rhs.f90 index 76abe24..0291dac 100644 --- a/src/atmosphere/photochem_atmosphere_rhs.f90 +++ b/src/atmosphere/photochem_atmosphere_rhs.f90 @@ -14,7 +14,7 @@ #:for TYPE1, NAME in TYPES_NAMES subroutine dochem_${NAME}$(self, neqs, nsp, np, nsl, nq, nz, trop_ind, nrT, usol, density, rx_rates, & gas_sat_den, molecules_per_particle, & - H2O_sat_mix, H2O_rh, rainout_rates, bg_density, & + H2O_sat_mix, H2O_rh, rainout_rates, scale_height, wfall, bg_density, & densities, xp, xl, rhs) use photochem_enum, only: CondensingParticle use photochem_common, only: chempl, chempl_sl @@ -31,13 +31,13 @@ subroutine dochem_${NAME}$(self, neqs, nsp, np, nsl, nq, nz, trop_ind, nrT, usol real(dp), intent(in) :: gas_sat_den(np,nz) real(dp), intent(in) :: molecules_per_particle(np,nz) real(dp), intent(in) :: H2O_sat_mix(nz), H2O_rh(trop_ind) - real(dp), intent(in) :: rainout_rates(nq, trop_ind) + real(dp), intent(in) :: rainout_rates(nq, trop_ind), scale_height(nz), wfall(np,nz) real(dp), intent(in) :: bg_density(nz) ${TYPE1}$, intent(inout) :: densities(nsp+1,nz), xp(nz), xl(nz) ${TYPE1}$, intent(inout) :: rhs(neqs) - real(dp) :: H2O_cold_trap - ${TYPE1}$ :: dn_gas_dt, dn_particle_dt, cond_rate + real(dp) :: H2O_cold_trap, cond_rate0 + ${TYPE1}$ :: rh, df_gas_dt, df_particle_dt, cond_rate integer :: i, ii, j, k, kk type(PhotochemData), pointer :: dat type(PhotochemVars), pointer :: var @@ -46,8 +46,9 @@ subroutine dochem_${NAME}$(self, neqs, nsp, np, nsl, nq, nz, trop_ind, nrT, usol var => self%var #:if NAME == 'dual' - dn_gas_dt = dual(size(usol(1,1)%der)) - dn_particle_dt = dual(size(usol(1,1)%der)) + rh = dual(size(usol(1,1)%der)) + df_gas_dt = dual(size(usol(1,1)%der)) + df_particle_dt = dual(size(usol(1,1)%der)) cond_rate = dual(size(usol(1,1)%der)) #:endif do j = 1,var%nz @@ -144,30 +145,49 @@ subroutine dochem_${NAME}$(self, neqs, nsp, np, nsl, nq, nz, trop_ind, nrT, usol ii = dat%particle_gas_phase_ind(i) ! index of gas phase kk = ii + (j - 1) * dat%nq ! gas phase rhs index k = i + (j - 1) * dat%nq ! particle rhs index + + ! compute the relative humidity + rh = max(densities(ii,j)/gas_sat_den(i,j),small_real) ! if the gas phase is super-saturated - if (densities(ii,j) >= gas_sat_den(i,j)*var%condensation_rate(2,i)) then - ! compute condensation rate - cond_rate = damp_condensation_rate(var%condensation_rate(1,i), & - var%condensation_rate(2,i), & - var%condensation_rate(3,i), & - densities(ii,j)/gas_sat_den(i,j)) + if (rh > var%cond_params(i)%RHc) then + ! Condensation occurs + + ! Condensation rate is based on how rapidly gases are mixing + ! in the atmosphere. We also smooth the rate to prevent stiffness. + cond_rate0 = var%cond_params(i)%k_cond*(var%edd(j)/scale_height(j)**2.0_dp) + cond_rate = damp_condensation_rate(cond_rate0, & + var%cond_params(i)%RHc, & + (1.0_dp + var%cond_params(i)%smooth_factor)*var%cond_params(i)%RHc, & + rh) - ! rate the gas molecules are going to particles - ! in molecules/cm3/s - dn_gas_dt = - cond_rate* & - (densities(ii,j) - gas_sat_den(i,j)*var%condensation_rate(2,i)) - ! add to rhs vector, convert to change in mixing ratio/s - rhs(kk) = rhs(kk) + dn_gas_dt/density(j) + ! Rate gases are being destroyed (1/s) + df_gas_dt = - cond_rate*usol(ii,j) + rhs(kk) = rhs(kk) + df_gas_dt - ! rate of particle production from gas condensing - ! in particles/cm3/s - dn_particle_dt = - dn_gas_dt/density(j) - ! add to rhs vector, convert to moles/cm3/s - rhs(k) = rhs(k) + dn_particle_dt - else - ! particles don't change! - rhs(k) = rhs(k) + 0.0_dp + ! Rate particles are being produced (1/s) + df_particle_dt = - df_gas_dt + rhs(k) = rhs(k) + df_particle_dt + + elseif (rh <= var%cond_params(i)%RHc .and. var%evaporation) then + ! Evaporation occurs + + ! Evaporation rate is based on how rapidly particles are falling + ! in the atmosphere. We also smooth the rate to prevent stiffness. + cond_rate0 = var%cond_params(i)%k_evap*(wfall(i,j)/scale_height(j)) + cond_rate = damp_condensation_rate(cond_rate0, & + 1.0_dp/var%cond_params(i)%RHc, & + (1.0_dp + var%cond_params(i)%smooth_factor)/var%cond_params(i)%RHc, & + 1.0_dp/rh) + + ! Rate gases are being produced (1/s) + df_gas_dt = cond_rate*usol(i,j) + rhs(kk) = rhs(kk) + df_gas_dt + + ! Rate particles are being destroyed (1/s) + df_particle_dt = - df_gas_dt + rhs(k) = rhs(k) + df_particle_dt + endif endif enddo @@ -219,14 +239,14 @@ subroutine diffusion_coefficients(dat, var, den, mubar, & binary_diffusion_param => default_binary_diffusion_param endif - ! eddy diffusion. particles and gases + ! eddy diffusion. gases only ! middle do i = 2,var%nz-1 eddav_p = sqrt(var%edd(i)*var%edd(i+1)) eddav_m = sqrt(var%edd(i)*var%edd(i-1)) denav_p = sqrt(den(i)*den(i+1)) denav_m = sqrt(den(i)*den(i-1)) - do j = 1,dat%nq + do j = dat%ng_1,dat%nq DU(j,i) = (eddav_p*denav_p)/(den(i)*var%dz(i)**2.0_dp) DL(j,i) = (eddav_m*denav_m)/(den(i)*var%dz(i)**2.0_dp) DD(j,i) = - DU(j,i) - DL(j,i) @@ -237,12 +257,24 @@ subroutine diffusion_coefficients(dat, var, den, mubar, & eddav_m = sqrt(var%edd(var%nz)*var%edd(var%nz-1)) denav_p = sqrt(den(1)*den(2)) denav_m = sqrt(den(var%nz)*den(var%nz-1)) - do j = 1,dat%nq + do j = dat%ng_1,dat%nq DU(j,1) = (eddav_p*denav_p)/(den(1)*var%dz(1)**2.0_dp) DD(j,1) = - DU(j,1) DL(j,var%nz) = (eddav_m*denav_m)/(den(var%nz)*var%dz(var%nz)**2.0_dp) DD(j,var%nz) = - DL(j,var%nz) enddo + + ! We do not include eddy diffusion for particles + do j = 1,var%nz + do i = 1,dat%np + DU(i,j) = 0.0_dp + DL(i,j) = 0.0_dp + DD(i,j) = 0.0_dp + ADU(i,j) = 0.0_dp + ADL(i,j) = 0.0_dp + ADD(i,j) = 0.0_dp + enddo + enddo ! Molecular diffusion. Only gas species ! middle @@ -503,7 +535,7 @@ module subroutine prep_all_background_gas(self, usol_in, err) wrk%DU, wrk%DD, wrk%DL, wrk%ADU, wrk%ADL, wrk%ADD, & wrk%wfall, wrk%VH2_esc, wrk%VH_esc) - wrk%surface_scale_height = (k_boltz*var%temperature(1)*N_avo)/(wrk%mubar(1)*var%grav(1)) + wrk%scale_height = (k_boltz*var%temperature(:)*N_avo)/(wrk%mubar(:)*var%grav(:)) !!! H and H2 escape wrk%upper_veff_copy = var%upper_veff @@ -608,7 +640,7 @@ module subroutine rhs_background_gas(self, neqs, tn, usol_flat, rhs, err) call dochem(self, var%neqs, dat%nsp, dat%np, dat%nsl, dat%nq, var%nz, & var%trop_ind, dat%nrT, wrk%usol, wrk%density, wrk%rx_rates, & wrk%gas_sat_den, wrk%molecules_per_particle, & - wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%densities(dat%nsp,:), & + wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%scale_height, wrk%wfall, wrk%densities(dat%nsp,:), & wrk%densities, wrk%xp, wrk%xl, rhs) ! Extra functions specifying production or destruction @@ -652,7 +684,7 @@ module subroutine rhs_background_gas(self, neqs, tn, usol_flat, rhs, err) elseif (var%lowerboundcond(i) == MosesBC) then rhs(i) = rhs(i) + wrk%DU(i,1)*wrk%usol(i,2) + wrk%ADU(i,1)*wrk%usol(i,2) & + wrk%DD(i,1)*wrk%usol(i,1) + wrk%ADD(i,1)*wrk%usol(i,1) & - - (var%edd(1)/wrk%surface_scale_height)*wrk%usol(i,1)/var%dz(1) + - (var%edd(1)/wrk%scale_height(1))*wrk%usol(i,1)/var%dz(1) endif enddo @@ -740,7 +772,7 @@ subroutine fcn(x_, f_) call dochem(self, var%neqs, dat%nsp, dat%np, dat%nsl, dat%nq, var%nz, var%trop_ind, dat%nrT, & usol_, wrk%density, wrk%rx_rates, & wrk%gas_sat_den, wrk%molecules_per_particle, & - wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%densities(dat%nsp,:), & + wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%scale_height, wrk%wfall, wrk%densities(dat%nsp,:), & densities, xp, xl, f_) end subroutine end subroutine @@ -802,7 +834,7 @@ module subroutine jac_background_gas(self, lda_neqs, neqs, tn, usol_flat, jac, e call dochem(self, var%neqs, dat%nsp, dat%np, dat%nsl, dat%nq, var%nz, var%trop_ind, dat%nrT, & wrk%usol, wrk%density, wrk%rx_rates, & wrk%gas_sat_den, wrk%molecules_per_particle, & - wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%densities(dat%nsp,:), & + wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%scale_height, wrk%wfall, wrk%densities(dat%nsp,:), & densities, xp, xl, rhs) !$omp parallel private(i, j, k, m, mm, usol_perturb, R, densities, xl, xp, rhs_perturb) usol_perturb = wrk%usol @@ -816,7 +848,7 @@ module subroutine jac_background_gas(self, lda_neqs, neqs, tn, usol_flat, jac, e call dochem(self, var%neqs, dat%nsp, dat%np, dat%nsl, dat%nq, var%nz, var%trop_ind, dat%nrT, & usol_perturb, wrk%density, wrk%rx_rates, & wrk%gas_sat_den, wrk%molecules_per_particle, & - wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%densities(dat%nsp,:), & + wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%scale_height, wrk%wfall, wrk%densities(dat%nsp,:), & densities, xp, xl, rhs_perturb) do m = 1,dat%nq @@ -884,7 +916,7 @@ module subroutine jac_background_gas(self, lda_neqs, neqs, tn, usol_flat, jac, e elseif (var%lowerboundcond(i) == MosesBC) then djac(dat%ku,i+dat%nq) = wrk%DU(i,1) + wrk%ADU(i,1) djac(dat%kd,i) = djac(dat%kd,i) + wrk%DD(i,1) + wrk%ADD(i,1) - & - (var%edd(1)/wrk%surface_scale_height)/var%dz(1) + (var%edd(1)/wrk%scale_height(1))/var%dz(1) ! elseif (var%lowerboundcond(i) == MixDependentFluxBC) then ! djac(dat%ku,i+dat%nq) = wrk%DU(i,1) + wrk%ADU(i,1) ! djac(dat%kd,i) = djac(dat%kd,i) + wrk%DD(i,1) + wrk%ADD(i,1) & @@ -943,7 +975,7 @@ module subroutine right_hand_side_chem(self, usol, rhs, err) call dochem(self, var%neqs, dat%nsp, dat%np, dat%nsl, dat%nq, var%nz, & var%trop_ind, dat%nrT, wrk%usol, wrk%density, wrk%rx_rates, & wrk%gas_sat_den, wrk%molecules_per_particle, & - wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%densities(dat%nsp,:), & + wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%scale_height, wrk%wfall, wrk%densities(dat%nsp,:), & wrk%densities, wrk%xp, wrk%xl, rhs) end subroutine diff --git a/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 b/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 index bcc9064..08bbce8 100644 --- a/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 +++ b/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 @@ -27,7 +27,7 @@ module subroutine equilibrium_climate(self, usol_den, molecules_per_particle, T_ #:for TYPE1, NAME in TYPES_NAMES subroutine dochem_${NAME}$(self, usol, rx_rates, & gas_sat_den, molecules_per_particle, & - H2O_sat_mix, H2O_rh, rainout_rates, & + H2O_sat_mix, H2O_rh, rainout_rates, scale_height, wfall, & density, mix, densities, xp, xl, rhs) use photochem_enum, only: CondensingParticle use photochem_common, only: chempl, chempl_sl @@ -42,14 +42,14 @@ subroutine dochem_${NAME}$(self, usol, rx_rates, & real(dp), intent(in) :: gas_sat_den(:,:) real(dp), intent(in) :: molecules_per_particle(:,:) real(dp), intent(in) :: H2O_sat_mix(:), H2O_rh(:) - real(dp), intent(in) :: rainout_rates(:,:) + real(dp), intent(in) :: rainout_rates(:,:), scale_height(:), wfall(:,:) real(dp), intent(in) :: density(:) ${TYPE1}$, intent(inout) :: mix(:,:) ${TYPE1}$, intent(inout) :: densities(:,:), xp(:), xl(:) ${TYPE1}$, intent(inout) :: rhs(:) ! neqs - real(dp) :: H2O_cold_trap - ${TYPE1}$ :: dn_gas_dt, dn_particle_dt, cond_rate + real(dp) :: H2O_cold_trap, cond_rate0 + ${TYPE1}$ :: rh, dn_gas_dt, dn_particle_dt, cond_rate integer :: i, ii, j, k, kk type(PhotochemData), pointer :: dat @@ -59,6 +59,7 @@ subroutine dochem_${NAME}$(self, usol, rx_rates, & var => self%var #:if NAME == 'dual' + rh = dual(size(usol(1,1)%der)) dn_gas_dt = dual(size(usol(1,1)%der)) dn_particle_dt = dual(size(usol(1,1)%der)) cond_rate = dual(size(usol(1,1)%der)) @@ -144,6 +145,12 @@ subroutine dochem_${NAME}$(self, usol, rx_rates, & rhs(k) = rhs(k) + (xp(j) - xl(j)) enddo enddo + + ! 4 numbers for each condensing species. + ! condensation rate coefficient ~ 100.0 + ! evaporation rate coefficient ~ 10.0 + ! RH of condensation ~ 1.0 + ! RH smoothing factor ~ 0.1 ! particle condensation do j = 1,var%nz @@ -153,30 +160,48 @@ subroutine dochem_${NAME}$(self, usol, rx_rates, & ii = dat%particle_gas_phase_ind(i) ! index of gas phase kk = ii + (j - 1) * dat%nq ! gas phase rhs index k = i + (j - 1) * dat%nq ! particle rhs index + + ! compute the relative humidity + rh = max(usol(ii,j)/gas_sat_den(i,j),small_real) + + if (rh > var%cond_params(i)%RHc) then + ! Condensation occurs + + ! Condensation rate is based on how rapidly gases are mixing + ! in the atmosphere. We also smooth the rate to prevent stiffness. + cond_rate0 = var%cond_params(i)%k_cond*(var%edd(j)/scale_height(j)**2.0_dp) + cond_rate = damp_condensation_rate(cond_rate0, & + var%cond_params(i)%RHc, & + (1.0_dp + var%cond_params(i)%smooth_factor)*var%cond_params(i)%RHc, & + rh) - ! if the gas phase is super-saturated - if (densities(ii,j) >= gas_sat_den(i,j)*var%condensation_rate(2,i)) then - ! compute condensation rate - cond_rate = damp_condensation_rate(var%condensation_rate(1,i), & - var%condensation_rate(2,i), & - var%condensation_rate(3,i), & - densities(ii,j)/gas_sat_den(i,j)) - - ! rate the gas molecules are going to particles - ! in molecules/cm3/s - dn_gas_dt = - cond_rate* & - (densities(ii,j) - gas_sat_den(i,j)*var%condensation_rate(2,i)) - ! add to rhs vector, convert to change in mixing ratio/s + ! Rate gases are being destroyed (molecules/cm^3/s) + dn_gas_dt = - cond_rate*usol(ii,j) rhs(kk) = rhs(kk) + dn_gas_dt - ! rate of particle production from gas condensing - ! in particles/cm3/s + ! Rate particles are being produced (molecules/cm^3/s) dn_particle_dt = - dn_gas_dt - ! add to rhs vector, convert to moles/cm3/s rhs(k) = rhs(k) + dn_particle_dt - else - ! particles don't change! - rhs(k) = rhs(k) + 0.0_dp + + elseif (rh <= var%cond_params(i)%RHc .and. var%evaporation) then + ! Evaporation occurs + + ! Evaporation rate is based on how rapidly particles are falling + ! in the atmosphere. We also smooth the rate to prevent stiffness. + cond_rate0 = var%cond_params(i)%k_evap*(wfall(i,j)/scale_height(j)) + cond_rate = damp_condensation_rate(cond_rate0, & + 1.0_dp/var%cond_params(i)%RHc, & + (1.0_dp + var%cond_params(i)%smooth_factor)/var%cond_params(i)%RHc, & + 1.0_dp/rh) + + ! Rate gases are being produced (molecules/cm^3/s) + dn_gas_dt = cond_rate*usol(i,j) + rhs(kk) = rhs(kk) + dn_gas_dt + + ! Rate particles are being destroyed (molecules/cm^3/s) + dn_particle_dt = - dn_gas_dt + rhs(k) = rhs(k) + dn_particle_dt + endif endif enddo @@ -289,39 +314,51 @@ subroutine diffusion_coefficients_evo(dat, var, den, mubar, & ADD(i,j) = ADL(i,j) enddo - ! particles (eddy diffusion) - ! middle - do j = 2,var%nz-1 - do i = 1,dat%np - ! diffusion - DU(i,j) = eddav(j)/(var%dz(j)**2.0_dp) - DL(i,j) = eddav(j-1)/(var%dz(j)**2.0_dp) - DD(i,j) = - DU(i,j) - DL(i,j) + ! ! particles (eddy diffusion) + ! ! middle + ! do j = 2,var%nz-1 + ! do i = 1,dat%np + ! ! diffusion + ! DU(i,j) = eddav(j)/(var%dz(j)**2.0_dp) + ! DL(i,j) = eddav(j-1)/(var%dz(j)**2.0_dp) + ! DD(i,j) = - DU(i,j) - DL(i,j) + + ! ! advection + ! ADU(i,j) = gamma_i_part_av(i,j)/(2.0_dp*var%dz(j)) + ! ADL(i,j) = - gamma_i_part_av(i,j-1)/(2.0_dp*var%dz(j)) + ! ADD(i,j) = ADU(i,j) + ADL(i,j) + ! enddo + ! enddo + ! ! lower boundary + ! j = 1 + ! do i = 1,dat%np + ! DU(i,j) = eddav(j)/(var%dz(j)**2.0_dp) + ! DD(i,j) = - DU(i,j) + + ! ADU(i,j) = gamma_i_part_av(i,j)/(2.0_dp*var%dz(j)) + ! ADD(i,j) = ADU(i,j) + ! enddo + ! ! upper boundary + ! j = var%nz + ! do i = 1,dat%np + ! DL(i,j) = eddav(j-1)/(var%dz(j)**2.0_dp) + ! DD(i,j) = - DL(i,j) + + ! ADL(i,j) = - gamma_i_part_av(i,j-1)/(2.0_dp*var%dz(j)) + ! ADD(i,j) = ADL(i,j) + ! enddo - ! advection - ADU(i,j) = gamma_i_part_av(i,j)/(2.0_dp*var%dz(j)) - ADL(i,j) = - gamma_i_part_av(i,j-1)/(2.0_dp*var%dz(j)) - ADD(i,j) = ADU(i,j) + ADL(i,j) + ! We do not include eddy diffusion for particles + do j = 1,var%nz + do i = 1,dat%np + DU(i,j) = 0.0_dp + DL(i,j) = 0.0_dp + DD(i,j) = 0.0_dp + ADU(i,j) = 0.0_dp + ADL(i,j) = 0.0_dp + ADD(i,j) = 0.0_dp enddo enddo - ! lower boundary - j = 1 - do i = 1,dat%np - DU(i,j) = eddav(j)/(var%dz(j)**2.0_dp) - DD(i,j) = - DU(i,j) - - ADU(i,j) = gamma_i_part_av(i,j)/(2.0_dp*var%dz(j)) - ADD(i,j) = ADU(i,j) - enddo - ! upper boundary - j = var%nz - do i = 1,dat%np - DL(i,j) = eddav(j-1)/(var%dz(j)**2.0_dp) - DD(i,j) = - DL(i,j) - - ADL(i,j) = - gamma_i_part_av(i,j-1)/(2.0_dp*var%dz(j)) - ADD(i,j) = ADL(i,j) - enddo ! particles (falling) ! middle @@ -578,8 +615,8 @@ module subroutine prep_all_evo_gas(self, usol_in, err) call diffusion_coefficients_evo(dat, var, wrk%density, wrk%mubar, & wrk%DU, wrk%DD, wrk%DL, wrk%ADU, wrk%ADL, wrk%ADD, wrk%wfall, wrk%VH2_esc, wrk%VH_esc) - wrk%surface_scale_height = (k_boltz*var%temperature(1)*N_avo)/(wrk%mubar(1)*var%grav(1)) - + wrk%scale_height = (k_boltz*var%temperature(:)*N_avo)/(wrk%mubar(:)*var%grav(:)) + !!! H and H2 escape wrk%upper_veff_copy = var%upper_veff wrk%lower_vdep_copy = var%lower_vdep @@ -678,7 +715,7 @@ module subroutine rhs_evo_gas(self, neqs, tn, usol_flat, rhs, err) call dochem(self, wrk%usol, wrk%rx_rates, & wrk%gas_sat_den, wrk%molecules_per_particle, & - wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, & + wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%scale_height, wrk%wfall, & wrk%density, wrk%mix, wrk%densities, wrk%xp, wrk%xl, rhs) ! Extra functions specifying production or destruction @@ -722,7 +759,7 @@ module subroutine rhs_evo_gas(self, neqs, tn, usol_flat, rhs, err) elseif (var%lowerboundcond(i) == MosesBC) then rhs(i) = rhs(i) + wrk%DU(i,1)*wrk%usol(i,2) + wrk%ADU(i,1)*wrk%usol(i,2) & + wrk%DD(i,1)*wrk%usol(i,1) + wrk%ADD(i,1)*wrk%usol(i,1) & - - (var%edd(1)/wrk%surface_scale_height)*wrk%usol(i,1)/var%dz(1) + - (var%edd(1)/wrk%scale_height(1))*wrk%usol(i,1)/var%dz(1) endif enddo @@ -811,7 +848,7 @@ subroutine fcn(x_, f_) call initialize_dual_array(f_, blocksize) call dochem(self, usol_, wrk%rx_rates, & wrk%gas_sat_den, wrk%molecules_per_particle, & - wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, & + wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%scale_height, wrk%wfall, & wrk%density, mix, densities, xp, xl, f_) end subroutine end subroutine @@ -868,7 +905,7 @@ module subroutine jac_evo_gas(self, lda_neqs, neqs, usol_flat, jac, err) ! compute chemistry contribution to jacobian using forward differences call dochem(self, wrk%usol, wrk%rx_rates, & wrk%gas_sat_den, wrk%molecules_per_particle, & - wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, & + wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%scale_height, wrk%wfall, & wrk%density, mix, densities, xp, xl, rhs) !$omp parallel private(i, j, k, m, mm, usol_perturb, R, mix, densities, xl, xp, rhs_perturb) @@ -882,7 +919,7 @@ module subroutine jac_evo_gas(self, lda_neqs, neqs, usol_flat, jac, err) call dochem(self, usol_perturb, wrk%rx_rates, & wrk%gas_sat_den, wrk%molecules_per_particle, & - wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, & + wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%scale_height, wrk%wfall, & wrk%density, mix, densities, xp, xl, rhs_perturb) do m = 1,dat%nq @@ -951,7 +988,7 @@ module subroutine jac_evo_gas(self, lda_neqs, neqs, usol_flat, jac, err) elseif (var%lowerboundcond(i) == MosesBC) then djac(dat%ku,i+dat%nq) = wrk%DU(i,1) + wrk%ADU(i,1) djac(dat%kd,i) = djac(dat%kd,i) + wrk%DD(i,1) + wrk%ADD(i,1) - & - (var%edd(1)/wrk%surface_scale_height)/var%dz(1) + (var%edd(1)/wrk%scale_height(1))/var%dz(1) endif enddo @@ -1012,7 +1049,7 @@ module subroutine right_hand_side_chem(self, usol, rhs, err) call dochem(self, wrk%usol, wrk%rx_rates, & wrk%gas_sat_den, wrk%molecules_per_particle, & - wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, & + wrk%H2O_sat_mix, wrk%H2O_rh, wrk%rainout_rates, wrk%scale_height, wrk%wfall, & wrk%density, wrk%mix, wrk%densities, wrk%xp, wrk%xl, rhs) end subroutine diff --git a/src/input/photochem_input_read.f90 b/src/input/photochem_input_read.f90 index 6f254db..df2ed3b 100644 --- a/src/input/photochem_input_read.f90 +++ b/src/input/photochem_input_read.f90 @@ -633,7 +633,6 @@ subroutine unpack_settings(infile, s, dat, var, err) character(:), allocatable, intent(out) :: err integer :: j, i, ind(1), io - logical, allocatable :: particle_checklist(:) !!!!!!!!!!!!!!!!!!!!!!! !!! atmosphere-grid !!! @@ -760,6 +759,15 @@ subroutine unpack_settings(infile, s, dat, var, err) return endif endif + if ((dat%fix_water_in_trop .or. dat%water_cond) .and. dat%there_are_particles) then + ind = findloc(dat%particle_gas_phase,'H2O') + if (ind(1) /= 0) then + err = 'IOError: Either "fix-water-in-troposphere" or "water-condensation" is turned on'// & + ' in the settings file, but the reaction mechanism already implements H2O condensation'// & + ' via a particle.' + return + endif + endif if (dat%fix_water_in_trop) then @@ -800,52 +808,8 @@ subroutine unpack_settings(infile, s, dat, var, err) !!!!!!!!!!!!!!!!! !!! particles !!! !!!!!!!!!!!!!!!!! - if (dat%there_are_particles) then - if (any(dat%particle_formation_method == CondensingParticle)) then - ! then we need rate data - - if (.not. allocated(s%con_names)) then - err = 'settings file "'//trim(infile)//'" does not contain'// & - ' any particle condensation rate data.' - return - endif - - allocate(particle_checklist(dat%np)) - allocate(var%condensation_rate(3,dat%np)) - particle_checklist = .false. - - do i = 1,size(s%con_names) - - ind = findloc(dat%particle_names,trim(s%con_names(i))) - if (particle_checklist(ind(1))) then - err = "IOError: particle "//trim(s%con_names(i))//" in the settings"// & - " file is listed more than once" - return - endif - if (ind(1) == 0) then - err = "IOError: particle "//trim(s%con_names(i))//" in the settings"// & - " file isn't in the list of particles in the reaction mechanism file" - return - else - particle_checklist(ind(1)) = .true. - endif - - var%condensation_rate(1,ind(1)) = s%con(i)%A - var%condensation_rate(2,ind(1)) = s%con(i)%rhc - var%condensation_rate(3,ind(1)) = s%con(i)%rh0 - - enddo - - do i = 1,dat%np - if (dat%particle_formation_method(i) == CondensingParticle .and. .not. particle_checklist(i)) then - err = 'IOError: Particle '//trim(dat%particle_names(i))// & - ' does not have any condensation rate data in the file '//trim(infile) - return - endif - enddo - - endif - endif + ! Condensation rate parameters. size is zero when there are no particles + allocate(var%cond_params(dat%np)) !!!!!!!!!!!!!!!!!!!!!!!!!!! !!! boundary-conditions !!! diff --git a/src/photochem_types.f90 b/src/photochem_types.f90 index 5d80f83..9fb9fa1 100644 --- a/src/photochem_types.f90 +++ b/src/photochem_types.f90 @@ -22,12 +22,6 @@ module photochem_types ! make a giant IO object !!! Settings !!! !!!!!!!!!!!!!!!! - type :: SettingsCondensationRate - real(dp) :: A - real(dp) :: rhc - real(dp) :: rh0 - end type - type :: SettingsBC integer :: bc_type real(dp) :: vel @@ -78,10 +72,6 @@ module photochem_types ! make a giant IO object real(dp) :: trop_alt real(dp) :: H2O_condensation_rate(3) - ! particles - character(s_str_len), allocatable :: con_names(:) - type(SettingsCondensationRate), allocatable :: con(:) - ! boundary-conditions type(SettingsBC), allocatable :: ubcs(:) type(SettingsBC), allocatable :: lbcs(:) @@ -380,6 +370,15 @@ subroutine time_dependent_rate_fcn(tn, nz, rate) integer :: LH !! H index end type + + !> Condensation parameters + type :: CondensationParameters + real(dp) :: k_cond = 100.0_dp !! rate coefficient for condensation + real(dp) :: k_evap = 10.0_dp !! rate coefficient for evaporation + real(dp) :: RHc = 1.0_dp !! RH where condensation occurs + real(dp) :: smooth_factor = 0.1_dp !! A factor that smooths condensation/evaporation + !! rate to prevents stiffness + end type type :: PhotochemVars ! PhotochemVars contains information that can change between @@ -434,8 +433,10 @@ subroutine time_dependent_rate_fcn(tn, nz, rate) procedure(time_dependent_flux_fcn), nopass, pointer :: photon_flux_fcn => null() ! particles - ! condensation rate of particles - real(dp), allocatable :: condensation_rate(:,:) ! (3,np) + !> Parameters describing condensation and evaporation rates and + !> the RH needed for condensation + type(CondensationParameters), allocatable :: cond_params(:) ! (np) + logical :: evaporation = .true. !! If true, then evaporation occurs. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! set AFTER file read-in !!! @@ -578,7 +579,7 @@ subroutine time_dependent_rate_fcn(tn, nz, rate) real(dp) :: VH_esc ! other real(dp), allocatable :: sum_usol(:) !! (nz) - real(dp) :: surface_scale_height + real(dp), allocatable :: scale_height(:) real(dp), allocatable :: wfall(:,:) real(dp), allocatable :: gas_sat_den(:,:) real(dp), allocatable :: molecules_per_particle(:,:) @@ -657,6 +658,7 @@ subroutine init_PhotochemWrk(self, nsp, np, nq, nz, nrT, kj, nw) deallocate(self%upper_veff_copy) deallocate(self%lower_vdep_copy) deallocate(self%sum_usol) + deallocate(self%scale_height) deallocate(self%wfall) deallocate(self%gas_sat_den) deallocate(self%molecules_per_particle) @@ -690,6 +692,7 @@ subroutine init_PhotochemWrk(self, nsp, np, nq, nz, nrT, kj, nw) allocate(self%upper_veff_copy(nq)) allocate(self%lower_vdep_copy(nq)) allocate(self%sum_usol(nz)) + allocate(self%scale_height(nz)) allocate(self%wfall(np,nz)) allocate(self%gas_sat_den(np,nz)) allocate(self%molecules_per_particle(np,nz)) diff --git a/src/photochem_types_create.f90 b/src/photochem_types_create.f90 index a701a50..4a54a0f 100644 --- a/src/photochem_types_create.f90 +++ b/src/photochem_types_create.f90 @@ -276,49 +276,6 @@ function unpack_PhotoSettings(root, filename, err) result(s) endif - !!!!!!!!!!!!!!!!! - !!! particles !!! - !!!!!!!!!!!!!!!!! - - list => root%get_list('particles', .false., error = io_err) - - if (associated(list)) then - allocate(s%con_names(list%size())) - allocate(s%con(size(s%con_names))) - - j = 1 - item => list%first - do while (associated(item)) - select type (e => item%node) - class is (type_dictionary) - s%con_names(j) = trim(e%get_string('name',error = io_err)) - if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif - - tmp2 => e%get_dictionary('condensation-rate',.true.,error = io_err) - if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif - - s%con(j)%A = tmp2%get_real('A',error = io_err) - if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif - s%con(j)%rhc = tmp2%get_real('rhc',error = io_err) - if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif - s%con(j)%rh0 = tmp2%get_real('rh0',error = io_err) - if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif - if (s%con(j)%rh0 <= s%con(j)%rhc) then - err = 'IOError: Rate constant "rh0" for '//trim(s%con_names(j)) & - //' condensation must be > "rhc". See '//trim(filename) - return - endif - - class default - err = "IOError: Particle settings must be a list of dictionaries." - return - end select - j = j + 1 - item => item%next - end do - - endif - !!!!!!!!!!!!!!!!!!!!!!!!!!! !!! boundary-conditions !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!