00001
00002
00003
00004
00005 !> @file
00006 !!
00007 !! Computation of the evaporation rate
00008 !! (buoyant-diffusion approach by Ingersoll).
00009 !!
00010 !! @section Copyright
00011 !!
00012 !! Copyright 2010, 2011 Ralf Greve, Bjoern Grieger, Oliver J. Stenzel
00013 !!
00014 !! @section License
00015 !!
00016 !! This file is part of MAIC-2.
00017 !!
00018 !! MAIC-2 is free software: you can redistribute it and/or modify
00019 !! it under the terms of the GNU General Public License as published by
00020 !! the Free Software Foundation, either version 3 of the License, or
00021 !! (at your option) any later version.
00022 !!
00023 !! MAIC-2 is distributed in the hope that it will be useful,
00024 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
00025 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00026 !! GNU General Public License for more details.
00027 !!
00028 !! You should have received a copy of the GNU General Public License
00029 !! along with MAIC-2. If not, see <http://www.gnu.org/licenses/>.
00030 !<
00031 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00032
00033
00034 !> Computation of the evaporation rate
00035 !! (buoyant-diffusion approach by Ingersoll).
00036 !<------------------------------------------------------------------------------
00037 module evaporation
00038
00039 use maic2_types
00040
00041 implicit none
00042 real(dp), private :: evap_fact, gamma_reg,
00043 R_univ, mol_w, mol_c, diff_w_c, visc_c, g
00044
00045 contains
00046
00047
00048 !> Setting of parameters.
00049 !<------------------------------------------------------------------------------
00050 subroutine setevappar_bd(evap_fact_para, gamma_reg_para, &
00051 R_univ_para, mol_w_para, mol_c_para, &
00052 diff_w_c_para, visc_c_para, g_para)
00053
00054 implicit none
00055
00056 real(dp), optional :: evap_fact_para, gamma_reg_para,
00057 R_univ_para, mol_w_para, mol_c_para,
00058 diff_w_c_para, visc_c_para, g_para
00059
00060 if ( present(evap_fact_para) ) then
00061 evap_fact = evap_fact_para
00062 else
00063 evap_fact = 1.0_dp
00064 end if
00065
00066 if ( present(gamma_reg_para) ) then
00067 gamma_reg = gamma_reg_para
00068 else
00069 gamma_reg = 1.11e+11_dp
00070 end if
00071
00072 if ( present(R_univ_para) ) then
00073 R_univ = R_univ_para
00074 else
00075 R_univ = 8.314_dp
00076 end if
00077
00078 if ( present(mol_w_para) ) then
00079 mol_w = mol_w_para
00080 else
00081 mol_w = 1.802e-02_dp
00082 end if
00083
00084 if ( present(mol_c_para) ) then
00085 mol_c = mol_c_para
00086 else
00087 mol_c = 4.401e-02_dp
00088 end if
00089
00090 if ( present(diff_w_c_para) ) then
00091 diff_w_c = diff_w_c_para
00092 else
00093 diff_w_c = 1.4e-03_dp
00094 end if
00095
00096 if ( present(visc_c_para) ) then
00097 visc_c = visc_c_para
00098 else
00099 visc_c = 6.93e-04_dp
00100 end if
00101
00102 if ( present(g_para) ) then
00103 g = g_para
00104 else
00105 g = 3.7_dp
00106 end if
00107
00108 end subroutine setevappar_bd
00109
00110
00111 !> Computation of evaporation.
00112 !<------------------------------------------------------------------------------
00113 subroutine getevap_bd(temp, temp_amp, p, H, evap)
00114
00115 implicit none
00116
00117 real(dp) :: temp(:), temp_amp(:), p(:), H(:), evap(:)
00118
00119 integer(i4b) :: i, hr, n
00120 real(dp) :: inv_visc_c_2, one_third, one_eighth
00121 real(dp) :: p_sat, delta_eta, rho,
00122 delta_rho_rho_1, delta_rho_rho_2, delta_rho_rho,
00123 inv_length, inv_gamma_reg,
00124 temp_hr, evap_hr
00125 real(dp) :: sat_pressure
00126 real(dp), dimension(8) :: cos_factor
00127 real(dp), parameter :: pi = 3.141592653589793_dp
00128
00129 n = size(temp)
00130
00131 inv_visc_c_2 = 1.0_dp/visc_c**2
00132 inv_gamma_reg = 1.0_dp/gamma_reg
00133 one_third = 1.0_dp/3.0_dp
00134 one_eighth = 1.0_dp/8.0_dp
00135 cos_factor = (/(cos(2.0_dp*pi*one_eighth*real(hr,dp)),hr=1,8)/)
00136
00137 do i = 1, n
00138
00139 evap(i) = 0.0_dp
00140
00141 do hr=1, 8
00142
00143 temp_hr = temp(i) - temp_amp(i)*cos_factor(hr)
00144 sat_pressure = p_sat(temp_hr)
00145 delta_eta = mol_w*sat_pressure/(mol_c*p(i))
00146 rho = mol_c*p(i)/(R_univ*temp_hr)
00147
00148 delta_rho_rho_1 = (mol_c-mol_w)*sat_pressure
00149 delta_rho_rho_2 = mol_c*p(i)-(mol_c-mol_w)*sat_pressure
00150
00151 if (delta_rho_rho_1 <= delta_rho_rho_2) then
00152
00153 delta_rho_rho = delta_rho_rho_1/delta_rho_rho_2
00154 else
00155 delta_rho_rho = 1.0_dp
00156 end if
00157
00158 inv_length = (delta_rho_rho*g*inv_visc_c_2)**one_third
00159
00160 evap_hr = evap_fact * 0.17_dp*delta_eta*rho*diff_w_c*inv_length
00161
00162 evap(i) = evap(i) + one_eighth*evap_hr
00163
00164 end do
00165
00166 if (H(i) < 0.0_dp) evap(i) = evap(i) * exp(inv_gamma_reg*H(i))
00167
00168 end do
00169
00170 end subroutine getevap_bd
00171
00172 end module evaporation
00173