Changeset 3458 for palm/trunk/SOURCE
- Timestamp:
- Oct 30, 2018 2:51:23 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 10 edited
-
Makefile (modified) (2 diffs)
-
chem_emissions_mod.f90 (modified) (42 diffs)
-
chem_modules.f90 (modified) (4 diffs)
-
chemistry_model_mod.f90 (modified) (14 diffs)
-
date_and_time_mod.f90 (modified) (11 diffs)
-
flow_statistics.f90 (modified) (2 diffs)
-
init_3d_model.f90 (modified) (2 diffs)
-
netcdf_data_input_mod.f90 (modified) (4 diffs)
-
palm.f90 (modified) (3 diffs)
-
prognostic_equations.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r3448 r3458 25 25 # ----------------- 26 26 # $Id$ 27 # from chemistry branch r3443, banzhafs, Russo, forkel, basit: 28 # radiation_model_mod added to chemistry_model_mod 29 # Added some missing dependencies and replaced blanks with tabs 30 # Added chemistry emission module 31 # chemistry_model_mod added to flow_statistics 32 # 33 # 3448 2018-10-29 18:14:31Z kanani 27 34 # Adjustment of biometeorology dependencies 28 35 # … … 856 863 modules.o \ 857 864 netcdf_data_input_mod.o \ 865 radiation_model_mod.o \ 858 866 surface_mod.o \ 859 867 user_module.o -
palm/trunk/SOURCE/chem_emissions_mod.f90
r3373 r3458 27 27 ! ----------------- 28 28 ! $Id$ 29 ! from chemistry branch r3443, banzhafs, Russo: 30 ! Additional correction for index of input file of pre-processed mode 31 ! Removed atomic and molecular weights as now available in chem_modules and 32 ! added par_emis_time_factor (formerly in netcdf_data_input_mod) 33 ! Added correction for index of input file of pre-processed mode 34 ! Added a condition for default mode necessary for information read from ncdf file 35 ! in pre-processed and default mode 36 ! Correction of multiplying factor necessary for scaling emission values in time 37 ! Introduced correction for scaling units in the case of DEFAULT emission mode 38 ! 39 ! 3373 2018-10-18 15:25:56Z kanani 29 40 ! Fix wrong location of __netcdf directive 30 41 ! … … 40 51 ! Authors: 41 52 ! -------- 42 ! Emmanuele Russo Fu-Berlin43 ! Sabine Banzhaf FU-Berlin44 ! Martijn Schaap FU-Berlin, TNO Utrecht53 ! @author Emmanuele Russo (Fu-Berlin) 54 ! @author Sabine Banzhaf (FU-Berlin) 55 ! @author Martijn Schaap (FU-Berlin, TNO Utrecht) 45 56 ! 46 57 ! Description: … … 49 60 !> 50 61 !> @todo Check_parameters routine should be developed: for now it includes just one condition 51 !> @todo Add option for capital or small letters in the matching routine52 62 !> @todo Use of Restart files not contempled at the moment 63 !> @todo revise indices of files read from the netcdf: preproc_emission_data and expert_emission_data 64 !> @todo for now emission data may be passed on a singular vertical level: need to be more flexible 53 65 !> @note <Enter notes on the module> 54 66 !> @bug <Enter known bugs here> … … 61 73 62 74 USE control_parameters, & 63 ONLY: initializing_actions, end_time, message_string,&75 ONLY: initializing_actions, end_time, message_string, & 64 76 intermediate_timestep_count, dt_3d 65 77 … … 67 79 68 80 USE kinds 81 82 #if defined ( __netcdf ) 69 83 70 84 USE netcdf_data_input_mod, & 71 85 ONLY: chem_emis_att_type, chem_emis_val_type 72 86 73 #if defined ( __netcdf )74 87 USE NETCDF 88 75 89 #endif 76 90 77 91 USE date_and_time_mod, & 78 ONLY: time_default_indices, month_of_year, day_of_month, hour_of_day, & 92 ONLY: time_default_indices, time_preprocessed_indices, & 93 month_of_year, day_of_month, hour_of_day, & 79 94 index_mm, index_dd, index_hh 80 95 … … 92 107 93 108 CHARACTER (LEN=80) :: filename_emis !> Variable for the name of the netcdf input file 94 CHARACTER (LEN=80), ALLOCATABLE, DIMENSION(:) :: string_values !> Output of string variables read from netcdf input file 95 96 INTEGER(iwp) :: n_dims !> Number of dimensions of input variable netcdf file 109 97 110 INTEGER(iwp) :: i !> index 1st selected dimension (some dims are not spatial) 98 111 INTEGER(iwp) :: j !> index 2nd selected dimension … … 103 116 INTEGER(iwp) :: z_start !> Index to start read variable from netcdf in additional dims 104 117 INTEGER(iwp) :: z_end !> Index to end read variable from netcdf in additional dims 105 INTEGER(iwp), ALLOCATABLE, DIMENSION(:) :: id_dims !> id of dimension of selected variable netcdf input106 INTEGER(iwp), ALLOCATABLE, DIMENSION(:) :: dum_var !> value of variable read from netcdf input107 INTEGER(iwp) :: id_ncfile !> id netcdf file108 INTEGER(iwp) :: errno !> error number NF90_???? function109 INTEGER(iwp) :: id_var !< variable id110 118 INTEGER(iwp) :: dt_emis !> Time Step Emissions 111 119 INTEGER(iwp) :: len_index !> length of index (used for several indices) … … 124 132 !> Dobson units: 125 133 REAL, PARAMETER :: Dobs = 2.68668e16 ! (mlc/cm2) / DU 126 127 ! molar weights of components128 129 !> atomic weights:130 REAL, PARAMETER :: xm_H = 1.00790e-3 ! kg/mol131 REAL, PARAMETER :: xm_N = 14.00670e-3 ! kg/mol132 REAL, PARAMETER :: xm_C = 12.01115e-3 ! kg/mol133 REAL, PARAMETER :: xm_S = 32.06400e-3 ! kg/mol134 REAL, PARAMETER :: xm_O = 15.99940e-3 ! kg/mol135 REAL, PARAMETER :: xm_F = 18.99840e-3 ! kg/mol136 REAL, PARAMETER :: xm_Na = 22.98977e-3 ! kg/mol137 REAL, PARAMETER :: xm_Cl = 35.45300e-3 ! kg/mol138 REAL, PARAMETER :: xm_Rn222 = 222.00000e-3 ! kg/mol139 REAL, PARAMETER :: xm_Pb210 = 210.00000e-3 ! kg/mol140 REAL, PARAMETER :: xm_Ca = 40.07800e-3 ! kg/mol141 REAL, PARAMETER :: xm_K = 39.09800e-3 ! kg/mol142 REAL, PARAMETER :: xm_Mg = 24.30500e-3 ! kg/mol143 REAL, PARAMETER :: xm_Pb = 207.20000e-3 ! kg/mol144 REAL, PARAMETER :: xm_Cd = 112.41000e-3 ! kg/mol145 REAL, PARAMETER :: xm_Rh = 102.90550e-3 ! kg/mol146 147 !> molecule weights:148 REAL, PARAMETER :: xm_h2o = xm_H * 2 + xm_O ! kg/mol149 REAL, PARAMETER :: xm_o3 = xm_O * 3 ! kg/mol150 REAL, PARAMETER :: xm_N2O5 = xm_N * 2 + xm_O * 5 ! kg/mol151 REAL, PARAMETER :: xm_HNO3 = xm_H + xm_N + xm_O * 3 ! kg/mol152 REAL, PARAMETER :: xm_NH4 = xm_N + xm_H * 4 ! kg/mol153 REAL, PARAMETER :: xm_SO4 = xm_S + xm_O * 4 ! kg/mol154 REAL, PARAMETER :: xm_NO3 = xm_N + xm_O * 3 ! kg/mol155 REAL, PARAMETER :: xm_CO2 = xm_C + xm_O * 2 ! kg/mol156 157 !> mass of air158 REAL, PARAMETER :: xm_air = 28.964e-3 ! kg/mol159 REAL, PARAMETER :: xmair = 28.94 ! g/mol; old name!160 134 161 135 !> sesalt composition: … … 263 237 !-- Public Variables 264 238 265 PUBLIC id_ncfile,errno,id_dims,n_dims,dum_var, con_factor, &266 len_index,len_index_voc,len_index_pm, string_values 239 PUBLIC con_factor, len_index,len_index_voc,len_index_pm 240 267 241 CONTAINS 268 242 … … 275 249 276 250 277 !TBD: Where Should we put the call to chem_emissions_check_parameters? In chem_init or in check_parameters?251 !TBD: Where Should we put the call to chem_emissions_check_parameters? In chem_init or in check_parameters? 278 252 279 253 IMPLICIT NONE … … 548 522 IF (len_index_voc>0) THEN 549 523 ALLOCATE(match_spec_voc_model(len_index_voc)) !> contains indices of the VOC model species 550 ALLOCATE(match_spec_voc_input(len_index_voc)) !> In input there is only one value for VOCs in the DEFAULT mode. This array contains the indices of the different values of VOC compositions of the input variable VOC_composition 524 ALLOCATE(match_spec_voc_input(len_index_voc)) !> In input there is only one value for VOCs in the DEFAULT mode. 525 ! This array contains the indices of the different values of VOC compositions of the input variable VOC_composition 551 526 ENDIF 552 527 … … 661 636 ' DO NOT MATCH' // & 662 637 ' model chemical species' // & 663 ' Chemistry Emissions routine is not called' !TBD: IMPORTANT:Add a condition in chem init that sets the emission_output_required to false when len_index==0 so that the chem_emissions modules are not used638 ' Chemistry Emissions routine is not called' 664 639 CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 ) 665 640 … … 724 699 DO ispec = 1 , len_index 725 700 726 IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND. & 727 emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN 728 729 message_string = 'PARAMETERIZED emissions mode selected:' // & 701 IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND. emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN 702 703 message_string = 'PARAMETERIZED emissions mode selected:' // & 730 704 ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES' // & 731 705 ' but values of scaling factors for street types' // & … … 780 754 SUBROUTINE chem_emissions_init(emt_att,emt,nspec_out) 781 755 782 756 #if defined( __netcdf ) 783 757 784 758 USE surface_mod, & … … 799 773 800 774 INTEGER(iwp) :: ispec !> Index to go through the emission chemical species 801 #if defined( __netcdf ) 775 802 776 803 777 !-- Actions for initial runs : TBD: needs to be updated … … 854 828 SELECT CASE(TRIM(mode_emis)) !TBD: Add the option for CApital or small letters 855 829 830 831 !> PRE-PROCESSED case 832 CASE ("PRE-PROCESSED") 833 834 IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out)) 835 836 CALL location_message( 'emis_distribution array allocated in PRE-PROCESSED mode', .FALSE. ) 837 838 !> Calculate the values of the emissions at the first time step 839 CALL chem_emissions_setup(emt_att,emt,nspec_out) 840 856 841 !> Default case 857 842 CASE ("DEFAULT") … … 864 849 CALL chem_emissions_setup(emt_att,emt,nspec_out) 865 850 866 867 !> PRE-PROCESSED case 868 CASE ("PRE-PROCESSED") 869 870 IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out)) 871 872 CALL location_message( 'emis_distribution array allocated in PRE-PROCESSED mode', .FALSE. ) 873 874 !> Calculate the values of the emissions at the first time step 875 CALL chem_emissions_setup(emt_att,emt,nspec_out) 876 851 !> PARAMETERIZED case 877 852 CASE ("PARAMETERIZED") 878 853 879 854 CALL location_message( 'emis_distribution array allocated in PARAMETERIZED mode', .FALSE. ) 880 855 881 ! For now for PAR and DEF values only, first vertical level of emis_distribution is allocated, while for EXPall.856 ! For now for PAR and DEF values only, first vertical level of emis_distribution is allocated, while for PRE-PROCESSED all. 882 857 IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(1,0:ny,0:nx,nspec_out)) 883 858 … … 919 894 IMPLICIT NONE 920 895 921 896 #if defined( __netcdf ) 922 897 923 898 !--- IN/OUT … … 949 924 REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: emis 950 925 926 REAL(wp), DIMENSION(24) :: par_emis_time_factor !< time factors 927 ! for the parameterized mode: these are fixed for each hour 928 ! of a single day. 951 929 REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: conv_to_ratio !> factor used for converting input 952 930 ! to adimensional concentration ratio … … 964 942 ! --- const ------------------------------- 965 943 !-CONVERSION FACTORS: TIME 944 ! number of sec per hour = 3600 945 REAL, PARAMETER :: s_per_hour = 3600.0 ! (s)/(hour) 946 ! number of sec per day = 86400 947 REAL, PARAMETER :: s_per_day = 86400.0 ! (s)/(day) 966 948 ! number of hours in a year of 365 days: 967 REAL, PARAMETER :: hour_per_year = 8760.0 !> TBD: What about leap years? 949 REAL, PARAMETER :: hour_per_year = 8760.0 !> TBD: What about leap years? 950 ! number of hours in a day: 951 REAL, PARAMETER :: hour_per_day = 24.0 952 968 953 ! conversion from hours to seconds (s/hour) = 1/3600.0 ~ 0.2777778 969 REAL, PARAMETER :: hour_to_s = 0.0002777778! (hour)/(s)954 REAL, PARAMETER :: hour_to_s = 1/s_per_hour ! (hour)/(s) 970 955 ! conversion from day to seconds (s/day) = 1/86400.0 ~ 1.157407e-05 971 REAL, PARAMETER :: day_to_s = 1 .157407e-05! (day)/(s)956 REAL, PARAMETER :: day_to_s = 1/s_per_day ! (day)/(s) 972 957 ! conversion from year to sec (s/year) = 1/31536000.0 ~ 3.170979e-08 973 REAL, PARAMETER :: year_to_s = 3.170979e-08! (year)/(s)958 REAL, PARAMETER :: year_to_s = 1/(s_per_hour*hour_per_year) ! (year)/(s) 974 959 975 960 !-CONVERSION FACTORS: WEIGHT … … 984 969 REAL(wp), PARAMETER :: ratio2ppm = 1.0e06_wp 985 970 !------------------------------------------------------ 986 #if defined( __netcdf ) 971 987 972 IF ( emission_output_required ) THEN 988 973 989 !> Set emis_dt !TBD: this is the same as dt_chem. We should consider the fact that dt_emis should be the timestep of input emissions or better defined, the timestep at which the emission routines are calle : for now one hour. It should be made changeable.974 !> Set emis_dt !TBD: this is the same as dt_chem. We should consider the fact that dt_emis should be the timestep of input emissions or better defined, the timestep at which the emission routines are called: for now one hour. It should be made changeable. 990 975 991 976 IF ( call_chem_at_all_substeps ) THEN … … 1093 1078 1094 1079 message_string = 'No Units conversion required for units of chemistry emissions' // & 1095 ' of the PARAMETERIZED mode: units have to be provided in mole/m**2/s ' 1080 ' of the PARAMETERIZED mode: units have to be provided in' // & 1081 ' micromole/m**2/day for GASES and' // & 1082 ' kg/m**2/day for PMs' 1096 1083 CALL message( 'chem_emissions_setup', 'CM0447', 0, 0, 0, 6, 0 ) 1097 1084 … … 1110 1097 ! V/N=RT/P 1111 1098 1112 !> m**3/Nmole (J/mol)*K^-1 K Pa1099 !> m**3/Nmole (J/mol)*K^-1 K Pa 1113 1100 conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) ) 1114 1101 ENDDO … … 1125 1112 !> PRE-PROCESSED MODE 1126 1113 IF (TRIM(mode_emis)=="PRE-PROCESSED") THEN 1114 1115 !> Update time indices 1116 CALL time_preprocessed_indices(index_hh) 1127 1117 1128 1118 CALL location_message( 'PRE-PROCESSED MODE: No time-factor specification required', .FALSE. ) … … 1194 1184 ' 24 values for every day of the year ', .FALSE. ) 1195 1185 1186 !Assign Constant Values of time factors, diurnal time profile for traffic sector: 1187 par_emis_time_factor( : ) = (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039, 0.056, 0.053, 0.051, 0.051, 0.052, 0.055, & 1188 0.059, 0.061, 0.064, 0.067, 0.069, 0.069, 0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /) 1189 1196 1190 !> in this case allocate time factor each hour in a day 1197 1191 IF (.NOT. ALLOCATED(time_factor)) ALLOCATE(time_factor(1)) … … 1200 1194 index_hh=hour_of_day 1201 1195 1202 time_factor(1) =emt_att%par_emis_time_factor(index_hh)1196 time_factor(1) = par_emis_time_factor(index_hh) 1203 1197 1204 1198 ENDIF … … 1212 1206 DO ispec=1,nspec_out 1213 1207 1214 !> Values are still micromoles/(m**2*s). Units are not given by the user for this mode and are always micromoles (or microgramsfor PMs)1215 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=surface_csflux(match_spec_input(ispec))*time_factor(1) 1208 !> Values are still micromoles/(m**2*s). Units in this case are always micromoles/m**2*day (or kilograms/m**2*day for PMs) 1209 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=surface_csflux(match_spec_input(ispec))*time_factor(1)*hour_to_s 1216 1210 1217 1211 ENDDO … … 1223 1217 DO ispec=1,nspec_out !> nspec_out represents the number of species in common between 1224 1218 ! the emission input data and the chemistry mechanism used 1225 1226 emis_distribution( nzb:nzt+1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))%&1227 preproc_emission_data(index_hh,nzb:nzt+1,nys:nyn,nxl:nxr)* &1228 con_factor1229 1219 1220 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))% & 1221 preproc_emission_data(index_hh,1,nys+1:nyn+1,nxl+1:nxr+1)* & 1222 con_factor 1223 1230 1224 ENDDO 1231 1225 … … 1249 1243 !TBD: The consideration of dt_emis of the input data is still missing. Basically the emissions could be provided every 10, 30 minutes and not always at one hour. This should be eventually solved in the date_and_time mode routine. 1250 1244 1245 !> the time factors are 24 for each day. When multiplied by a daily value, they allow to have an hourly value. Then to convert it to seconds, we still have to divide this value by 3600. 1246 ! So given any units, we convert them to seconds and finally multiply them by 24 ((value/sec)*(24*3600)=value/day ---- (value/day)*time_factor=value/hour ---(value/hour)/(3600)=value/sec ) 1247 ! ((value/sec)*(24*3600)*time_factor)/3600=24*(value/sec)*time_factor 1248 1251 1249 !> NOX Compositions 1252 1250 IF (TRIM(spc_names(match_spec_model(ispec)))=="NO") THEN 1253 !> Kg/m2 kg/m2*s 1254 delta_emis(nys:nyn,nxl:nxr) = & 1255 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor 1251 !> Kg/m2*s kg/m2*s 1252 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor*hour_per_day 1256 1253 1257 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1258 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1254 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1259 1255 1260 1256 ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="NO2") THEN 1261 1257 1262 delta_emis(nys:nyn,nxl:nxr) = & 1263 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,2)*con_factor 1264 1265 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1266 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1258 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,2)*con_factor*hour_per_day 1259 1260 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1267 1261 1268 1262 !> SOX Compositions 1269 1263 1270 1264 ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO2") THEN 1271 1272 delta_emis(nys:nyn,nxl:nxr) = & 1273 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor 1274 1275 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1276 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1265 !> Kg/m2*s kg/m2*s 1266 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor*hour_per_day 1267 1268 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1277 1269 1278 1270 ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO4") THEN 1279 1280 delta_emis(nys:nyn,nxl:nxr) = & 1281 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor 1282 1283 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1284 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1271 !> Kg/m2*s kg/m2*s 1272 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor*hour_per_day 1273 1274 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1285 1275 1286 1276 … … 1292 1282 DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,1)) 1293 1283 1294 delta_emis(nys:nyn,nxl:nxr) = & 1295 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,1)*con_factor 1284 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,1)*con_factor*hour_per_day 1296 1285 1297 1286 1298 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1299 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1287 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1300 1288 1301 1289 ENDDO … … 1307 1295 DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,2)) 1308 1296 1309 delta_emis(nys:nyn,nxl:nxr) = & 1310 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,2)*con_factor 1297 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,2)*con_factor*hour_per_day 1311 1298 1312 1299 1313 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1314 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1300 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1315 1301 1316 1302 ENDDO … … 1322 1308 DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,3)) 1323 1309 1324 delta_emis(nys:nyn,nxl:nxr) = & 1325 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,3)*con_factor 1310 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,3)*con_factor*hour_per_day 1326 1311 1327 1312 1328 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1329 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1313 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1330 1314 1331 1315 ENDDO … … 1339 1323 IF (TRIM(spc_names(match_spec_model(ispec)))==TRIM(emt_att%voc_name(ivoc))) THEN 1340 1324 1341 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)* & 1342 emt_att%voc_comp(icat,match_spec_voc_input(ivoc))*con_factor 1343 1344 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1345 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1325 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)* & 1326 emt_att%voc_comp(icat,match_spec_voc_input(ivoc))*con_factor*hour_per_day 1327 1328 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1346 1329 1347 1330 ENDIF … … 1352 1335 ELSE 1353 1336 1354 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*con_factor 1355 1356 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1357 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1337 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*con_factor*hour_per_day 1338 1339 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1358 1340 1359 1341 ENDIF ! IF (spc_names==) … … 1396 1378 DO ispec=1,nspec_out 1397 1379 1398 !> PMs are already in mass units:micrograms:they have to be converted to kilograms 1399 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" & 1400 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1380 !> PMs are already in mass units: kilograms 1381 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1401 1382 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1402 1383 … … 1412 1393 !> Other Species: inputs are micromoles: they have to be converted in moles 1413 1394 ! ppm/s *m *kg/m^3 1414 surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))* &1395 surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))* & 1415 1396 ! micromoles/(m^2*s) 1416 emis_distribution(1,j,i,ispec) * &1417 ! m **3/Nmole1418 conv_to_ratio(k,j,i)* &1397 emis_distribution(1,j,i,ispec) * & 1398 ! m^3/Nmole 1399 conv_to_ratio(k,j,i)* & 1419 1400 ! kg/m^3 1420 1401 rho_air(k) … … 1426 1407 1427 1408 1428 ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. & 1429 street_type_f%var(j,i) < main_street_id ) THEN 1409 ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. street_type_f%var(j,i) < main_street_id ) THEN 1430 1410 1431 1411 !> Cycle over already matched species … … 1433 1413 1434 1414 !> PMs are already in mass units: micrograms 1435 IF ( TRIM(spc_names(match_spec_model(ispec)))=="PM1" & 1436 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1415 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1437 1416 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1438 1417 1439 1418 ! kg/(m^2*s) *kg/m^3 1440 surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) * &1419 surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) * & 1441 1420 ! kg/(m^2*s) 1442 1421 emis_distribution(1,j,i,ispec)* & … … 1448 1427 !>Other Species: inputs are micromoles 1449 1428 ! ppm/s *m *kg/m^3 1450 surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * &1429 surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * & 1451 1430 ! micromoles/(m^2*s) 1452 emis_distribution(1,j,i,ispec) * &1453 ! m **3/Nmole1454 conv_to_ratio(k,j,i)* &1431 emis_distribution(1,j,i,ispec) * & 1432 ! m^3/Nmole 1433 conv_to_ratio(k,j,i)* & 1455 1434 ! kg/m^3 1456 1435 rho_air(k) … … 1486 1465 j = surf_def_h(0)%j(m) 1487 1466 1488 !> Distinguish between PMs (no needing conversion in ppms), 1489 ! VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and 1490 ! other Species (still expressed in Kg/(m**2*s) at this point) 1491 1492 !> PMs 1493 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" & 1494 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1495 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1467 IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN 1468 1469 1470 !> Distinguish between PMs (no needing conversion in ppms), 1471 ! VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and 1472 ! other Species (still expressed in Kg/(m**2*s) at this point) 1473 1474 !> PMs 1475 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1476 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1496 1477 1497 ! kg/(m^2*s) *kg/m^3 kg/(m^2*s) 1498 surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)* & 1499 ! kg/m^3 1500 rho_air(k) 1501 1502 1503 ELSE 1504 1505 !> VOCs 1506 IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN 1507 ! ( ppm/s) * m * kg/m^3 mole/(m**2/s) 1508 surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1509 ! m**3/mole ppm 1510 conv_to_ratio(k,j,i)*ratio2ppm * & 1511 ! kg/m^3 1512 rho_air(k) 1513 1514 1515 !> OTHER SPECIES 1478 ! kg/(m^2*s) *kg/m^3 kg/(m^2*s) 1479 surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)* & 1480 ! kg/m^3 1481 rho_air(nzb) 1482 1483 1516 1484 ELSE 1517 1485 1518 ! ( ppm/s )*m * kg/m^3 kg/(m**2/s) 1519 surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1520 ! mole/Kg 1521 (1/emt_att%xm(ispec))* & 1522 ! m**3/mole ppm 1523 conv_to_ratio(k,j,i)*ratio2ppm* & 1524 ! kg/m^3 1525 rho_air(k) 1526 1486 !> VOCs 1487 IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN 1488 ! ( ppm/s) * m * kg/m^3 mole/(m^2/s) 1489 surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1490 ! m^3/mole ppm 1491 conv_to_ratio(nzb,j,i)*ratio2ppm * & 1492 ! kg/m^3 1493 rho_air(nzb) 1494 1495 1496 !> OTHER SPECIES 1497 ELSE 1498 1499 ! ( ppm/s )*m * kg/m^3 kg/(m^2/s) 1500 surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1501 ! mole/Kg 1502 (1/emt_att%xm(ispec))* & 1503 ! m^3/mole ppm 1504 conv_to_ratio(nzb,j,i)*ratio2ppm* & 1505 ! kg/m^3 1506 rho_air(nzb) 1507 1508 1509 ENDIF 1527 1510 1528 1511 ENDIF … … 1543 1526 k = surf_lsm_h%k(m) 1544 1527 1545 !> Distinguish between PMs (no needing conversion in ppms), 1546 ! VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and 1547 ! other Species (still expressed in Kg/(m**2*s) at this point) 1548 1549 !> PMs 1550 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" & 1551 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1552 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1553 1554 ! kg/(m^2*s) * kg/m^3 kg/(m^2*s) 1555 surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1556 ! kg/m^3 1557 rho_air(k) 1558 1559 ELSE 1560 1561 !> VOCs 1562 IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN 1563 ! ( ppm/s) * m * kg/m^3 mole/(m**2/s) 1564 surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1565 ! m**3/mole ppm 1566 conv_to_ratio(k,j,i)*ratio2ppm* & 1567 ! kg/m^3 1568 rho_air(k) 1569 1570 1571 !> OTHER SPECIES 1528 IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN 1529 1530 !> Distinguish between PMs (no needing conversion in ppms), 1531 ! VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and 1532 ! other Species (still expressed in Kg/(m**2*s) at this point) 1533 1534 !> PMs 1535 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1536 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1537 1538 ! kg/(m^2*s) * kg/m^3 kg/(m^2*s) 1539 surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1540 ! kg/m^3 1541 rho_air(k) 1542 1572 1543 ELSE 1573 1544 1574 ! ( ppm/s) * m * kg/m^3 kg/(m**2/s) 1575 surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1576 ! mole/Kg 1577 (1/emt_att%xm(ispec))* & 1578 ! m**3/mole ppm 1579 conv_to_ratio(k,j,i)*ratio2ppm* & 1580 ! kg/m^3 1581 rho_air(k) 1545 !> VOCs 1546 IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN 1547 ! ( ppm/s) * m * kg/m^3 mole/(m^2/s) 1548 surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1549 ! m^3/mole ppm 1550 conv_to_ratio(k,j,i)*ratio2ppm* & 1551 ! kg/m^3 1552 rho_air(k) 1553 1554 !> OTHER SPECIES 1555 ELSE 1556 1557 ! ( ppm/s) * m * kg/m^3 kg/(m^2/s) 1558 surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1559 ! mole/Kg 1560 (1/emt_att%xm(ispec))* & 1561 ! m^3/mole ppm 1562 conv_to_ratio(k,j,i)*ratio2ppm* & 1563 ! kg/m^3 1564 rho_air(k) 1582 1565 1566 ENDIF 1567 1583 1568 ENDIF 1584 1569 1585 1570 ENDIF 1571 1586 1572 ENDDO 1587 1573 … … 1598 1584 k = surf_usm_h%k(m) 1599 1585 1600 1601 !> PMs 1602 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" & 1603 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1604 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1605 1606 ! kg/(m^2*s) *kg/m^3 kg/(m^2*s) 1607 surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)* & 1608 ! kg/m^3 1609 rho_air(k) 1610 1611 1612 ELSE 1613 1614 !> VOCs 1615 IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN 1616 ! ( ppm/s) * m * kg/m^3 mole/(m**2/s) 1617 surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1618 ! m**3/mole ppm 1619 conv_to_ratio(k,j,i)*ratio2ppm* & 1620 ! kg/m^3 1586 IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN 1587 1588 !> PMs 1589 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1590 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1591 1592 ! kg/(m^2*s) *kg/m^3 kg/(m^2*s) 1593 surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)* & 1594 ! kg/m^3 1621 1595 rho_air(k) 1622 1596 1623 !> OTHER SPECIES1597 1624 1598 ELSE 1625 1599 1626 1627 ! ( ppm/s) * m * kg/m^3 kg/(m**2/s) 1628 surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1629 ! mole/Kg 1630 (1/emt_att%xm(ispec))* & 1631 ! m**3/mole ppm 1632 conv_to_ratio(k,j,i)*ratio2ppm* & 1633 ! kg/m^3 1634 rho_air(k) 1635 1600 !> VOCs 1601 IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN 1602 ! ( ppm/s) * m * kg/m^3 mole/(m^2/s) 1603 surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1604 ! m^3/mole ppm 1605 conv_to_ratio(k,j,i)*ratio2ppm* & 1606 ! kg/m^3 1607 rho_air(k) 1608 1609 !> OTHER SPECIES 1610 ELSE 1611 1612 1613 ! ( ppm/s) * m * kg/m^3 kg/(m^2/s) 1614 surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) * & 1615 ! mole/Kg 1616 (1/emt_att%xm(ispec))* & 1617 ! m^3/mole ppm 1618 conv_to_ratio(k,j,i)*ratio2ppm* & 1619 ! kg/m^3 1620 rho_air(k) 1621 1622 1623 ENDIF 1636 1624 1637 1625 ENDIF -
palm/trunk/SOURCE/chem_modules.f90
r3298 r3458 27 27 ! ----------------- 28 28 ! $Id$ 29 ! from chemistry branch r3443: 30 ! ?? 31 ! 32 ! 3298 2018-10-02 12:21:11Z kanani 29 33 ! - Minor formatting 30 34 ! - Introduced Variables for Chemistry emissions Module (Russo) … … 43 47 ! @author Farah Kanani-Suehring 44 48 ! @author Basit Khan 49 ! @author Sabine Banzhaf 50 ! @author Emmanuele Russo 45 51 ! 46 52 !------------------------------------------------------------------------------! … … 88 94 LOGICAL :: do_emis = .FALSE. !< Flag for turning on chemistry emissions 89 95 LOGICAL :: cs_pr_namelist_found = .FALSE. !< Namelist parameter: Names of t 96 LOGICAL :: do_depo = .FALSE. !< namelist parameter for activation of deposition calculation 90 97 91 98 … … 161 168 ! matching between the model and the input files 162 169 170 171 !-- Selected atomic/molecular weights: 172 173 REAL, PARAMETER :: xm_H = 1.00790e-3 ! kg/mol 174 REAL, PARAMETER :: xm_N = 14.00670e-3 ! kg/mol 175 REAL, PARAMETER :: xm_C = 12.01115e-3 ! kg/mol 176 REAL, PARAMETER :: xm_S = 32.06400e-3 ! kg/mol 177 REAL, PARAMETER :: xm_O = 15.99940e-3 ! kg/mol 178 REAL, PARAMETER :: xm_F = 18.99840e-3 ! kg/mol 179 REAL, PARAMETER :: xm_Na = 22.98977e-3 ! kg/mol 180 REAL, PARAMETER :: xm_Cl = 35.45300e-3 ! kg/mol 181 REAL, PARAMETER :: xm_Rn222 = 222.00000e-3 ! kg/mol 182 REAL, PARAMETER :: xm_Pb210 = 210.00000e-3 ! kg/mol 183 REAL, PARAMETER :: xm_Ca = 40.07800e-3 ! kg/mol 184 REAL, PARAMETER :: xm_K = 39.09800e-3 ! kg/mol 185 REAL, PARAMETER :: xm_Mg = 24.30500e-3 ! kg/mol 186 REAL, PARAMETER :: xm_Pb = 207.20000e-3 ! kg/mol 187 REAL, PARAMETER :: xm_Cd = 112.41000e-3 ! kg/mol 188 189 REAL, PARAMETER :: xm_h2o = xm_H * 2 + xm_O ! kg/mol 190 REAL, PARAMETER :: xm_o3 = xm_O * 3 ! kg/mol 191 REAL, PARAMETER :: xm_N2O5 = xm_N * 2 + xm_O * 5 ! kg/mol 192 REAL, PARAMETER :: xm_HNO3 = xm_H + xm_N + xm_O * 3 ! kg/mol 193 REAL, PARAMETER :: xm_NH4 = xm_N + xm_H * 4 ! kg/mol 194 REAL, PARAMETER :: xm_SO4 = xm_S + xm_O * 4 ! kg/mol 195 REAL, PARAMETER :: xm_NO3 = xm_N + xm_O * 3 ! kg/mol 196 REAL, PARAMETER :: xm_CO2 = xm_C + xm_O * 2 ! kg/mol 197 198 ! mass of air 199 REAL, PARAMETER :: xm_air = 28.964e-3 ! kg/mol 200 201 ! dummy weight, used for complex molecules: 202 REAL, PARAMETER :: xm_dummy = 1000.0e-3 ! kg/mol 203 204 163 205 SAVE 164 206 END MODULE chem_modules -
palm/trunk/SOURCE/chemistry_model_mod.f90
r3449 r3458 27 27 ! ----------------- 28 28 ! $Id$ 29 ! from chemistry branch r3443, banzhafs, basit: 30 ! replace surf_lsm_h%qv1(m) by q(k,j,i) for mixing ratio in chem_depo 31 ! bug fix in chem_depo: allow different surface fractions for one 32 ! surface element and set lai to zero for non vegetated surfaces 33 ! bug fixed in chem_data_output_2d 34 ! bug fix in chem_depo subroutine 35 ! added code on deposition of gases and particles 36 ! removed cs_profile_name from chem_parin 37 ! bug fixed in output profiles and code cleaned 38 ! 39 ! 3449 2018-10-29 19:36:56Z suehring 29 40 ! additional output - merged from branch resler 30 41 ! … … 151 162 ! @author Klaus Ketelsen 152 163 ! @author Basit Khan 164 ! @author Sabine Banzhaf 153 165 ! 154 166 ! … … 197 209 initializing_actions, message_string, & 198 210 omega, tsc, intermediate_timestep_count, & 199 intermediate_timestep_count_max, 211 intermediate_timestep_count_max, max_pr_user, & 200 212 timestep_scheme, use_prescribed_profile_data, ws_scheme_sca 201 213 USE arrays_3d, ONLY: exner, hyp, pt, q, ql, rdf_sc, tend, zu … … 271 283 ! 'fastj' (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED 272 284 285 286 !----------------------------------------------------------------------- 287 ! Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010) 288 ! 289 INTEGER(iwp), PARAMETER :: nlu_dep = 15 ! Number of DEPAC landuse classes (lu's) 290 INTEGER(iwp), PARAMETER :: ncmp = 10 ! Number of DEPAC gas components 291 !------------------------------------------------------------------------ 292 ! DEPAC landuse classes as defined in LOTOS-EUROS model v2.1 293 ! 294 INTEGER(iwp) :: ilu_grass = 1 295 INTEGER(iwp) :: ilu_arable = 2 296 INTEGER(iwp) :: ilu_permanent_crops = 3 297 INTEGER(iwp) :: ilu_coniferous_forest = 4 298 INTEGER(iwp) :: ilu_deciduous_forest = 5 299 INTEGER(iwp) :: ilu_water_sea = 6 300 INTEGER(iwp) :: ilu_urban = 7 301 INTEGER(iwp) :: ilu_other = 8 302 INTEGER(iwp) :: ilu_desert = 9 303 INTEGER(iwp) :: ilu_ice = 10 304 INTEGER(iwp) :: ilu_savanna = 11 305 INTEGER(iwp) :: ilu_tropical_forest = 12 306 INTEGER(iwp) :: ilu_water_inland = 13 307 INTEGER(iwp) :: ilu_mediterrean_scrub = 14 308 INTEGER(iwp) :: ilu_semi_natural_veg = 15 309 ! 310 !-------------------------------------------------------------------------- 311 ! NH3/SO2 ratio regimes: 312 INTEGER, PARAMETER :: iratns_low = 1 ! low ratio NH3/SO2 313 INTEGER, PARAMETER :: iratns_high = 2 ! high ratio NH3/SO2 314 INTEGER, PARAMETER :: iratns_very_low = 3 ! very low ratio NH3/SO2 315 ! Default: 316 INTEGER, PARAMETER :: iratns_default = iratns_low 317 ! 318 ! Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) naar W m-2 319 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: alpha =(/0.009,0.009, 0.009,0.006,0.006, -999., -999.,0.009,-999.,-999.,0.009,0.006,-999.,0.009,0.008/)*4.57 320 ! 321 ! Set temperatures per land use for F_temp 322 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Tmin = (/12.0, 12.0, 12.0, 0.0, 0.0, -999., -999., 12.0, -999., -999., 12.0, 0.0, -999., 12.0, 8.0/) 323 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Topt = (/26.0, 26.0, 26.0, 18.0, 20.0, -999., -999., 26.0, -999., -999., 26.0, 20.0, -999., 26.0, 24.0 /) 324 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Tmax = (/40.0, 40.0, 40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999., 40.0, 35.0, -999., 40.0, 39.0 /) 325 ! 326 ! Set F_min: 327 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: F_min =(/0.01, 0.01, 0.01, 0.1, 0.1, -999., -999.,0.01, -999.,-999.,0.01,0.1,-999.,0.01, 0.04/) 328 329 ! Set maximal conductance (m/s) 330 ! (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from mmol O3/m2/s to m/s 331 ! Could be refined to a function of T and P. in Jones 332 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: g_max =(/270., 300., 300., 140., 150., -999., -999.,270., -999.,-999.,270., 150.,-999.,300., 422./)/41000 333 ! 334 ! Set max, min for vapour pressure deficit vpd; 335 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: vpd_max =(/1.3, 0.9, 0.9, 0.5, 1.0, -999., -999.,1.3, -999.,-999.,1.3,1.0, -999.,0.9, 2.8/) 336 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: vpd_min =(/3.0, 2.8, 2.8, 3.0, 3.25, -999., -999.,3.0, -999.,-999.,3.0,3.25, -999.,2.8, 4.5/) 337 ! 338 ! 339 !------------------------------------------------------------------------ 340 341 273 342 PUBLIC nest_chemistry 274 343 PUBLIC nspec … … 357 426 END INTERFACE chem_wrd_local 358 427 428 INTERFACE chem_depo 429 MODULE PROCEDURE chem_depo 430 END INTERFACE chem_depo 431 432 INTERFACE drydepos_gas_depac 433 MODULE PROCEDURE drydepos_gas_depac 434 END INTERFACE drydepos_gas_depac 435 436 INTERFACE rc_special 437 MODULE PROCEDURE rc_special 438 END INTERFACE rc_special 439 440 INTERFACE rc_gw 441 MODULE PROCEDURE rc_gw 442 END INTERFACE rc_gw 443 444 INTERFACE rw_so2 445 MODULE PROCEDURE rw_so2 446 END INTERFACE rw_so2 447 448 INTERFACE rw_nh3_sutton 449 MODULE PROCEDURE rw_nh3_sutton 450 END INTERFACE rw_nh3_sutton 451 452 INTERFACE rw_constant 453 MODULE PROCEDURE rw_constant 454 END INTERFACE rw_constant 455 456 INTERFACE rc_gstom 457 MODULE PROCEDURE rc_gstom 458 END INTERFACE rc_gstom 459 460 INTERFACE rc_gstom_emb 461 MODULE PROCEDURE rc_gstom_emb 462 END INTERFACE rc_gstom_emb 463 464 INTERFACE par_dir_diff 465 MODULE PROCEDURE par_dir_diff 466 END INTERFACE par_dir_diff 467 468 INTERFACE rc_get_vpd 469 MODULE PROCEDURE rc_get_vpd 470 END INTERFACE rc_get_vpd 471 472 INTERFACE rc_gsoil_eff 473 MODULE PROCEDURE rc_gsoil_eff 474 END INTERFACE rc_gsoil_eff 475 476 INTERFACE rc_rinc 477 MODULE PROCEDURE rc_rinc 478 END INTERFACE rc_rinc 479 480 INTERFACE rc_rctot 481 MODULE PROCEDURE rc_rctot 482 END INTERFACE rc_rctot 483 484 INTERFACE rc_comp_point_rc_eff 485 MODULE PROCEDURE rc_comp_point_rc_eff 486 END INTERFACE rc_comp_point_rc_eff 487 488 INTERFACE drydepo_aero_zhang_vd 489 MODULE PROCEDURE drydepo_aero_zhang_vd 490 END INTERFACE drydepo_aero_zhang_vd 491 492 INTERFACE get_rb_cell 493 MODULE PROCEDURE get_rb_cell 494 END INTERFACE get_rb_cell 495 496 359 497 360 498 PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, & … … 364 502 chem_init_profiles, chem_integrate, chem_parin, & 365 503 chem_prognostic_equations, chem_rrd_local, & 366 chem_statistics, chem_swap_timelevel, chem_wrd_local 504 chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo 367 505 368 506 … … 838 976 ! It is possible to plant PM10 and PM25 into the gasphase chemistry code 839 977 ! as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): 840 ! set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)978 ! set unit to kilograms per m**3 for PM10 and PM25 (PM2.5) 841 979 IF (spec_name(1:2) == 'PM') THEN 842 980 unit = 'kg m-3' … … 852 990 853 991 854 RETURN992 RETURN 855 993 END SUBROUTINE chem_check_data_output 856 994 ! … … 1072 1210 DO k = nzb_do, nzt_do 1073 1211 local_pf(i,j,k) = MERGE( & 1074 chem_species(lsp)%conc (k,j,i),&1212 chem_species(lsp)%conc_av(k,j,i), & 1075 1213 REAL( fill_value, KIND = wp ), & 1076 1214 BTEST( wall_flags_0(k,j,i), 0 ) ) … … 1104 1242 USE surface_mod 1105 1243 1106 !USE chem_modules, ONLY: nvar1107 1244 1108 1245 IMPLICIT NONE … … 1239 1376 1240 1377 spec_name = TRIM( variable(4:) ) 1241 !av == 01242 1378 1243 1379 DO lsp=1,nspec … … 1870 2006 cs_name, & 1871 2007 cs_profile, & 1872 cs_profile_name, &1873 2008 cs_surface, & 1874 2009 decycle_chem_lr, & 1875 2010 decycle_chem_ns, & 1876 decycle_method, & 2011 decycle_method, & 2012 do_depo, & 1877 2013 emiss_factor_main, & 1878 2014 emiss_factor_side, & … … 2365 2501 DO lpr = 1, cs_pr_count 2366 2502 2367 sums_l(k,pr_palm+ lpr,tn) = sums_l(k,pr_palm+lpr,tn) + &2503 sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) + & 2368 2504 chem_species(cs_pr_index(lpr))%conc(k,j,i) * & 2369 2505 rmask(j,i,sr) * & … … 2432 2568 2433 2569 END SUBROUTINE chem_wrd_local 2570 2571 2572 2573 !-------------------------------------------------------------------------------! 2574 ! 2575 ! Description: 2576 ! ------------ 2577 !> Subroutine to calculate the deposition of gases and PMs. For now deposition 2578 !> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT 2579 !> considered. The deposition of particles is derived following Zhang et al., 2580 !> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010). 2581 !> 2582 !> @TODO: Consider deposition on vertical surfaces 2583 !> @TODO: Consider overlaying horizontal surfaces 2584 !> @TODO: Check error messages 2585 !-------------------------------------------------------------------------------! 2586 2587 SUBROUTINE chem_depo (i,j) 2588 2589 2590 USE control_parameters, & 2591 ONLY: dt_3d, intermediate_timestep_count, latitude 2592 2593 USE arrays_3d, & 2594 ONLY: dzw, rho_air_zw 2595 2596 USE date_and_time_mod, & 2597 ONLY: day_of_year 2598 2599 USE surface_mod, & 2600 ONLY: surf_lsm_h, surf_usm_h, surf_type, ind_veg_wall, & 2601 ind_pav_green, ind_wat_win 2602 2603 USE radiation_model_mod, & 2604 ONLY: zenith 2605 2606 2607 2608 IMPLICIT NONE 2609 2610 INTEGER,INTENT(IN) :: i,j 2611 2612 INTEGER(iwp) :: k ! matching k to surface m at i,j 2613 INTEGER(iwp) :: lsp ! running index for chem spcs. 2614 INTEGER(iwp) :: lu ! running index for landuse classes 2615 INTEGER(iwp) :: lu_palm ! index of PALM LSM vegetation_type at current surface element 2616 INTEGER(iwp) :: lup_palm ! index of PALM LSM pavement_type at current surface element 2617 INTEGER(iwp) :: luw_palm ! index of PALM LSM water_type at current surface element 2618 INTEGER(iwp) :: luu_palm ! index of PALM USM walls/roofs at current surface element 2619 INTEGER(iwp) :: lug_palm ! index of PALM USM green walls/roofs at current surface element 2620 INTEGER(iwp) :: lud_palm ! index of PALM USM windows at current surface element 2621 INTEGER(iwp) :: lu_dep ! matching DEPAC LU to lu_palm 2622 INTEGER(iwp) :: lup_dep ! matching DEPAC LU to lup_palm 2623 INTEGER(iwp) :: luw_dep ! matching DEPAC LU to luw_palm 2624 INTEGER(iwp) :: luu_dep ! matching DEPAC LU to luu_palm 2625 INTEGER(iwp) :: lug_dep ! matching DEPAC LU to lug_palm 2626 INTEGER(iwp) :: lud_dep ! matching DEPAC LU to lud_palm 2627 INTEGER(iwp) :: m ! index for horizontal surfaces 2628 2629 INTEGER(iwp) :: pspec ! running index 2630 INTEGER(iwp) :: i_pspec ! index for matching depac gas component 2631 2632 2633 ! Vegetation !Match to DEPAC 2634 INTEGER(iwp) :: ind_lu_user = 0 ! --> ERROR 2635 INTEGER(iwp) :: ind_lu_b_soil = 1 ! --> ilu_desert 2636 INTEGER(iwp) :: ind_lu_mixed_crops = 2 ! --> ilu_arable 2637 INTEGER(iwp) :: ind_lu_s_grass = 3 ! --> ilu_grass 2638 INTEGER(iwp) :: ind_lu_ev_needle_trees = 4 ! --> ilu_coniferous_forest 2639 INTEGER(iwp) :: ind_lu_de_needle_trees = 5 ! --> ilu_coniferous_forest 2640 INTEGER(iwp) :: ind_lu_ev_broad_trees = 6 ! --> ilu_tropical_forest 2641 INTEGER(iwp) :: ind_lu_de_broad_trees = 7 ! --> ilu_deciduous_forest 2642 INTEGER(iwp) :: ind_lu_t_grass = 8 ! --> ilu_grass 2643 INTEGER(iwp) :: ind_lu_desert = 9 ! --> ilu_desert 2644 INTEGER(iwp) :: ind_lu_tundra = 10 ! --> ilu_other 2645 INTEGER(iwp) :: ind_lu_irr_crops = 11 ! --> ilu_arable 2646 INTEGER(iwp) :: ind_lu_semidesert = 12 ! --> ilu_other 2647 INTEGER(iwp) :: ind_lu_ice = 13 ! --> ilu_ice 2648 INTEGER(iwp) :: ind_lu_marsh = 14 ! --> ilu_other 2649 INTEGER(iwp) :: ind_lu_ev_shrubs = 15 ! --> ilu_mediterrean_scrub 2650 INTEGER(iwp) :: ind_lu_de_shrubs = 16 ! --> ilu_mediterrean_scrub 2651 INTEGER(iwp) :: ind_lu_mixed_forest = 17 ! --> ilu_coniferous_forest (ave(decid+conif)) 2652 INTEGER(iwp) :: ind_lu_intrup_forest = 18 ! --> ilu_other (ave(other+decid)) 2653 2654 ! Water 2655 INTEGER(iwp) :: ind_luw_user = 0 ! --> ERROR 2656 INTEGER(iwp) :: ind_luw_lake = 1 ! --> ilu_water_inland 2657 INTEGER(iwp) :: ind_luw_river = 2 ! --> ilu_water_inland 2658 INTEGER(iwp) :: ind_luw_ocean = 3 ! --> ilu_water_sea 2659 INTEGER(iwp) :: ind_luw_pond = 4 ! --> ilu_water_inland 2660 INTEGER(iwp) :: ind_luw_fountain = 5 ! --> ilu_water_inland 2661 2662 ! Pavement 2663 INTEGER(iwp) :: ind_lup_user = 0 ! --> ERROR 2664 INTEGER(iwp) :: ind_lup_asph_conc = 1 ! --> ilu_desert 2665 INTEGER(iwp) :: ind_lup_asph = 2 ! --> ilu_desert 2666 INTEGER(iwp) :: ind_lup_conc = 3 ! --> ilu_desert 2667 INTEGER(iwp) :: ind_lup_sett = 4 ! --> ilu_desert 2668 INTEGER(iwp) :: ind_lup_pav_stones = 5 ! --> ilu_desert 2669 INTEGER(iwp) :: ind_lup_cobblest = 6 ! --> ilu_desert 2670 INTEGER(iwp) :: ind_lup_metal = 7 ! --> ilu_desert 2671 INTEGER(iwp) :: ind_lup_wood = 8 ! --> ilu_desert 2672 INTEGER(iwp) :: ind_lup_gravel = 9 ! --> ilu_desert 2673 INTEGER(iwp) :: ind_lup_f_gravel = 10 ! --> ilu_desert 2674 INTEGER(iwp) :: ind_lup_pebblest = 11 ! --> ilu_desert 2675 INTEGER(iwp) :: ind_lup_woodchips = 12 ! --> ilu_desert 2676 INTEGER(iwp) :: ind_lup_tartan = 13 ! --> ilu_desert 2677 INTEGER(iwp) :: ind_lup_art_turf = 14 ! --> ilu_desert 2678 INTEGER(iwp) :: ind_lup_clay = 15 ! --> ilu_desert 2679 2680 2681 2682 !-- Particle parameters according to the respective aerosol classes (PM25, PM10) 2683 INTEGER(iwp) :: ind_p_size = 1 !< index for partsize in particle_pars 2684 INTEGER(iwp) :: ind_p_dens = 2 !< index for rhopart in particle_pars 2685 INTEGER(iwp) :: ind_p_slip = 3 !< index for slipcor in particle_pars 2686 2687 INTEGER(iwp) :: part_type 2688 2689 INTEGER(iwp) :: nwet ! wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow 2690 2691 REAL(kind=wp) :: dt_chem 2692 REAL(kind=wp) :: dh 2693 REAL(kind=wp) :: inv_dh 2694 REAL(kind=wp) :: dt_dh 2695 2696 REAL(kind=wp) :: dens ! Density at lowest layer at i,j 2697 REAL(kind=wp) :: r_aero_lu ! aerodynamic resistance (s/m) at current surface element 2698 REAL(kind=wp) :: ustar_lu ! ustar at current surface element 2699 REAL(kind=wp) :: z0h_lu ! roughness length for heat at current surface element 2700 REAL(kind=wp) :: glrad ! rad_sw_in at current surface element 2701 REAL(kind=wp) :: ppm_to_ugm3 ! conversion factor 2702 REAL(kind=wp) :: relh ! relative humidity at current surface element 2703 REAL(kind=wp) :: lai ! leaf area index at current surface element 2704 REAL(kind=wp) :: sai ! surface area index at current surface element assumed to be lai + 1 2705 2706 REAL(kind=wp) :: slinnfac 2707 REAL(kind=wp) :: visc ! Viscosity 2708 REAL(kind=wp) :: vs ! Sedimentation velocity 2709 REAL(kind=wp) :: vd_lu ! deposition velocity (m/s) 2710 REAL(kind=wp) :: Rs ! Sedimentaion resistance (s/m) 2711 REAL(kind=wp) :: Rb ! quasi-laminar boundary layer resistance (s/m) 2712 REAL(kind=wp) :: Rc_tot ! total canopy resistance Rc (s/m) 2713 2714 REAL(kind=wp) :: catm ! gasconc. at i,j,k=1 in ug/m3 2715 REAL(kind=wp) :: diffc ! Diffusivity 2716 2717 2718 REAL(kind=wp),DIMENSION(nspec) :: bud_now_lu ! budget for LSM vegetation type at current surface element 2719 REAL(kind=wp),DIMENSION(nspec) :: bud_now_lup ! budget for LSM pavement type at current surface element 2720 REAL(kind=wp),DIMENSION(nspec) :: bud_now_luw ! budget for LSM water type at current surface element 2721 REAL(kind=wp),DIMENSION(nspec) :: bud_now_luu ! budget for USM walls/roofs at current surface element 2722 REAL(kind=wp),DIMENSION(nspec) :: bud_now_lug ! budget for USM green surfaces at current surface element 2723 REAL(kind=wp),DIMENSION(nspec) :: bud_now_lud ! budget for USM windows at current surface element 2724 REAL(kind=wp),DIMENSION(nspec) :: bud_now ! overall budget at current surface element 2725 REAL(kind=wp),DIMENSION(nspec) :: cc ! concentration i,j,k 2726 REAL(kind=wp),DIMENSION(nspec) :: ccomp_tot ! total compensation point (ug/m3) (= 0 for species that don't have a compensation) 2727 ! For now kept to zero for all species! 2728 2729 REAL(kind=wp) :: ttemp ! temperatur at i,j,k 2730 REAL(kind=wp) :: ts ! surface temperatur in degrees celsius 2731 REAL(kind=wp) :: qv_tmp ! surface mixing ratio at current surface element 2732 2733 2734 !-- Particle parameters (PM10 (1), PM25 (2)) 2735 !-- partsize (diameter in m) , rhopart (density in kg/m3), slipcor (slip correction factor DIMENSIONless, Seinfeld and Pandis 2006, Table 9.3) 2736 REAL(wp), DIMENSION(1:3,1:2), PARAMETER :: particle_pars = RESHAPE( (/ & 2737 8.0e-6_wp, 1.14e3_wp, 1.016_wp, & ! 1 2738 0.7e-6_wp, 1.14e3_wp, 1.082_wp & ! 2 2739 /), (/ 3, 2 /) ) 2740 2741 2742 LOGICAL :: match_lsm ! flag indicating natural-type surface 2743 LOGICAL :: match_usm ! flag indicating urban-type surface 2744 2745 2746 2747 !-- List of names of possible tracers 2748 CHARACTER(len=*), PARAMETER :: pspecnames(69) = (/ & 2749 'NO2 ', & ! NO2 2750 'NO ', & ! NO 2751 'O3 ', & ! O3 2752 'CO ', & ! CO 2753 'form ', & ! FORM 2754 'ald ', & ! ALD 2755 'pan ', & ! PAN 2756 'mgly ', & ! MGLY 2757 'par ', & ! PAR 2758 'ole ', & ! OLE 2759 'eth ', & ! ETH 2760 'tol ', & ! TOL 2761 'cres ', & ! CRES 2762 'xyl ', & ! XYL 2763 'SO4a_f ', & ! SO4a_f 2764 'SO2 ', & ! SO2 2765 'HNO2 ', & ! HNO2 2766 'CH4 ', & ! CH4 2767 'NH3 ', & ! NH3 2768 'NO3 ', & ! NO3 2769 'OH ', & ! OH 2770 'HO2 ', & ! HO2 2771 'N2O5 ', & ! N2O5 2772 'SO4a_c ', & ! SO4a_c 2773 'NH4a_f ', & ! NH4a_f 2774 'NO3a_f ', & ! NO3a_f 2775 'NO3a_c ', & ! NO3a_c 2776 'C2O3 ', & ! C2O3 2777 'XO2 ', & ! XO2 2778 'XO2N ', & ! XO2N 2779 'cro ', & ! CRO 2780 'HNO3 ', & ! HNO3 2781 'H2O2 ', & ! H2O2 2782 'iso ', & ! ISO 2783 'ispd ', & ! ISPD 2784 'to2 ', & ! TO2 2785 'open ', & ! OPEN 2786 'terp ', & ! TERP 2787 'ec_f ', & ! EC_f 2788 'ec_c ', & ! EC_c 2789 'pom_f ', & ! POM_f 2790 'pom_c ', & ! POM_c 2791 'ppm_f ', & ! PPM_f 2792 'ppm_c ', & ! PPM_c 2793 'na_ff ', & ! Na_ff 2794 'na_f ', & ! Na_f 2795 'na_c ', & ! Na_c 2796 'na_cc ', & ! Na_cc 2797 'na_ccc ', & ! Na_ccc 2798 'dust_ff ', & ! dust_ff 2799 'dust_f ', & ! dust_f 2800 'dust_c ', & ! dust_c 2801 'dust_cc ', & ! dust_cc 2802 'dust_ccc ', & ! dust_ccc 2803 'tpm10 ', & ! tpm10 2804 'tpm25 ', & ! tpm25 2805 'tss ', & ! tss 2806 'tdust ', & ! tdust 2807 'tc ', & ! tc 2808 'tcg ', & ! tcg 2809 'tsoa ', & ! tsoa 2810 'tnmvoc ', & ! tnmvoc 2811 'SOxa ', & ! SOxa 2812 'NOya ', & ! NOya 2813 'NHxa ', & ! NHxa 2814 'NO2_obs ', & ! NO2_obs 2815 'tpm10_biascorr', & ! tpm10_biascorr 2816 'tpm25_biascorr', & ! tpm25_biascorr 2817 'O3_biascorr ' /) ! o3_biascorr 2818 2819 2820 2821 ! tracer mole mass: 2822 REAL, PARAMETER :: specmolm(69) = (/ & 2823 xm_O * 2 + xm_N, & ! NO2 2824 xm_O + xm_N, & ! NO 2825 xm_O * 3, & ! O3 2826 xm_C + xm_O, & ! CO 2827 xm_H * 2 + xm_C + xm_O, & ! FORM 2828 xm_H * 3 + xm_C * 2 + xm_O, & ! ALD 2829 xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, & ! PAN 2830 xm_H * 4 + xm_C * 3 + xm_O * 2, & ! MGLY 2831 xm_H * 3 + xm_C, & ! PAR 2832 xm_H * 3 + xm_C * 2, & ! OLE 2833 xm_H * 4 + xm_C * 2, & ! ETH 2834 xm_H * 8 + xm_C * 7, & ! TOL 2835 xm_H * 8 + xm_C * 7 + xm_O, & ! CRES 2836 xm_H * 10 + xm_C * 8, & ! XYL 2837 xm_S + xm_O * 4, & ! SO4a_f 2838 xm_S + xm_O * 2, & ! SO2 2839 xm_H + xm_O * 2 + xm_N, & ! HNO2 2840 xm_H * 4 + xm_C, & ! CH4 2841 xm_H * 3 + xm_N, & ! NH3 2842 xm_O * 3 + xm_N, & ! NO3 2843 xm_H + xm_O, & ! OH 2844 xm_H + xm_O * 2, & ! HO2 2845 xm_O * 5 + xm_N * 2, & ! N2O5 2846 xm_S + xm_O * 4, & ! SO4a_c 2847 xm_H * 4 + xm_N, & ! NH4a_f 2848 xm_O * 3 + xm_N, & ! NO3a_f 2849 xm_O * 3 + xm_N, & ! NO3a_c 2850 xm_C * 2 + xm_O * 3, & ! C2O3 2851 xm_dummy, & ! XO2 2852 xm_dummy, & ! XO2N 2853 xm_dummy, & ! CRO 2854 xm_H + xm_O * 3 + xm_N, & ! HNO3 2855 xm_H * 2 + xm_O * 2, & ! H2O2 2856 xm_H * 8 + xm_C * 5, & ! ISO 2857 xm_dummy, & ! ISPD 2858 xm_dummy, & ! TO2 2859 xm_dummy, & ! OPEN 2860 xm_H * 16 + xm_C * 10, & ! TERP 2861 xm_dummy, & ! EC_f 2862 xm_dummy, & ! EC_c 2863 xm_dummy, & ! POM_f 2864 xm_dummy, & ! POM_c 2865 xm_dummy, & ! PPM_f 2866 xm_dummy, & ! PPM_c 2867 xm_Na, & ! Na_ff 2868 xm_Na, & ! Na_f 2869 xm_Na, & ! Na_c 2870 xm_Na, & ! Na_cc 2871 xm_Na, & ! Na_ccc 2872 xm_dummy, & ! dust_ff 2873 xm_dummy, & ! dust_f 2874 xm_dummy, & ! dust_c 2875 xm_dummy, & ! dust_cc 2876 xm_dummy, & ! dust_ccc 2877 xm_dummy, & ! tpm10 2878 xm_dummy, & ! tpm25 2879 xm_dummy, & ! tss 2880 xm_dummy, & ! tdust 2881 xm_dummy, & ! tc 2882 xm_dummy, & ! tcg 2883 xm_dummy, & ! tsoa 2884 xm_dummy, & ! tnmvoc 2885 xm_dummy, & ! SOxa 2886 xm_dummy, & ! NOya 2887 xm_dummy, & ! NHxa 2888 xm_O * 2 + xm_N, & ! NO2_obs 2889 xm_dummy, & ! tpm10_biascorr 2890 xm_dummy, & ! tpm25_biascorr 2891 xm_O * 3 /) ! o3_biascorr 2892 2893 2894 2895 ! Initialize surface element m 2896 m = 0 2897 k = 0 2898 ! LSM or USM surface present at i,j: 2899 ! Default surfaces are NOT considered for deposition 2900 match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i) 2901 match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i) 2902 2903 2904 !!!!!!! 2905 ! LSM ! 2906 !!!!!!! 2907 2908 IF ( match_lsm ) THEN 2909 2910 ! Get surface element information at i,j: 2911 m = surf_lsm_h%start_index(j,i) 2912 2913 ! Get corresponding k 2914 k = surf_lsm_h%k(m) 2915 2916 ! Get needed variables for surface element m 2917 ustar_lu = surf_lsm_h%us(m) 2918 z0h_lu = surf_lsm_h%z0h(m) 2919 r_aero_lu = surf_lsm_h%r_a(m) 2920 glrad = surf_lsm_h%rad_sw_in(m) 2921 lai = surf_lsm_h%lai(m) 2922 sai = lai + 1 2923 2924 ! For small grid spacing neglect R_a 2925 IF (dzw(k) <= 1.0) THEN 2926 r_aero_lu = 0.0_wp 2927 ENDIF 2928 2929 !Initialize lu's 2930 lu_palm = 0 2931 lu_dep = 0 2932 lup_palm = 0 2933 lup_dep = 0 2934 luw_palm = 0 2935 luw_dep = 0 2936 2937 !Initialize budgets 2938 bud_now_lu = 0.0_wp 2939 bud_now_lup = 0.0_wp 2940 bud_now_luw = 0.0_wp 2941 2942 2943 ! Get land use for i,j and assign to DEPAC lu 2944 IF (surf_lsm_h%frac(ind_veg_wall,m) > 0) THEN 2945 lu_palm = surf_lsm_h%vegetation_type(m) 2946 IF (lu_palm == ind_lu_user) THEN 2947 message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation' 2948 CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 ) 2949 ELSEIF (lu_palm == ind_lu_b_soil) THEN 2950 lu_dep = 9 2951 ELSEIF (lu_palm == ind_lu_mixed_crops) THEN 2952 lu_dep = 2 2953 ELSEIF (lu_palm == ind_lu_s_grass) THEN 2954 lu_dep = 1 2955 ELSEIF (lu_palm == ind_lu_ev_needle_trees) THEN 2956 lu_dep = 4 2957 ELSEIF (lu_palm == ind_lu_de_needle_trees) THEN 2958 lu_dep = 4 2959 ELSEIF (lu_palm == ind_lu_ev_broad_trees) THEN 2960 lu_dep = 12 2961 ELSEIF (lu_palm == ind_lu_de_broad_trees) THEN 2962 lu_dep = 5 2963 ELSEIF (lu_palm == ind_lu_t_grass) THEN 2964 lu_dep = 1 2965 ELSEIF (lu_palm == ind_lu_desert) THEN 2966 lu_dep = 9 2967 ELSEIF (lu_palm == ind_lu_tundra ) THEN 2968 lu_dep = 8 2969 ELSEIF (lu_palm == ind_lu_irr_crops) THEN 2970 lu_dep = 2 2971 ELSEIF (lu_palm == ind_lu_semidesert) THEN 2972 lu_dep = 8 2973 ELSEIF (lu_palm == ind_lu_ice) THEN 2974 lu_dep = 10 2975 ELSEIF (lu_palm == ind_lu_marsh) THEN 2976 lu_dep = 8 2977 ELSEIF (lu_palm == ind_lu_ev_shrubs) THEN 2978 lu_dep = 14 2979 ELSEIF (lu_palm == ind_lu_de_shrubs ) THEN 2980 lu_dep = 14 2981 ELSEIF (lu_palm == ind_lu_mixed_forest) THEN 2982 lu_dep = 4 2983 ELSEIF (lu_palm == ind_lu_intrup_forest) THEN 2984 lu_dep = 8 2985 ENDIF 2986 ENDIF 2987 2988 IF (surf_lsm_h%frac(ind_pav_green,m) > 0) THEN 2989 lup_palm = surf_lsm_h%pavement_type(m) 2990 IF (lup_palm == ind_lup_user) THEN 2991 message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation' 2992 CALL message( 'chem_depo', 'CM0452', 1, 2, 0, 6, 0 ) 2993 ELSEIF (lup_palm == ind_lup_asph_conc) THEN 2994 lup_dep = 9 2995 ELSEIF (lup_palm == ind_lup_asph) THEN 2996 lup_dep = 9 2997 ELSEIF (lup_palm == ind_lup_conc) THEN 2998 lup_dep = 9 2999 ELSEIF (lup_palm == ind_lup_sett) THEN 3000 lup_dep = 9 3001 ELSEIF (lup_palm == ind_lup_pav_stones) THEN 3002 lup_dep = 9 3003 ELSEIF (lup_palm == ind_lup_cobblest) THEN 3004 lup_dep = 9 3005 ELSEIF (lup_palm == ind_lup_metal) THEN 3006 lup_dep = 9 3007 ELSEIF (lup_palm == ind_lup_wood) THEN 3008 lup_dep = 9 3009 ELSEIF (lup_palm == ind_lup_gravel) THEN 3010 lup_dep = 9 3011 ELSEIF (lup_palm == ind_lup_f_gravel) THEN 3012 lup_dep = 9 3013 ELSEIF (lup_palm == ind_lup_pebblest) THEN 3014 lup_dep = 9 3015 ELSEIF (lup_palm == ind_lup_woodchips) THEN 3016 lup_dep = 9 3017 ELSEIF (lup_palm == ind_lup_tartan) THEN 3018 lup_dep = 9 3019 ELSEIF (lup_palm == ind_lup_art_turf) THEN 3020 lup_dep = 9 3021 ELSEIF (lup_palm == ind_lup_clay) THEN 3022 lup_dep = 9 3023 ENDIF 3024 ENDIF 3025 3026 IF (surf_lsm_h%frac(ind_wat_win,m) > 0) THEN 3027 luw_palm = surf_lsm_h%water_type(m) 3028 IF (luw_palm == ind_luw_user) THEN 3029 message_string = 'No water type defined. Please define water type to enable deposition calculation' 3030 CALL message( 'chem_depo', 'CM0453', 1, 2, 0, 6, 0 ) 3031 ELSEIF (luw_palm == ind_luw_lake) THEN 3032 luw_dep = 13 3033 ELSEIF (luw_palm == ind_luw_river) THEN 3034 luw_dep = 13 3035 ELSEIF (luw_palm == ind_luw_ocean) THEN 3036 luw_dep = 6 3037 ELSEIF (luw_palm == ind_luw_pond) THEN 3038 luw_dep = 13 3039 ELSEIF (luw_palm == ind_luw_fountain) THEN 3040 luw_dep = 13 3041 ENDIF 3042 ENDIF 3043 3044 3045 3046 ! Set wetness indicator to dry or wet for lsm vegetation or pavement 3047 IF (surf_lsm_h%c_liq(m) > 0) THEN 3048 nwet = 1 3049 ELSE 3050 nwet = 0 3051 ENDIF 3052 3053 ! Compute length of time step 3054 IF ( call_chem_at_all_substeps ) THEN 3055 dt_chem = dt_3d * weight_pres(intermediate_timestep_count) 3056 ELSE 3057 dt_chem = dt_3d 3058 ENDIF 3059 3060 3061 dh = dzw(k) 3062 inv_dh = 1.0_wp / dh 3063 dt_dh = dt_chem/dh 3064 3065 ! Concentration at i,j,k 3066 DO lsp = 1, nspec 3067 cc(lsp) = chem_species(lsp)%conc(k,j,i) 3068 ENDDO 3069 3070 3071 ! Temperature column at i,j 3072 ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp 3073 3074 ! Surface temperature in degrees Celcius: 3075 ts = ttemp - 273.15 ! in degrees celcius 3076 3077 ! Viscosity of air in lowest layer 3078 visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0) 3079 3080 ! Air density in lowest layer 3081 dens = rho_air_zw(k) 3082 3083 ! Calculate relative humidity from specific humidity for DEPAC 3084 qv_tmp = max(q(k,j,i),0.0_wp) 3085 relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(k) ) 3086 3087 3088 3089 ! Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget 3090 ! for each surface fraction. Then derive overall budget taking into account the surface fractions. 3091 3092 ! Vegetation 3093 IF (surf_lsm_h%frac(ind_veg_wall,m) > 0) THEN 3094 3095 3096 slinnfac = 1.0_wp 3097 3098 ! Get vd 3099 DO lsp = 1, nvar 3100 !Initialize 3101 vs = 0.0_wp 3102 vd_lu = 0.0_wp 3103 Rs = 0.0_wp 3104 Rb = 0.0_wp 3105 Rc_tot = 0.0_wp 3106 IF (spc_names(lsp) == 'PM10') THEN 3107 part_type = 1 3108 ! sedimentation velocity in lowest layer 3109 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3110 particle_pars(ind_p_size, part_type), & 3111 particle_pars(ind_p_slip, part_type), & 3112 visc) 3113 3114 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 3115 vs, & 3116 particle_pars(ind_p_size, part_type), & 3117 particle_pars(ind_p_slip, part_type), & 3118 nwet, ttemp, dens, visc, & 3119 lu_dep, & 3120 r_aero_lu, ustar_lu) 3121 3122 bud_now_lu(lsp) = - cc(lsp) * & 3123 (1.0 - exp(-vd_lu*dt_dh ))*dh 3124 3125 3126 ELSEIF (spc_names(lsp) == 'PM25') THEN 3127 part_type = 2 3128 ! sedimentation velocity in lowest layer: 3129 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3130 particle_pars(ind_p_size, part_type), & 3131 particle_pars(ind_p_slip, part_type), & 3132 visc) 3133 3134 3135 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 3136 vs, & 3137 particle_pars(ind_p_size, part_type), & 3138 particle_pars(ind_p_slip, part_type), & 3139 nwet, ttemp, dens, visc, & 3140 lu_dep , & 3141 r_aero_lu, ustar_lu) 3142 3143 bud_now_lu(lsp) = - cc(lsp) * & 3144 (1.0 - exp(-vd_lu*dt_dh ))*dh 3145 3146 3147 ELSE 3148 3149 ! GASES 3150 3151 ! Read spc_name of current species for gas parameter 3152 3153 IF (ANY(pspecnames(:) == spc_names(lsp))) THEN 3154 i_pspec = 0 3155 DO pspec = 1, 69 3156 IF (pspecnames(pspec) == spc_names(lsp)) THEN 3157 i_pspec = pspec 3158 END IF 3159 ENDDO 3160 3161 ELSE 3162 ! Species for now not deposited 3163 CYCLE 3164 ENDIF 3165 3166 ! Factor used for conversion from ppb to ug/m3 : 3167 ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 3168 ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 3169 ! c 1e-9 xm_tracer 1e9 / xm_air dens 3170 ! thus: 3171 ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 3172 ! Use density at lowest layer: 3173 3174 ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 3175 3176 ! atmospheric concentration in DEPAC is requested in ug/m3: 3177 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3178 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 3179 3180 ! Diffusivity for DEPAC relevant gases 3181 ! Use default value 3182 diffc = 0.11e-4 3183 ! overwrite with known coefficients of diffusivity from Massman (1998) 3184 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 3185 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 3186 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 3187 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 3188 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 3189 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 3190 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 3191 3192 3193 ! Get quasi-laminar boundary layer resistance Rb: 3194 CALL get_rb_cell( (lu_dep == ilu_water_sea) .or. (lu_dep == ilu_water_inland), & 3195 z0h_lu, ustar_lu, diffc, & 3196 Rb) 3197 3198 ! Get Rc_tot 3199 CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & 3200 relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3201 r_aero_lu , Rb) 3202 3203 3204 ! Calculate budget 3205 IF (Rc_tot .le. 0.0) THEN 3206 3207 bud_now_lu(lsp) = 0.0_wp 3208 3209 ELSE 3210 3211 ! Compute exchange velocity for current lu: 3212 vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) 3213 3214 bud_now_lu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & 3215 (1.0 - exp(-vd_lu*dt_dh ))*dh 3216 ENDIF 3217 3218 ENDIF 3219 ENDDO 3220 ENDIF 3221 3222 3223 ! Pavement 3224 IF (surf_lsm_h%frac(ind_pav_green,m) > 0) THEN 3225 3226 3227 ! No vegetation on pavements: 3228 lai = 0.0_wp 3229 sai = 0.0_wp 3230 3231 slinnfac = 1.0_wp 3232 3233 ! Get vd 3234 DO lsp = 1, nvar 3235 !Initialize 3236 vs = 0.0_wp 3237 vd_lu = 0.0_wp 3238 Rs = 0.0_wp 3239 Rb = 0.0_wp 3240 Rc_tot = 0.0_wp 3241 IF (spc_names(lsp) == 'PM10') THEN 3242 part_type = 1 3243 ! sedimentation velocity in lowest layer: 3244 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3245 particle_pars(ind_p_size, part_type), & 3246 particle_pars(ind_p_slip, part_type), & 3247 visc) 3248 3249 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 3250 vs, & 3251 particle_pars(ind_p_size, part_type), & 3252 particle_pars(ind_p_slip, part_type), & 3253 nwet, ttemp, dens, visc, & 3254 lup_dep, & 3255 r_aero_lu, ustar_lu) 3256 3257 bud_now_lup(lsp) = - cc(lsp) * & 3258 (1.0 - exp(-vd_lu*dt_dh ))*dh 3259 3260 3261 ELSEIF (spc_names(lsp) == 'PM25') THEN 3262 part_type = 2 3263 ! sedimentation velocity in lowest layer: 3264 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3265 particle_pars(ind_p_size, part_type), & 3266 particle_pars(ind_p_slip, part_type), & 3267 visc) 3268 3269 3270 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 3271 vs, & 3272 particle_pars(ind_p_size, part_type), & 3273 particle_pars(ind_p_slip, part_type), & 3274 nwet, ttemp, dens, visc, & 3275 lup_dep, & 3276 r_aero_lu, ustar_lu) 3277 3278 bud_now_lup(lsp) = - cc(lsp) * & 3279 (1.0 - exp(-vd_lu*dt_dh ))*dh 3280 3281 3282 ELSE 3283 3284 ! GASES 3285 3286 ! Read spc_name of current species for gas parameter 3287 3288 IF (ANY(pspecnames(:) == spc_names(lsp))) THEN 3289 i_pspec = 0 3290 DO pspec = 1, 69 3291 IF (pspecnames(pspec) == spc_names(lsp)) THEN 3292 i_pspec = pspec 3293 END IF 3294 ENDDO 3295 3296 ELSE 3297 ! Species for now not deposited 3298 CYCLE 3299 ENDIF 3300 3301 ! Factor used for conversion from ppb to ug/m3 : 3302 ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 3303 ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 3304 ! c 1e-9 xm_tracer 1e9 / xm_air dens 3305 ! thus: 3306 ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 3307 ! Use density at lowest layer: 3308 3309 ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 3310 3311 ! atmospheric concentration in DEPAC is requested in ug/m3: 3312 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3313 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 3314 3315 ! Diffusivity for DEPAC relevant gases 3316 ! Use default value 3317 diffc = 0.11e-4 3318 ! overwrite with known coefficients of diffusivity from Massman (1998) 3319 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 3320 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 3321 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 3322 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 3323 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 3324 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 3325 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 3326 3327 3328 ! Get quasi-laminar boundary layer resistance Rb: 3329 CALL get_rb_cell( (lup_dep == ilu_water_sea) .or. (lup_dep == ilu_water_inland), & 3330 z0h_lu, ustar_lu, diffc, & 3331 Rb) 3332 3333 3334 ! Get Rc_tot 3335 CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & 3336 relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3337 r_aero_lu , Rb) 3338 3339 3340 ! Calculate budget 3341 IF (Rc_tot .le. 0.0) THEN 3342 3343 bud_now_lup(lsp) = 0.0_wp 3344 3345 ELSE 3346 3347 ! Compute exchange velocity for current lu: 3348 vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) 3349 3350 bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & 3351 (1.0 - exp(-vd_lu*dt_dh ))*dh 3352 ENDIF 3353 3354 3355 ENDIF 3356 ENDDO 3357 ENDIF 3358 3359 3360 ! Water 3361 IF (surf_lsm_h%frac(ind_wat_win,m) > 0) THEN 3362 3363 3364 ! No vegetation on water: 3365 lai = 0.0_wp 3366 sai = 0.0_wp 3367 3368 slinnfac = 1.0_wp 3369 3370 ! Get vd 3371 DO lsp = 1, nvar 3372 !Initialize 3373 vs = 0.0_wp 3374 vd_lu = 0.0_wp 3375 Rs = 0.0_wp 3376 Rb = 0.0_wp 3377 Rc_tot = 0.0_wp 3378 IF (spc_names(lsp) == 'PM10') THEN 3379 part_type = 1 3380 ! sedimentation velocity in lowest layer: 3381 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3382 particle_pars(ind_p_size, part_type), & 3383 particle_pars(ind_p_slip, part_type), & 3384 visc) 3385 3386 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 3387 vs, & 3388 particle_pars(ind_p_size, part_type), & 3389 particle_pars(ind_p_slip, part_type), & 3390 nwet, ttemp, dens, visc, & 3391 luw_dep, & 3392 r_aero_lu, ustar_lu) 3393 3394 bud_now_luw(lsp) = - cc(lsp) * & 3395 (1.0 - exp(-vd_lu*dt_dh ))*dh 3396 3397 3398 ELSEIF (spc_names(lsp) == 'PM25') THEN 3399 part_type = 2 3400 ! sedimentation velocity in lowest layer: 3401 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3402 particle_pars(ind_p_size, part_type), & 3403 particle_pars(ind_p_slip, part_type), & 3404 visc) 3405 3406 3407 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 3408 vs, & 3409 particle_pars(ind_p_size, part_type), & 3410 particle_pars(ind_p_slip, part_type), & 3411 nwet, ttemp, dens, visc, & 3412 luw_dep, & 3413 r_aero_lu, ustar_lu) 3414 3415 bud_now_luw(lsp) = - cc(lsp) * & 3416 (1.0 - exp(-vd_lu*dt_dh ))*dh 3417 3418 3419 ELSE 3420 3421 ! GASES 3422 3423 ! Read spc_name of current species for gas PARAMETER 3424 3425 IF (ANY(pspecnames(:) == spc_names(lsp))) THEN 3426 i_pspec = 0 3427 DO pspec = 1, 69 3428 IF (pspecnames(pspec) == spc_names(lsp)) THEN 3429 i_pspec = pspec 3430 END IF 3431 ENDDO 3432 3433 ELSE 3434 ! Species for now not deposited 3435 CYCLE 3436 ENDIF 3437 3438 ! Factor used for conversion from ppb to ug/m3 : 3439 ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 3440 ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 3441 ! c 1e-9 xm_tracer 1e9 / xm_air dens 3442 ! thus: 3443 ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 3444 ! Use density at lowest layer: 3445 3446 ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 3447 3448 ! atmospheric concentration in DEPAC is requested in ug/m3: 3449 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3450 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 3451 3452 ! Diffusivity for DEPAC relevant gases 3453 ! Use default value 3454 diffc = 0.11e-4 3455 ! overwrite with known coefficients of diffusivity from Massman (1998) 3456 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 3457 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 3458 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 3459 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 3460 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 3461 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 3462 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 3463 3464 3465 ! Get quasi-laminar boundary layer resistance Rb: 3466 CALL get_rb_cell( (luw_dep == ilu_water_sea) .or. (luw_dep == ilu_water_inland), & 3467 z0h_lu, ustar_lu, diffc, & 3468 Rb) 3469 3470 ! Get Rc_tot 3471 CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & 3472 relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3473 r_aero_lu , Rb) 3474 3475 3476 ! Calculate budget 3477 IF (Rc_tot .le. 0.0) THEN 3478 3479 bud_now_luw(lsp) = 0.0_wp 3480 3481 ELSE 3482 3483 ! Compute exchange velocity for current lu: 3484 vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) 3485 3486 bud_now_luw(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & 3487 (1.0 - exp(-vd_lu*dt_dh ))*dh 3488 ENDIF 3489 3490 ENDIF 3491 ENDDO 3492 ENDIF 3493 3494 3495 bud_now = 0.0_wp 3496 ! Calculate budget for surface m and adapt concentration 3497 DO lsp = 1, nspec 3498 3499 3500 bud_now(lsp) = surf_lsm_h%frac(ind_veg_wall,m)*bud_now_lu(lsp) + & 3501 surf_lsm_h%frac(ind_pav_green,m)*bud_now_lup(lsp) + & 3502 surf_lsm_h%frac(ind_wat_win,m)*bud_now_luw(lsp) 3503 3504 ! Compute new concentration, add contribution of all exchange processes: 3505 cc(lsp) = cc(lsp) + bud_now(lsp)*inv_dh 3506 3507 chem_species(lsp)%conc(k,j,i) = max(0.0_wp, cc(lsp)) 3508 3509 ENDDO 3510 3511 ENDIF 3512 3513 3514 3515 !!!!!!! 3516 ! USM ! 3517 !!!!!!! 3518 3519 IF ( match_usm ) THEN 3520 3521 ! Get surface element information at i,j: 3522 m = surf_usm_h%start_index(j,i) 3523 3524 k = surf_usm_h%k(m) 3525 3526 ! Get needed variables for surface element m 3527 ustar_lu = surf_usm_h%us(m) 3528 z0h_lu = surf_usm_h%z0h(m) 3529 r_aero_lu = surf_usm_h%r_a(m) 3530 glrad = surf_usm_h%rad_sw_in(m) 3531 lai = surf_usm_h%lai(m) 3532 sai = lai + 1 3533 3534 ! For small grid spacing neglect R_a 3535 IF (dzw(k) <= 1.0) THEN 3536 r_aero_lu = 0.0_wp 3537 ENDIF 3538 3539 !Initialize lu's 3540 luu_palm = 0 3541 luu_dep = 0 3542 lug_palm = 0 3543 lug_dep = 0 3544 lud_palm = 0 3545 lud_dep = 0 3546 3547 !Initialize budgets 3548 bud_now_luu = 0.0_wp 3549 bud_now_lug = 0.0_wp 3550 bud_now_lud = 0.0_wp 3551 3552 3553 ! Get land use for i,j and assign to DEPAC lu 3554 IF (surf_usm_h%frac(ind_pav_green,m) > 0) THEN 3555 ! For green urban surfaces (e.g. green roofs 3556 ! assume LU short grass 3557 lug_palm = ind_lu_s_grass 3558 IF (lug_palm == ind_lu_user) THEN 3559 message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation' 3560 CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 ) 3561 ELSEIF (lug_palm == ind_lu_b_soil) THEN 3562 lug_dep = 9 3563 ELSEIF (lug_palm == ind_lu_mixed_crops) THEN 3564 lug_dep = 2 3565 ELSEIF (lug_palm == ind_lu_s_grass) THEN 3566 lug_dep = 1 3567 ELSEIF (lug_palm == ind_lu_ev_needle_trees) THEN 3568 lug_dep = 4 3569 ELSEIF (lug_palm == ind_lu_de_needle_trees) THEN 3570 lug_dep = 4 3571 ELSEIF (lug_palm == ind_lu_ev_broad_trees) THEN 3572 lug_dep = 12 3573 ELSEIF (lug_palm == ind_lu_de_broad_trees) THEN 3574 lug_dep = 5 3575 ELSEIF (lug_palm == ind_lu_t_grass) THEN 3576 lug_dep = 1 3577 ELSEIF (lug_palm == ind_lu_desert) THEN 3578 lug_dep = 9 3579 ELSEIF (lug_palm == ind_lu_tundra ) THEN 3580 lug_dep = 8 3581 ELSEIF (lug_palm == ind_lu_irr_crops) THEN 3582 lug_dep = 2 3583 ELSEIF (lug_palm == ind_lu_semidesert) THEN 3584 lug_dep = 8 3585 ELSEIF (lug_palm == ind_lu_ice) THEN 3586 lug_dep = 10 3587 ELSEIF (lug_palm == ind_lu_marsh) THEN 3588 lug_dep = 8 3589 ELSEIF (lug_palm == ind_lu_ev_shrubs) THEN 3590 lug_dep = 14 3591 ELSEIF (lug_palm == ind_lu_de_shrubs ) THEN 3592 lug_dep = 14 3593 ELSEIF (lug_palm == ind_lu_mixed_forest) THEN 3594 lug_dep = 4 3595 ELSEIF (lug_palm == ind_lu_intrup_forest) THEN 3596 lug_dep = 8 3597 ENDIF 3598 ENDIF 3599 3600 IF (surf_usm_h%frac(ind_veg_wall,m) > 0) THEN 3601 ! For walls in USM assume concrete walls/roofs, 3602 ! assumed LU class desert as also assumed for 3603 ! pavements in LSM 3604 luu_palm = ind_lup_conc 3605 IF (luu_palm == ind_lup_user) THEN 3606 message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation' 3607 CALL message( 'chem_depo', 'CM0455', 1, 2, 0, 6, 0 ) 3608 ELSEIF (luu_palm == ind_lup_asph_conc) THEN 3609 luu_dep = 9 3610 ELSEIF (luu_palm == ind_lup_asph) THEN 3611 luu_dep = 9 3612 ELSEIF (luu_palm == ind_lup_conc) THEN 3613 luu_dep = 9 3614 ELSEIF (luu_palm == ind_lup_sett) THEN 3615 luu_dep = 9 3616 ELSEIF (luu_palm == ind_lup_pav_stones) THEN 3617 luu_dep = 9 3618 ELSEIF (luu_palm == ind_lup_cobblest) THEN 3619 luu_dep = 9 3620 ELSEIF (luu_palm == ind_lup_metal) THEN 3621 luu_dep = 9 3622 ELSEIF (luu_palm == ind_lup_wood) THEN 3623 luu_dep = 9 3624 ELSEIF (luu_palm == ind_lup_gravel) THEN 3625 luu_dep = 9 3626 ELSEIF (luu_palm == ind_lup_f_gravel) THEN 3627 luu_dep = 9 3628 ELSEIF (luu_palm == ind_lup_pebblest) THEN 3629 luu_dep = 9 3630 ELSEIF (luu_palm == ind_lup_woodchips) THEN 3631 luu_dep = 9 3632 ELSEIF (luu_palm == ind_lup_tartan) THEN 3633 luu_dep = 9 3634 ELSEIF (luu_palm == ind_lup_art_turf) THEN 3635 luu_dep = 9 3636 ELSEIF (luu_palm == ind_lup_clay) THEN 3637 luu_dep = 9 3638 ENDIF 3639 ENDIF 3640 3641 IF (surf_usm_h%frac(ind_wat_win,m) > 0) THEN 3642 ! For windows in USM assume metal as this is 3643 ! as close as we get, assumed LU class desert 3644 ! as also assumed for pavements in LSM 3645 lud_palm = ind_lup_metal 3646 IF (lud_palm == ind_lup_user) THEN 3647 message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation' 3648 CALL message( 'chem_depo', 'CM0456', 1, 2, 0, 6, 0 ) 3649 ELSEIF (lud_palm == ind_lup_asph_conc) THEN 3650 lud_dep = 9 3651 ELSEIF (lud_palm == ind_lup_asph) THEN 3652 lud_dep = 9 3653 ELSEIF (lud_palm == ind_lup_conc) THEN 3654 lud_dep = 9 3655 ELSEIF (lud_palm == ind_lup_sett) THEN 3656 lud_dep = 9 3657 ELSEIF (lud_palm == ind_lup_pav_stones) THEN 3658 lud_dep = 9 3659 ELSEIF (lud_palm == ind_lup_cobblest) THEN 3660 lud_dep = 9 3661 ELSEIF (lud_palm == ind_lup_metal) THEN 3662 lud_dep = 9 3663 ELSEIF (lud_palm == ind_lup_wood) THEN 3664 lud_dep = 9 3665 ELSEIF (lud_palm == ind_lup_gravel) THEN 3666 lud_dep = 9 3667 ELSEIF (lud_palm == ind_lup_f_gravel) THEN 3668 lud_dep = 9 3669 ELSEIF (lud_palm == ind_lup_pebblest) THEN 3670 lud_dep = 9 3671 ELSEIF (lud_palm == ind_lup_woodchips) THEN 3672 lud_dep = 9 3673 ELSEIF (lud_palm == ind_lup_tartan) THEN 3674 lud_dep = 9 3675 ELSEIF (lud_palm == ind_lup_art_turf) THEN 3676 lud_dep = 9 3677 ELSEIF (lud_palm == ind_lup_clay) THEN 3678 lud_dep = 9 3679 ENDIF 3680 ENDIF 3681 3682 3683 ! TODO: Activate these lines as soon as new ebsolver branch is merged: 3684 ! Set wetness indicator to dry or wet for usm vegetation or pavement 3685 !IF (surf_usm_h%c_liq(m) > 0) THEN 3686 ! nwet = 1 3687 !ELSE 3688 nwet = 0 3689 !ENDIF 3690 3691 ! Compute length of time step 3692 IF ( call_chem_at_all_substeps ) THEN 3693 dt_chem = dt_3d * weight_pres(intermediate_timestep_count) 3694 ELSE 3695 dt_chem = dt_3d 3696 ENDIF 3697 3698 3699 dh = dzw(k) 3700 inv_dh = 1.0_wp / dh 3701 dt_dh = dt_chem/dh 3702 3703 ! Concentration at i,j,k 3704 DO lsp = 1, nspec 3705 cc(lsp) = chem_species(lsp)%conc(k,j,i) 3706 ENDDO 3707 3708 ! Temperature at i,j,k 3709 ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp 3710 3711 ! Surface temperature in degrees Celcius: 3712 ts = ttemp - 273.15 ! in degrees celcius 3713 3714 ! Viscosity of air in lowest layer 3715 visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0) 3716 3717 ! Air density in lowest layer 3718 dens = rho_air_zw(k) 3719 3720 ! Calculate relative humidity from specific humidity for DEPAC 3721 qv_tmp = max(q(k,j,i),0.0_wp) 3722 relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(nzb) ) 3723 3724 3725 3726 ! Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget 3727 ! for each surface fraction. Then derive overall budget taking into account the surface fractions. 3728 3729 ! Walls/roofs 3730 IF (surf_usm_h%frac(ind_veg_wall,m) > 0) THEN 3731 3732 3733 ! No vegetation on non-green walls: 3734 lai = 0.0_wp 3735 sai = 0.0_wp 3736 3737 slinnfac = 1.0_wp 3738 3739 ! Get vd 3740 DO lsp = 1, nvar 3741 !Initialize 3742 vs = 0.0_wp 3743 vd_lu = 0.0_wp 3744 Rs = 0.0_wp 3745 Rb = 0.0_wp 3746 Rc_tot = 0.0_wp 3747 IF (spc_names(lsp) == 'PM10') THEN 3748 part_type = 1 3749 ! sedimentation velocity in lowest layer 3750 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3751 particle_pars(ind_p_size, part_type), & 3752 particle_pars(ind_p_slip, part_type), & 3753 visc) 3754 3755 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 3756 vs, & 3757 particle_pars(ind_p_size, part_type), & 3758 particle_pars(ind_p_slip, part_type), & 3759 nwet, ttemp, dens, visc, & 3760 luu_dep, & 3761 r_aero_lu, ustar_lu) 3762 3763 bud_now_luu(lsp) = - cc(lsp) * & 3764 (1.0 - exp(-vd_lu*dt_dh ))*dh 3765 3766 3767 ELSEIF (spc_names(lsp) == 'PM25') THEN 3768 part_type = 2 3769 ! sedimentation velocity in lowest layer: 3770 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3771 particle_pars(ind_p_size, part_type), & 3772 particle_pars(ind_p_slip, part_type), & 3773 visc) 3774 3775 3776 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 3777 vs, & 3778 particle_pars(ind_p_size, part_type), & 3779 particle_pars(ind_p_slip, part_type), & 3780 nwet, ttemp, dens, visc, & 3781 luu_dep , & 3782 r_aero_lu, ustar_lu) 3783 3784 bud_now_luu(lsp) = - cc(lsp) * & 3785 (1.0 - exp(-vd_lu*dt_dh ))*dh 3786 3787 3788 ELSE 3789 3790 ! GASES 3791 3792 ! Read spc_name of current species for gas parameter 3793 3794 IF (ANY(pspecnames(:) == spc_names(lsp))) THEN 3795 i_pspec = 0 3796 DO pspec = 1, 69 3797 IF (pspecnames(pspec) == spc_names(lsp)) THEN 3798 i_pspec = pspec 3799 END IF 3800 ENDDO 3801 3802 ELSE 3803 ! Species for now not deposited 3804 CYCLE 3805 ENDIF 3806 3807 ! Factor used for conversion from ppb to ug/m3 : 3808 ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 3809 ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 3810 ! c 1e-9 xm_tracer 1e9 / xm_air dens 3811 ! thus: 3812 ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 3813 ! Use density at lowest layer: 3814 3815 ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 3816 3817 ! atmospheric concentration in DEPAC is requested in ug/m3: 3818 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3819 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 3820 3821 ! Diffusivity for DEPAC relevant gases 3822 ! Use default value 3823 diffc = 0.11e-4 3824 ! overwrite with known coefficients of diffusivity from Massman (1998) 3825 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 3826 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 3827 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 3828 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 3829 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 3830 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 3831 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 3832 3833 3834 ! Get quasi-laminar boundary layer resistance Rb: 3835 CALL get_rb_cell( (luu_dep == ilu_water_sea) .or. (luu_dep == ilu_water_inland), & 3836 z0h_lu, ustar_lu, diffc, & 3837 Rb) 3838 3839 ! Get Rc_tot 3840 CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & 3841 relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3842 r_aero_lu , Rb) 3843 3844 3845 ! Calculate budget 3846 IF (Rc_tot .le. 0.0) THEN 3847 3848 bud_now_luu(lsp) = 0.0_wp 3849 3850 ELSE 3851 3852 ! Compute exchange velocity for current lu: 3853 vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) 3854 3855 bud_now_luu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & 3856 (1.0 - exp(-vd_lu*dt_dh ))*dh 3857 ENDIF 3858 3859 ENDIF 3860 ENDDO 3861 ENDIF 3862 3863 3864 ! Green usm surfaces 3865 IF (surf_usm_h%frac(ind_pav_green,m) > 0) THEN 3866 3867 3868 slinnfac = 1.0_wp 3869 3870 ! Get vd 3871 DO lsp = 1, nvar 3872 !Initialize 3873 vs = 0.0_wp 3874 vd_lu = 0.0_wp 3875 Rs = 0.0_wp 3876 Rb = 0.0_wp 3877 Rc_tot = 0.0_wp 3878 IF (spc_names(lsp) == 'PM10') THEN 3879 part_type = 1 3880 ! sedimentation velocity in lowest layer: 3881 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3882 particle_pars(ind_p_size, part_type), & 3883 particle_pars(ind_p_slip, part_type), & 3884 visc) 3885 3886 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 3887 vs, & 3888 particle_pars(ind_p_size, part_type), & 3889 particle_pars(ind_p_slip, part_type), & 3890 nwet, ttemp, dens, visc, & 3891 lug_dep, & 3892 r_aero_lu, ustar_lu) 3893 3894 bud_now_lug(lsp) = - cc(lsp) * & 3895 (1.0 - exp(-vd_lu*dt_dh ))*dh 3896 3897 3898 ELSEIF (spc_names(lsp) == 'PM25') THEN 3899 part_type = 2 3900 ! sedimentation velocity in lowest layer: 3901 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3902 particle_pars(ind_p_size, part_type), & 3903 particle_pars(ind_p_slip, part_type), & 3904 visc) 3905 3906 3907 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 3908 vs, & 3909 particle_pars(ind_p_size, part_type), & 3910 particle_pars(ind_p_slip, part_type), & 3911 nwet, ttemp, dens, visc, & 3912 lug_dep, & 3913 r_aero_lu, ustar_lu) 3914 3915 bud_now_lug(lsp) = - cc(lsp) * & 3916 (1.0 - exp(-vd_lu*dt_dh ))*dh 3917 3918 3919 ELSE 3920 3921 ! GASES 3922 3923 ! Read spc_name of current species for gas parameter 3924 3925 IF (ANY(pspecnames(:) == spc_names(lsp))) THEN 3926 i_pspec = 0 3927 DO pspec = 1, 69 3928 IF (pspecnames(pspec) == spc_names(lsp)) THEN 3929 i_pspec = pspec 3930 END IF 3931 ENDDO 3932 3933 ELSE 3934 ! Species for now not deposited 3935 CYCLE 3936 ENDIF 3937 3938 ! Factor used for conversion from ppb to ug/m3 : 3939 ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 3940 ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 3941 ! c 1e-9 xm_tracer 1e9 / xm_air dens 3942 ! thus: 3943 ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 3944 ! Use density at lowest layer: 3945 3946 ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 3947 3948 ! atmospheric concentration in DEPAC is requested in ug/m3: 3949 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3950 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 3951 3952 ! Diffusivity for DEPAC relevant gases 3953 ! Use default value 3954 diffc = 0.11e-4 3955 ! overwrite with known coefficients of diffusivity from Massman (1998) 3956 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 3957 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 3958 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 3959 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 3960 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 3961 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 3962 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 3963 3964 3965 ! Get quasi-laminar boundary layer resistance Rb: 3966 CALL get_rb_cell( (lug_dep == ilu_water_sea) .or. (lug_dep == ilu_water_inland), & 3967 z0h_lu, ustar_lu, diffc, & 3968 Rb) 3969 3970 3971 ! Get Rc_tot 3972 CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & 3973 relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3974 r_aero_lu , Rb) 3975 3976 3977 ! Calculate budget 3978 IF (Rc_tot .le. 0.0) THEN 3979 3980 bud_now_lug(lsp) = 0.0_wp 3981 3982 ELSE 3983 3984 ! Compute exchange velocity for current lu: 3985 vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) 3986 3987 bud_now_lug(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & 3988 (1.0 - exp(-vd_lu*dt_dh ))*dh 3989 ENDIF 3990 3991 3992 ENDIF 3993 ENDDO 3994 ENDIF 3995 3996 3997 ! Windows 3998 IF (surf_usm_h%frac(ind_wat_win,m) > 0) THEN 3999 4000 4001 ! No vegetation on windows: 4002 lai = 0.0_wp 4003 sai = 0.0_wp 4004 4005 slinnfac = 1.0_wp 4006 4007 ! Get vd 4008 DO lsp = 1, nvar 4009 !Initialize 4010 vs = 0.0_wp 4011 vd_lu = 0.0_wp 4012 Rs = 0.0_wp 4013 Rb = 0.0_wp 4014 Rc_tot = 0.0_wp 4015 IF (spc_names(lsp) == 'PM10') THEN 4016 part_type = 1 4017 ! sedimentation velocity in lowest layer: 4018 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 4019 particle_pars(ind_p_size, part_type), & 4020 particle_pars(ind_p_slip, part_type), & 4021 visc) 4022 4023 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 4024 vs, & 4025 particle_pars(ind_p_size, part_type), & 4026 particle_pars(ind_p_slip, part_type), & 4027 nwet, ttemp, dens, visc, & 4028 lud_dep, & 4029 r_aero_lu, ustar_lu) 4030 4031 bud_now_lud(lsp) = - cc(lsp) * & 4032 (1.0 - exp(-vd_lu*dt_dh ))*dh 4033 4034 4035 ELSEIF (spc_names(lsp) == 'PM25') THEN 4036 part_type = 2 4037 ! sedimentation velocity in lowest layer: 4038 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 4039 particle_pars(ind_p_size, part_type), & 4040 particle_pars(ind_p_slip, part_type), & 4041 visc) 4042 4043 4044 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 4045 vs, & 4046 particle_pars(ind_p_size, part_type), & 4047 particle_pars(ind_p_slip, part_type), & 4048 nwet, ttemp, dens, visc, & 4049 lud_dep, & 4050 r_aero_lu, ustar_lu) 4051 4052 bud_now_lud(lsp) = - cc(lsp) * & 4053 (1.0 - exp(-vd_lu*dt_dh ))*dh 4054 4055 4056 ELSE 4057 4058 ! GASES 4059 4060 ! Read spc_name of current species for gas PARAMETER 4061 4062 IF (ANY(pspecnames(:) == spc_names(lsp))) THEN 4063 i_pspec = 0 4064 DO pspec = 1, 69 4065 IF (pspecnames(pspec) == spc_names(lsp)) THEN 4066 i_pspec = pspec 4067 END IF 4068 ENDDO 4069 4070 ELSE 4071 ! Species for now not deposited 4072 CYCLE 4073 ENDIF 4074 4075 ! Factor used for conversion from ppb to ug/m3 : 4076 ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 4077 ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 4078 ! c 1e-9 xm_tracer 1e9 / xm_air dens 4079 ! thus: 4080 ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 4081 ! Use density at lowest layer: 4082 4083 ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 4084 4085 ! atmospheric concentration in DEPAC is requested in ug/m3: 4086 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 4087 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 4088 4089 ! Diffusivity for DEPAC relevant gases 4090 ! Use default value 4091 diffc = 0.11e-4 4092 ! overwrite with known coefficients of diffusivity from Massman (1998) 4093 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 4094 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 4095 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 4096 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 4097 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 4098 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 4099 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 4100 4101 4102 ! Get quasi-laminar boundary layer resistance Rb: 4103 CALL get_rb_cell( (lud_dep == ilu_water_sea) .or. (lud_dep == ilu_water_inland), & 4104 z0h_lu, ustar_lu, diffc, & 4105 Rb) 4106 4107 ! Get Rc_tot 4108 CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & 4109 relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 4110 r_aero_lu , Rb) 4111 4112 4113 ! Calculate budget 4114 IF (Rc_tot .le. 0.0) THEN 4115 4116 bud_now_lud(lsp) = 0.0_wp 4117 4118 ELSE 4119 4120 ! Compute exchange velocity for current lu: 4121 vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) 4122 4123 bud_now_lud(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & 4124 (1.0 - exp(-vd_lu*dt_dh ))*dh 4125 ENDIF 4126 4127 ENDIF 4128 ENDDO 4129 ENDIF 4130 4131 4132 bud_now = 0.0_wp 4133 ! Calculate budget for surface m and adapt concentration 4134 DO lsp = 1, nspec 4135 4136 4137 bud_now(lsp) = surf_usm_h%frac(ind_veg_wall,m)*bud_now_luu(lsp) + & 4138 surf_usm_h%frac(ind_pav_green,m)*bud_now_lug(lsp) + & 4139 surf_usm_h%frac(ind_wat_win,m)*bud_now_lud(lsp) 4140 4141 ! Compute new concentration, add contribution of all exchange processes: 4142 cc(lsp) = cc(lsp) + bud_now(lsp)*inv_dh 4143 4144 chem_species(lsp)%conc(k,j,i) = max(0.0_wp, cc(lsp)) 4145 4146 ENDDO 4147 4148 ENDIF 4149 4150 4151 4152 END SUBROUTINE chem_depo 4153 4154 4155 4156 !---------------------------------------------------------------------------------- 4157 ! 4158 ! DEPAC: 4159 ! Code of the DEPAC routine and corresponding subroutines below from the DEPAC 4160 ! module of the LOTOS-EUROS model (Manders et al., 2017) 4161 ! 4162 ! Original DEPAC routines by RIVM and TNO (2015), for Documentation see 4163 ! van Zanten et al., 2010. 4164 !--------------------------------------------------------------------------------- 4165 ! 4166 !---------------------------------------------------------------------------------- 4167 ! 4168 ! depac: compute total canopy (or surface) resistance Rc for gases 4169 !---------------------------------------------------------------------------------- 4170 4171 SUBROUTINE drydepos_gas_depac(compnam, day_of_year, lat, t, ust, glrad, sinphi, & 4172 rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, catm, diffc, & 4173 ra, rb) 4174 4175 4176 ! The last four rows of depac arguments are OPTIONAL: 4177 ! 4178 ! A. compute Rc_tot without compensation points (ccomp_tot will be zero): 4179 ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi]) 4180 ! 4181 ! B. compute Rc_tot with compensation points (used for LOTOS-EUROS): 4182 ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4183 ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water) 4184 ! 4185 ! C. compute effective Rc based on compensation points (used for OPS): 4186 ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4187 ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & 4188 ! ra, rb, rc_eff) 4189 ! X1. Extra (OPTIONAL) output variables: 4190 ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4191 ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & 4192 ! ra, rb, rc_eff, & 4193 ! gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out) 4194 ! X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves): 4195 ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4196 ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & 4197 ! ra, rb, rc_eff, & 4198 ! gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, & 4199 ! calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA) 4200 4201 4202 IMPLICIT NONE 4203 4204 CHARACTER(len=*) , INTENT(in) :: compnam ! component name 4205 ! 'HNO3','NO','NO2','O3','SO2','NH3' 4206 INTEGER(iwp) , INTENT(in) :: day_of_year ! day of year, 1 ... 365 (366) 4207 REAL(kind=wp) , INTENT(in) :: lat ! latitude Northern hemisphere (degrees) (DEPAC cannot be used for S. hemisphere) 4208 REAL(kind=wp) , INTENT(in) :: t ! temperature (C) 4209 ! NB discussion issue is temp T_2m or T_surf or T_leaf? 4210 REAL(kind=wp) , INTENT(in) :: ust ! friction velocity (m/s) 4211 REAL(kind=wp) , INTENT(in) :: glrad ! global radiation (W/m2) 4212 REAL(kind=wp) , INTENT(in) :: sinphi ! sin of solar elevation angle 4213 REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) 4214 REAL(kind=wp) , INTENT(in) :: lai ! one-sidedleaf area index (-) 4215 REAL(kind=wp) , INTENT(in) :: sai ! surface area index (-) (lai + branches and stems) 4216 REAL(kind=wp) , INTENT(in) :: diffc ! Diffusivity 4217 INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow 4218 INTEGER(iwp) , INTENT(in) :: lu ! land use type, lu = 1,...,nlu 4219 INTEGER(iwp) , INTENT(in) :: iratns ! index for NH3/SO2 ratio used for SO2; 4220 ! iratns = 1: low NH3/SO2 4221 ! iratns = 2: high NH3/SO2 4222 ! iratns = 3: very low NH3/SO2 4223 REAL(kind=wp) , INTENT(out) :: rc_tot ! total canopy resistance Rc (s/m) 4224 REAL(kind=wp) , INTENT(out) :: ccomp_tot ! total compensation point (ug/m3) [= 0 for species that don't have a compensation 4225 REAL(kind=wp) ,INTENT(in) :: p ! pressure (Pa) 4226 REAL(kind=wp), INTENT(in) :: catm ! actual atmospheric concentration (ug/m3) 4227 REAL(kind=wp), INTENT(in) :: ra ! aerodynamic resistance (s/m) 4228 REAL(kind=wp), INTENT(in) :: rb ! boundary layer resistance (s/m) 4229 4230 ! Local variables: 4231 REAL(kind=wp) :: laimax ! maximum leaf area index (-) 4232 LOGICAL :: ready ! Rc has been set 4233 ! = 1 -> constant Rc 4234 ! = 2 -> temperature dependent Rc 4235 REAL(kind=wp) :: gw ! external leaf conductance (m/s) 4236 REAL(kind=wp) :: gstom ! stomatal conductance (m/s) 4237 REAL(kind=wp) :: gsoil_eff ! effective soil conductance (m/s) 4238 REAL(kind=wp) :: gc_tot ! total canopy conductance (m/s) 4239 REAL(kind=wp) :: cw ! external leaf surface compensation point (ug/m3) 4240 REAL(kind=wp) :: cstom ! stomatal compensation point (ug/m3) 4241 REAL(kind=wp) :: csoil ! soil compensation point (ug/m3) 4242 4243 ! Vegetation indicators: 4244 LOGICAL :: LAI_present ! leaves are present for current land use type 4245 LOGICAL :: SAI_present ! vegetation is present for current land use type 4246 4247 ! Component number taken from component name, paramteres matched with include files 4248 INTEGER(iwp) :: icmp 4249 4250 ! component numbers: 4251 INTEGER(iwp), PARAMETER :: icmp_o3 = 1 4252 INTEGER(iwp), PARAMETER :: icmp_so2 = 2 4253 INTEGER(iwp), PARAMETER :: icmp_no2 = 3 4254 INTEGER(iwp), PARAMETER :: icmp_no = 4 4255 INTEGER(iwp), PARAMETER :: icmp_nh3 = 5 4256 INTEGER(iwp), PARAMETER :: icmp_co = 6 4257 INTEGER(iwp), PARAMETER :: icmp_no3 = 7 4258 INTEGER(iwp), PARAMETER :: icmp_hno3 = 8 4259 INTEGER(iwp), PARAMETER :: icmp_n2o5 = 9 4260 INTEGER(iwp), PARAMETER :: icmp_h2o2 = 10 4261 4262 4263 4264 ! Define component number 4265 SELECT CASE ( trim(compnam) ) 4266 4267 CASE ( 'O3', 'o3' ) 4268 icmp = icmp_o3 4269 4270 CASE ( 'SO2', 'so2' ) 4271 icmp = icmp_so2 4272 4273 CASE ( 'NO2', 'no2' ) 4274 icmp = icmp_no2 4275 4276 CASE ( 'NO', 'no' ) 4277 icmp = icmp_no 4278 4279 CASE ( 'NH3', 'nh3' ) 4280 icmp = icmp_nh3 4281 4282 CASE ( 'CO', 'co' ) 4283 icmp = icmp_co 4284 4285 CASE ( 'NO3', 'no3' ) 4286 icmp = icmp_no3 4287 4288 CASE ( 'HNO3', 'hno3' ) 4289 icmp = icmp_hno3 4290 4291 CASE ( 'N2O5', 'n2o5' ) 4292 icmp = icmp_n2o5 4293 4294 CASE ( 'H2O2', 'h2o2' ) 4295 icmp = icmp_h2o2 4296 4297 CASE default 4298 !Component not part of DEPAC --> not deposited 4299 RETURN 4300 4301 END SELECT 4302 4303 ! inititalize 4304 gw = 0. 4305 gstom = 0. 4306 gsoil_eff = 0. 4307 gc_tot = 0. 4308 cw = 0. 4309 cstom = 0. 4310 csoil = 0. 4311 4312 4313 ! Check whether vegetation is present (in that CASE the leaf or surface area index > 0): 4314 LAI_present = (lai .gt. 0.0) 4315 SAI_present = (sai .gt. 0.0) 4316 4317 ! Set Rc (i.e. rc_tot) in special cases: 4318 CALL rc_special(icmp,compnam,lu,t, nwet,rc_tot,ready,ccomp_tot) 4319 4320 ! set effective resistance equal to total resistance, this will be changed for compensation points 4321 !IF ( present(rc_eff) ) then 4322 ! rc_eff = rc_tot 4323 !END IF 4324 4325 ! IF Rc is not set: 4326 IF (.not. ready) then 4327 4328 ! External conductance: 4329 CALL rc_gw(compnam,iratns,t,rh,nwet,SAI_present,sai,gw) 4330 4331 ! Stomatal conductance: 4332 ! CALL rc_gstom(icmp,compnam,lu,LAI_present,lai,glrad,sinphi,t,rh,diffc,gstom,p,& 4333 ! smi=smi,calc_stom_o3flux=calc_stom_o3flux) 4334 CALL rc_gstom(icmp,compnam,lu,LAI_present,lai,glrad,sinphi,t,rh,diffc,gstom,p) 4335 ! Effective soil conductance: 4336 CALL rc_gsoil_eff(icmp,lu,sai,ust,nwet,t,gsoil_eff) 4337 4338 ! Total canopy conductance (gc_tot) and resistance Rc (rc_tot): 4339 CALL rc_rctot(gstom,gsoil_eff,gw,gc_tot,rc_tot) 4340 4341 ! Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3): 4342 ! CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,& 4343 ! LAI_present, SAI_present, & 4344 ! ccomp_tot, & 4345 ! catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, & 4346 ! c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, & 4347 ! tsea=tsea,cw=cw,cstom=cstom,csoil=csoil ) 4348 4349 ! Effective Rc based on compensation points: 4350 ! IF ( present(rc_eff) ) then 4351 ! check on required arguments: 4352 !IF ( (.not. present(catm)) .or. (.not. present(ra)) .or. (.not. present(rb)) ) then 4353 ! stop 'output argument rc_eff requires input arguments catm, ra and rb' 4354 !END IF 4355 ! compute rc_eff : 4356 !CALL rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff) 4357 !ENDIF 4358 ENDIF 4359 4360 END SUBROUTINE drydepos_gas_depac 4361 4362 4363 4364 !------------------------------------------------------------------- 4365 ! rc_special: compute total canopy resistance in special CASEs 4366 !------------------------------------------------------------------- 4367 SUBROUTINE rc_special(icmp,compnam,lu,t,nwet,rc_tot,ready,ccomp_tot) 4368 4369 INTEGER(iwp) , INTENT(in) :: icmp ! component index 4370 CHARACTER(len=*) , INTENT(in) :: compnam ! component name 4371 INTEGER(iwp) , INTENT(in) :: lu ! land use type, lu = 1,...,nlu 4372 REAL(kind=wp) , INTENT(in) :: t ! temperature (C) 4373 ! = 1 -> constant Rc 4374 ! = 2 -> temperature dependent Rc 4375 INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow 4376 REAL(kind=wp) , INTENT(out) :: rc_tot ! total canopy resistance Rc (s/m) 4377 LOGICAL , INTENT(out) :: ready ! Rc has been set 4378 REAL(kind=wp) , INTENT(out) :: ccomp_tot ! total compensation point (ug/m3) 4379 4380 ! rc_tot is not yet set: 4381 ready = .false. 4382 4383 ! Default compensation point in special CASEs = 0: 4384 ccomp_tot = 0.0 4385 4386 SELECT CASE(trim(compnam)) 4387 CASE('HNO3','N2O5','NO3','H2O2') 4388 ! No separate resistances for HNO3; just one total canopy resistance: 4389 IF (t .lt. -5.0 .and. nwet .eq. 9) then 4390 ! T < 5 C and snow: 4391 rc_tot = 50. 4392 ELSE 4393 ! all other circumstances: 4394 rc_tot = 10.0 4395 ENDIF 4396 ready = .true. 4397 4398 CASE('NO','CO') 4399 IF (lu .eq. ilu_water_sea .or. lu .eq. ilu_water_inland) then ! water 4400 rc_tot = 2000. 4401 ready = .true. 4402 ELSEIF (nwet .eq. 1) then ! wet 4403 rc_tot = 2000. 4404 ready = .true. 4405 ENDIF 4406 CASE('NO2','O3','SO2','NH3') 4407 ! snow surface: 4408 IF (nwet.eq.9) then 4409 !CALL rc_snow(ipar_snow(icmp),t,rc_tot) 4410 ready = .true. 4411 ENDIF 4412 CASE default 4413 message_string = 'Component '// trim(compnam) // ' not supported' 4414 CALL message( 'rc_special', 'CM0457', 1, 2, 0, 6, 0 ) 4415 END SELECT 4416 4417 END SUBROUTINE rc_special 4418 4419 4420 4421 !------------------------------------------------------------------- 4422 ! rc_gw: compute external conductance 4423 !------------------------------------------------------------------- 4424 SUBROUTINE rc_gw( compnam, iratns,t,rh,nwet, SAI_present, sai, gw ) 4425 4426 ! Input/output variables: 4427 CHARACTER(len=*) , INTENT(in) :: compnam ! component name 4428 INTEGER(iwp) , INTENT(in) :: iratns ! index for NH3/SO2 ratio; 4429 ! iratns = 1: low NH3/SO2 4430 ! iratns = 2: high NH3/SO2 4431 ! iratns = 3: very low NH3/SO2 4432 REAL(kind=wp) , INTENT(in) :: t ! temperature (C) 4433 REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) 4434 INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow 4435 LOGICAL , INTENT(in) :: SAI_present 4436 REAL(kind=wp) , INTENT(in) :: sai ! one-sided leaf area index (-) 4437 REAL(kind=wp) , INTENT(out) :: gw ! external leaf conductance (m/s) 4438 4439 SELECT CASE(trim(compnam)) 4440 4441 !CASE('HNO3','N2O5','NO3','H2O2') this routine is not CALLed for HNO3 4442 4443 CASE('NO2') 4444 CALL rw_constant(2000.0_wp,SAI_present,gw) 4445 4446 CASE('NO','CO') 4447 CALL rw_constant(-9999.0_wp,SAI_present,gw) ! see Erisman et al, 1994 section 3.2.3 4448 4449 CASE('O3') 4450 CALL rw_constant(2500.0_wp,SAI_present,gw) 4451 4452 CASE('SO2') 4453 CALL rw_so2( t, nwet, rh, iratns, SAI_present, gw ) 4454 4455 CASE('NH3') 4456 CALL rw_nh3_sutton(t,rh,SAI_present,gw) 4457 4458 ! conversion from leaf resistance to canopy resistance by multiplying with SAI: 4459 Gw = sai*gw 4460 4461 CASE default 4462 message_string = 'Component '// trim(compnam) // ' not supported' 4463 CALL message( 'rc_gw', 'CM0458', 1, 2, 0, 6, 0 ) 4464 END SELECT 4465 4466 END SUBROUTINE rc_gw 4467 4468 4469 4470 !------------------------------------------------------------------- 4471 ! rw_so2: compute external leaf conductance for SO2 4472 !------------------------------------------------------------------- 4473 SUBROUTINE rw_so2( t, nwet, rh, iratns, SAI_present, gw ) 4474 4475 ! Input/output variables: 4476 REAL(kind=wp) , INTENT(in) :: t ! temperature (C) 4477 INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow 4478 REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) 4479 INTEGER(iwp) , INTENT(in) :: iratns ! index for NH3/SO2 ratio; 4480 ! iratns = 1: low NH3/SO2 4481 ! iratns = 2: high NH3/SO2 4482 ! iratns = 3: very low NH3/SO2 4483 LOGICAL, INTENT(in) :: SAI_present 4484 REAL(kind=wp) , INTENT(out) :: gw ! external leaf conductance (m/s) 4485 4486 ! Variables from module: 4487 ! SAI_present: vegetation is present 4488 4489 ! Local variables: 4490 REAL(kind=wp) :: rw ! external leaf resistance (s/m) 4491 4492 ! Check IF vegetation present: 4493 IF (SAI_present) then 4494 4495 IF (nwet .eq. 0) then 4496 !-------------------------- 4497 ! dry surface 4498 !-------------------------- 4499 ! T > -1 C 4500 IF (t .gt. -1.0) then 4501 IF (rh .lt. 81.3) then 4502 rw = 25000*exp(-0.0693*rh) 4503 ELSE 4504 rw = 0.58e12*exp(-0.278*rh) + 10. 4505 ENDIF 4506 ELSE 4507 ! -5 C < T <= -1 C 4508 IF (t .gt. -5.0) then 4509 rw=200 4510 ELSE 4511 ! T <= -5 C 4512 rw=500 4513 ENDIF 4514 ENDIF 4515 ELSE 4516 !-------------------------- 4517 ! wet surface 4518 !-------------------------- 4519 rw = 10. !see Table 5, Erisman et al, 1994 Atm. Environment, 0 is impl. as 10 4520 ENDIF 4521 4522 ! very low NH3/SO2 ratio: 4523 IF (iratns == iratns_very_low ) rw = rw+50. 4524 4525 ! Conductance: 4526 gw = 1./rw 4527 ELSE 4528 ! no vegetation: 4529 gw = 0.0 4530 ENDIF 4531 4532 END SUBROUTINE rw_so2 4533 4534 4535 4536 !------------------------------------------------------------------- 4537 ! rw_nh3_sutton: compute external leaf conductance for NH3, 4538 ! following Sutton & Fowler, 1993 4539 !------------------------------------------------------------------- 4540 SUBROUTINE rw_nh3_sutton(tsurf,rh,SAI_present,gw) 4541 4542 ! Input/output variables: 4543 REAL(kind=wp) , INTENT(in) :: tsurf ! surface temperature (C) 4544 REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) 4545 LOGICAL, INTENT(in) :: SAI_present 4546 REAL(kind=wp) , INTENT(out) :: gw ! external leaf conductance (m/s) 4547 4548 ! Variables from module: 4549 ! SAI_present: vegetation is present 4550 4551 ! Local variables: 4552 REAL(kind=wp) :: rw ! external leaf resistance (s/m) 4553 REAL(kind=wp) :: sai_grass_haarweg ! surface area index at experimental site Haarweg 4554 4555 ! Fix SAI_grass at value valid for Haarweg data for which gamma_w parametrization is derived 4556 sai_grass_haarweg = 3.5 4557 4558 ! 100 - rh 4559 ! rw = 2.0 * exp(----------) 4560 ! 12 4561 4562 IF (SAI_present) then 4563 4564 ! External resistance according to Sutton & Fowler, 1993 4565 rw = 2.0 * exp((100.0 - rh)/12.0) 4566 rw = sai_grass_haarweg * rw 4567 4568 ! Frozen soil (from Depac v1): 4569 IF (tsurf .lt. 0) rw = 200 4570 4571 ! Conductance: 4572 gw = 1./rw 4573 ELSE 4574 ! no vegetation: 4575 gw = 0.0 4576 ENDIF 4577 4578 END SUBROUTINE rw_nh3_sutton 4579 4580 4581 4582 !------------------------------------------------------------------- 4583 ! rw_constant: compute constant external leaf conductance 4584 !------------------------------------------------------------------- 4585 SUBROUTINE rw_constant(rw_val,SAI_present,gw) 4586 4587 ! Input/output variables: 4588 REAL(kind=wp) , INTENT(in) :: rw_val ! constant value of Rw 4589 LOGICAL , INTENT(in) :: SAI_present 4590 REAL(kind=wp) , INTENT(out) :: gw ! wernal leaf conductance (m/s) 4591 4592 ! Variables from module: 4593 ! SAI_present: vegetation is present 4594 4595 ! Compute conductance: 4596 IF (SAI_present .and. .not.missing(rw_val)) then 4597 gw = 1./rw_val 4598 ELSE 4599 gw = 0. 4600 ENDIF 4601 4602 END SUBROUTINE rw_constant 4603 4604 4605 4606 !------------------------------------------------------------------- 4607 ! rc_gstom: compute stomatal conductance 4608 !------------------------------------------------------------------- 4609 SUBROUTINE rc_gstom( icmp, compnam, lu, LAI_present, lai, glrad, sinphi, t, rh, diffc, & 4610 gstom, & 4611 p) 4612 4613 ! input/output variables: 4614 INTEGER(iwp), INTENT(in) :: icmp ! component index 4615 CHARACTER(len=*), INTENT(in) :: compnam ! component name 4616 INTEGER(iwp), INTENT(in) :: lu ! land use type , lu = 1,...,nlu 4617 LOGICAL, INTENT(in) :: LAI_present 4618 REAL(kind=wp), INTENT(in) :: lai ! one-sided leaf area index 4619 REAL(kind=wp), INTENT(in) :: glrad ! global radiation (W/m2) 4620 REAL(kind=wp), INTENT(in) :: sinphi ! sin of solar elevation angle 4621 REAL(kind=wp), INTENT(in) :: t ! temperature (C) 4622 REAL(kind=wp), INTENT(in) :: rh ! relative humidity (%) 4623 REAL(kind=wp), INTENT(in) :: diffc ! diffusion coefficient of the gas involved 4624 REAL(kind=wp), INTENT(out):: gstom ! stomatal conductance (m/s) 4625 REAL(kind=wp), OPTIONAL,INTENT(in) :: p ! pressure (Pa) 4626 4627 4628 ! Local variables 4629 REAL(kind=wp) :: vpd ! vapour pressure deficit (kPa) 4630 4631 REAL(kind=wp), PARAMETER :: dO3 = 0.13e-4 ! diffusion coefficient of ozon (m2/s) 4632 4633 4634 SELECT CASE(trim(compnam)) 4635 4636 !CASE('HNO3','N2O5','NO3','H2O2') this routine is not CALLed for HNO3 4637 4638 CASE('NO','CO') 4639 ! for NO stomatal uptake is neglected: 4640 gstom = 0.0 4641 4642 CASE('NO2','O3','SO2','NH3') 4643 4644 ! IF vegetation present: 4645 IF (LAI_present) then 4646 4647 IF (glrad .gt. 0.0) then 4648 CALL rc_get_vpd(t,rh,vpd) 4649 CALL rc_gstom_emb( lu, glrad, t, vpd, LAI_present, lai, sinphi, gstom, p ) 4650 gstom = gstom*diffc/dO3 ! Gstom of Emberson is derived for ozone 4651 ELSE 4652 gstom = 0.0 4653 ENDIF 4654 ELSE 4655 ! no vegetation; zero conductance (infinite resistance): 4656 gstom = 0.0 4657 ENDIF 4658 4659 CASE default 4660 message_string = 'Component '// trim(compnam) // ' not supported' 4661 CALL message( 'rc_gstom', 'CM0459', 1, 2, 0, 6, 0 ) 4662 END SELECT 4663 4664 END SUBROUTINE rc_gstom 4665 4666 4667 4668 !------------------------------------------------------------------- 4669 ! rc_gstom_emb: stomatal conductance according to Emberson 4670 !------------------------------------------------------------------- 4671 SUBROUTINE rc_gstom_emb( lu, glrad, T, vpd, LAI_present, lai, sinp, & 4672 Gsto, p ) 4673 ! History 4674 ! Original code from Lotos-Euros, TNO, M. Schaap 4675 ! 2009-08, M.C. van Zanten, Rivm 4676 ! Updated and extended. 4677 ! 2009-09, Arjo Segers, TNO 4678 ! Limitted temperature influence to range to avoid 4679 ! floating point exceptions. 4680 ! 4681 ! Method 4682 ! 4683 ! Code based on Emberson et al, 2000, Env. Poll., 403-413 4684 ! Notation conform Unified EMEP Model Description Part 1, ch 8 4685 ! 4686 ! In the calculation of F_light the modification of L. Zhang 2001, AE to the PARshade and PARsun 4687 ! parametrizations of Norman 1982 are applied 4688 ! 4689 ! F_phen and F_SWP are set to 1 4690 ! 4691 ! Land use types DEPAC versus Emberson (Table 5.1, EMEP model description) 4692 ! DEPAC Emberson 4693 ! 1 = grass GR = grassland 4694 ! 2 = arable land TC = temperate crops ( LAI according to RC = rootcrops) 4695 ! 3 = permanent crops TC = temperate crops ( LAI according to RC = rootcrops) 4696 ! 4 = coniferous forest CF = tempareate/boREAL(kind=wp) coniferous forest 4697 ! 5 = deciduous forest DF = temperate/boREAL(kind=wp) deciduous forest 4698 ! 6 = water W = water 4699 ! 7 = urban U = urban 4700 ! 8 = other GR = grassland 4701 ! 9 = desert DE = desert 4702 ! 4703 4704 ! Emberson specific declarations 4705 4706 ! Input/output variables: 4707 INTEGER(iwp), INTENT(in) :: lu ! land use type, lu = 1,...,nlu 4708 REAL(kind=wp), INTENT(in) :: glrad ! global radiation (W/m2) 4709 REAL(kind=wp), INTENT(in) :: T ! temperature (C) 4710 REAL(kind=wp), INTENT(in) :: vpd ! vapour pressure deficit (kPa) 4711 LOGICAL, INTENT(in) :: LAI_present 4712 REAL(kind=wp), INTENT(in) :: lai ! one-sided leaf area index 4713 REAL(kind=wp), INTENT(in) :: sinp ! sin of solar elevation angle 4714 REAL(kind=wp), INTENT(out):: Gsto ! stomatal conductance (m/s) 4715 REAL(kind=wp), OPTIONAL, INTENT(in) :: p ! pressure (Pa) 4716 4717 ! Local variables: 4718 REAL(kind=wp) :: F_light, F_phen, F_temp, F_vpd, F_swp, bt, F_env 4719 REAL(kind=wp) :: pardir, pardiff, parshade, parsun, LAIsun, LAIshade, sinphi 4720 REAL(kind=wp) :: pres 4721 REAL(kind=wp), PARAMETER :: p_sealevel = 1.01325e5 ! Pa 4722 4723 4724 ! Check whether vegetation is present: 4725 IF (LAI_present) then 4726 4727 ! calculation of correction factors for stomatal conductance 4728 IF (sinp .le. 0.0) then 4729 sinphi = 0.0001 4730 ELSE 4731 sinphi = sinp 4732 END IF 4733 4734 ! ratio between actual and sea-level pressure is used 4735 ! to correct for height in the computation of par ; 4736 ! should not exceed sea-level pressure therefore ... 4737 IF ( present(p) ) then 4738 pres = min( p, p_sealevel ) 4739 ELSE 4740 pres = p_sealevel 4741 ENDIF 4742 4743 ! Direct and diffuse par, Photoactive (=visible) radiation: 4744 CALL par_dir_diff( glrad, sinphi, pres, p_sealevel, pardir, pardiff ) 4745 4746 ! par for shaded leaves (canopy averaged): 4747 parshade = pardiff*exp(-0.5*LAI**0.7)+0.07*pardir*(1.1-0.1*LAI)*exp(-sinphi) ! Norman, 1982 4748 IF (glrad .gt. 200 .and. LAI .gt. 2.5) then 4749 parshade = pardiff*exp(-0.5*LAI**0.8)+0.07*pardir*(1.1-0.1*LAI)*exp(-sinphi) ! Zhang et al., 2001 4750 END IF 4751 4752 ! par for sunlit leaves (canopy averaged): 4753 ! alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5 4754 parsun = pardir*0.5/sinphi + parshade ! Norman, 1982 4755 IF (glrad .gt. 200 .and. LAI .gt. 2.5) then 4756 parsun = pardir**0.8*0.5/sinphi + parshade ! Zhang et al., 2001 4757 END IF 4758 4759 ! leaf area index for sunlit and shaded leaves: 4760 IF (sinphi .gt. 0) then 4761 LAIsun = 2*sinphi*(1-exp(-0.5*LAI/sinphi )) 4762 LAIshade = LAI -LAIsun 4763 ELSE 4764 LAIsun = 0 4765 LAIshade = LAI 4766 END IF 4767 4768 F_light = (LAIsun*(1 - exp(-1.*alpha(lu)*parsun)) + LAIshade*(1 - exp(-1.*alpha(lu)*parshade)))/LAI 4769 4770 F_light = max(F_light,F_min(lu)) 4771 4772 ! temperature influence; only non-zero within range [Tmin,Tmax]: 4773 IF ( (Tmin(lu) < T) .and. (T < Tmax(lu)) ) then 4774 BT = (Tmax(lu)-Topt(lu))/(Topt(lu)-Tmin(lu)) 4775 F_temp = ((T-Tmin(lu))/(Topt(lu)-Tmin(lu))) * ((Tmax(lu)-T)/(Tmax(lu)-Topt(lu)))**BT 4776 ELSE 4777 F_temp = 0.0 4778 END IF 4779 F_temp = max(F_temp,F_min(lu)) 4780 4781 ! vapour pressure deficit influence 4782 F_vpd = min(1.,((1.-F_min(lu))*(vpd_min(lu)-vpd)/(vpd_min(lu)-vpd_max(lu)) + F_min(lu) )) 4783 F_vpd = max(F_vpd,F_min(lu)) 4784 4785 F_swp = 1. 4786 4787 ! influence of phenology on stom. conductance 4788 ! ignored for now in DEPAC since influence of F_phen on lu classes in use is negligible. 4789 ! When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore 4790 F_phen = 1. 4791 4792 ! evaluate total stomatal conductance 4793 F_env = F_temp*F_vpd*F_swp 4794 F_env = max(F_env,F_min(lu)) 4795 gsto = G_max(lu) * F_light * F_phen * F_env 4796 4797 ! gstom expressed per m2 leafarea; 4798 ! this is converted with LAI to m2 surface. 4799 Gsto = lai*gsto ! in m/s 4800 4801 ELSE 4802 GSto = 0.0 4803 ENDIF 4804 4805 END SUBROUTINE rc_gstom_emb 4806 4807 4808 4809 !------------------------------------------------------------------- 4810 ! par_dir_diff 4811 !------------------------------------------------------------------- 4812 SUBROUTINE par_dir_diff(glrad,sinphi,pres,pres_0,par_dir,par_diff) 4813 ! 4814 ! Weiss, A., Norman, J.M. (1985) Partitioning solar radiation into direct and 4815 ! diffuse, visible and near-infrared components. Agric. Forest Meteorol. 4816 ! 34, 205-213. 4817 ! 4818 ! From a SUBROUTINE obtained from Leiming Zhang (via Willem Asman), 4819 ! MeteoroLOGICAL Service of Canada 4820 ! e-mail: leiming.zhang@ec.gc.ca 4821 ! 4822 ! Leiming uses solar irradiance. This should be equal to global radiation and 4823 ! Willem Asman set it to global radiation 4824 4825 REAL(kind=wp), INTENT(in) :: glrad ! global radiation (W m-2) 4826 REAL(kind=wp), INTENT(in) :: sinphi ! sine of the solar elevation 4827 REAL(kind=wp), INTENT(in) :: pres ! actual pressure (to correct for height) (Pa) 4828 REAL(kind=wp), INTENT(in) :: pres_0 ! pressure at sea level (Pa) 4829 REAL(kind=wp), INTENT(out) :: par_dir ! PAR direct : visible (photoactive) direct beam radiation (W m-2) 4830 REAL(kind=wp), INTENT(out) :: par_diff ! PAR diffuse: visible (photoactive) diffuse radiation (W m-2) 4831 4832 ! fn = near-infrared direct beam fraction (DIMENSIONless) 4833 ! fv = PAR direct beam fraction (DIMENSIONless) 4834 ! ratio = ratio measured to potential solar radiation (DIMENSIONless) 4835 ! rdm = potential direct beam near-infrared radiation (W m-2); "potential" means clear-sky 4836 ! rdn = potential diffuse near-infrared radiation (W m-2) 4837 ! rdu = visible (PAR) direct beam radiation (W m-2) 4838 ! rdv = potential visible (PAR) diffuse radiation (W m-2) 4839 ! rn = near-infrared radiation (W m-2) 4840 ! rv = visible radiation (W m-2) 4841 ! ww = water absorption in the near infrared for 10 mm of precipitable water 4842 4843 REAL(kind=wp) :: rdu,rdv,ww,rdm,rdn,rv,rn,ratio,sv,fv 4844 4845 ! Calculate visible (PAR) direct beam radiation 4846 ! 600 W m-2 represents average amount of PAR (400-700 nm wavelength) 4847 ! at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2) 4848 rdu=600.*exp(-0.185*(pres/pres_0)/sinphi)*sinphi 4849 4850 ! Calculate potential visible diffuse radiation 4851 rdv=0.4*(600.- rdu)*sinphi 4852 4853 ! Calculate the water absorption in the-near infrared 4854 ww=1320*10**( -1.195+0.4459*log10(1./sinphi)-0.0345*(log10(1./sinphi))**2 ) 4855 4856 ! Calculate potential direct beam near-infrared radiation 4857 rdm=(720.*exp(-0.06*(pres/pres_0)/sinphi)-ww)*sinphi !720 = solar constant - 600 4858 4859 ! Calculate potential diffuse near-infrared radiation 4860 rdn=0.6*(720-rdm-ww)*sinphi 4861 4862 ! Compute visible and near-infrared radiation 4863 rv=max(0.1,rdu+rdv) 4864 rn=max(0.01,rdm+rdn) 4865 4866 ! Compute ratio between input global radiation and total radiation computed here 4867 ratio=min(0.9,glrad/(rv+rn)) 4868 4869 ! Calculate total visible radiation 4870 sv=ratio*rv 4871 4872 ! Calculate fraction of PAR in the direct beam 4873 fv=min(0.99, (0.9-ratio)/0.7) ! help variable 4874 fv=max(0.01,rdu/rv*(1.0-fv**0.6667)) ! fraction of PAR in the direct beam 4875 4876 ! Compute direct and diffuse parts of PAR 4877 par_dir=fv*sv 4878 par_diff=sv-par_dir 4879 4880 END SUBROUTINE par_dir_diff 4881 4882 !------------------------------------------------------------------- 4883 ! rc_get_vpd: get vapour pressure deficit (kPa) 4884 !------------------------------------------------------------------- 4885 SUBROUTINE rc_get_vpd(temp, relh,vpd) 4886 4887 IMPLICIT NONE 4888 4889 ! Input/output variables: 4890 REAL(kind=wp) , INTENT(in) :: temp ! temperature (C) 4891 REAL(kind=wp) , INTENT(in) :: relh ! relative humidity (%) 4892 REAL(kind=wp) , INTENT(out) :: vpd ! vapour pressure deficit (kPa) 4893 4894 ! Local variables: 4895 REAL(kind=wp) :: esat 4896 4897 ! fit PARAMETERs: 4898 REAL(kind=wp), PARAMETER :: a1 = 6.113718e-1 4899 REAL(kind=wp), PARAMETER :: a2 = 4.43839e-2 4900 REAL(kind=wp), PARAMETER :: a3 = 1.39817e-3 4901 REAL(kind=wp), PARAMETER :: a4 = 2.9295e-5 4902 REAL(kind=wp), PARAMETER :: a5 = 2.16e-7 4903 REAL(kind=wp), PARAMETER :: a6 = 3.0e-9 4904 4905 ! esat is saturation vapour pressure (kPa) at temp(C) following monteith(1973) 4906 esat = a1 + a2*temp + a3*temp**2 + a4*temp**3 + a5*temp**4 + a6*temp**5 4907 vpd = esat * (1-relh/100) 4908 4909 END SUBROUTINE rc_get_vpd 4910 4911 4912 4913 !------------------------------------------------------------------- 4914 ! rc_gsoil_eff: compute effective soil conductance 4915 !------------------------------------------------------------------- 4916 SUBROUTINE rc_gsoil_eff(icmp,lu,sai,ust,nwet,t,gsoil_eff) 4917 4918 ! Input/output variables: 4919 INTEGER(iwp) , INTENT(in) :: icmp ! component index 4920 INTEGER(iwp) , INTENT(in) :: lu ! land use type, lu = 1,...,nlu 4921 REAL(kind=wp) , INTENT(in) :: sai ! surface area index 4922 REAL(kind=wp) , INTENT(in) :: ust ! friction velocity (m/s) 4923 INTEGER(iwp) , INTENT(in) :: nwet ! index for wetness 4924 ! nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow 4925 ! N.B. this routine cannot be CALLed with nwet = 9, 4926 ! nwet = 9 should be handled outside this routine. 4927 REAL(kind=wp) , INTENT(in) :: t ! temperature (C) 4928 REAL(kind=wp) , INTENT(out) :: gsoil_eff ! effective soil conductance (m/s) 4929 4930 ! local variables: 4931 REAL(kind=wp) :: rinc ! in canopy resistance (s/m) 4932 REAL(kind=wp) :: rsoil_eff ! effective soil resistance (s/m) 4933 4934 4935 4936 ! Soil resistance (numbers matched with lu_classes and component numbers) 4937 REAL(kind=wp), PARAMETER :: rsoil(nlu_dep,ncmp) = reshape( (/ & 4938 ! grs ara crp cnf dec wat urb oth des ice sav trf wai med sem 4939 1000., 200., 200., 200., 200., 2000., 400., 1000., 2000., 2000., 1000., 200., 2000., 200., 400., & ! O3 4940 1000., 1000., 1000., 1000., 1000., 10., 1000., 1000., 1000., 500., 1000., 1000., 10., 1000., 1000., & ! SO2 4941 1000., 1000., 1000., 1000., 1000., 2000., 1000., 1000., 1000., 2000., 1000., 1000., 2000., 1000., 1000., & ! NO2 4942 -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., & ! NO 4943 100., 100., 100., 100., 100., 10., 100., 100., 100., 1000., 100., 100., 10., 100., 100., & ! NH3 4944 -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., & ! CO 4945 -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! NO3 4946 -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! HNO3 4947 -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! N2O5 4948 -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999. /),& ! H2O2 4949 (/nlu_dep,ncmp/) ) 4950 4951 ! o3 so2 no2 no nh3 co no3 hno3 n2o5 h2o2 4952 4953 4954 REAL(kind=wp), PARAMETER :: rsoil_wet(ncmp) = (/2000., 10., 2000., -999., 10., -999., -999., -999., -999., -999./) 4955 REAL(kind=wp), PARAMETER :: rsoil_frozen(ncmp) = (/2000., 500., 2000., -999., 1000., -999., -999., -999., -999., -999./) 4956 4957 4958 4959 ! Compute in canopy (in crop) resistance: 4960 CALL rc_rinc(lu,sai,ust,rinc) 4961 4962 ! Check for missing deposition path: 4963 IF (missing(rinc)) then 4964 rsoil_eff = -9999. 4965 ELSE 4966 ! Frozen soil (temperature below 0): 4967 IF (t .lt. 0.0) then 4968 IF (missing(rsoil_frozen(icmp))) then 4969 rsoil_eff = -9999. 4970 ELSE 4971 rsoil_eff = rsoil_frozen(icmp) + rinc 4972 ENDIF 4973 ELSE 4974 ! Non-frozen soil; dry: 4975 IF (nwet .eq. 0) then 4976 IF (missing(rsoil(lu,icmp))) then 4977 rsoil_eff = -9999. 4978 ELSE 4979 rsoil_eff = rsoil(lu,icmp) + rinc 4980 ENDIF 4981 4982 ! Non-frozen soil; wet: 4983 ELSEIF (nwet .eq. 1) then 4984 IF (missing(rsoil_wet(icmp))) then 4985 rsoil_eff = -9999. 4986 ELSE 4987 rsoil_eff = rsoil_wet(icmp) + rinc 4988 ENDIF 4989 ELSE 4990 message_string = 'nwet can only be 0 or 1' 4991 CALL message( 'rc_gsoil_eff', 'CM0460', 1, 2, 0, 6, 0 ) 4992 ENDIF 4993 ENDIF 4994 ENDIF 4995 4996 ! Compute conductance: 4997 IF (rsoil_eff .gt. 0.0) then 4998 gsoil_eff = 1./rsoil_eff 4999 ELSE 5000 gsoil_eff = 0.0 5001 ENDIF 5002 5003 END SUBROUTINE rc_gsoil_eff 5004 5005 5006 5007 !------------------------------------------------------------------- 5008 ! rc_rinc: compute in canopy (or in crop) resistance 5009 ! van Pul and Jacobs, 1993, BLM 5010 !------------------------------------------------------------------- 5011 SUBROUTINE rc_rinc(lu,sai,ust,rinc) 5012 5013 5014 ! Input/output variables: 5015 INTEGER(iwp) , INTENT(in) :: lu ! land use class, lu = 1, ..., nlu 5016 REAL(kind=wp) , INTENT(in) :: sai ! surface area index 5017 REAL(kind=wp) , INTENT(in) :: ust ! friction velocity (m/s) 5018 REAL(kind=wp) , INTENT(out) :: rinc ! in canopy resistance (s/m) 5019 5020 5021 ! b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 IF not applicable) 5022 ! h = vegetation height (m) grass arabl crops conIF decid water urba othr desr ice sav trf wai med semi 5023 REAL(kind=wp), DIMENSION(nlu_dep), PARAMETER :: b = (/ -999, 14, 14, 14, 14, -999, -999, -999, -999, -999, -999, 14, -999, 14, 14 /) 5024 REAL(kind=wp), DIMENSION(nlu_dep), PARAMETER :: h = (/ -999, 1, 1, 20, 20, -999, -999, -999, -999, -999, -999, 20, -999, 1 , 1 /) 5025 5026 ! Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0: 5027 IF (b(lu) .gt. 0.0) then 5028 5029 ! Check for u* > 0 (otherwise denominator = 0): 5030 IF (ust .gt. 0.0) then 5031 rinc = b(lu)*h(lu)*sai/ust 5032 ELSE 5033 rinc = 1000.0 5034 ENDIF 5035 ELSE 5036 IF (lu .eq. ilu_grass .or. lu .eq. ilu_other ) then 5037 rinc = -999. ! no deposition path for grass, other, and semi-natural 5038 ELSE 5039 rinc = 0.0 ! no in-canopy resistance 5040 ENDIF 5041 ENDIF 5042 5043 END SUBROUTINE rc_rinc 5044 5045 5046 5047 !------------------------------------------------------------------- 5048 ! rc_rctot: compute total canopy (or surface) resistance Rc 5049 !------------------------------------------------------------------- 5050 SUBROUTINE rc_rctot(gstom,gsoil_eff,gw,gc_tot,rc_tot) 5051 5052 ! Input/output variables: 5053 REAL(kind=wp), INTENT(in) :: gstom ! stomatal conductance (s/m) 5054 REAL(kind=wp), INTENT(in) :: gsoil_eff ! effective soil conductance (s/m) 5055 REAL(kind=wp), INTENT(in) :: gw ! external leaf conductance (s/m) 5056 REAL(kind=wp), INTENT(out) :: gc_tot ! total canopy conductance (m/s) 5057 REAL(kind=wp), INTENT(out) :: rc_tot ! total canopy resistance Rc (s/m) 5058 5059 ! Total conductance: 5060 gc_tot = gstom + gsoil_eff + gw 5061 5062 ! Total resistance (note: gw can be negative, but no total emission allowed here): 5063 IF (gc_tot .le. 0.0 .or. gw .lt. 0.0) then 5064 rc_tot = -9999. 5065 ELSE 5066 rc_tot = 1./gc_tot 5067 ENDIF 5068 5069 END SUBROUTINE rc_rctot 5070 5071 5072 5073 !------------------------------------------------------------------- 5074 ! rc_comp_point_rc_eff: calculate the effective resistance Rc 5075 ! based on one or more compensation points 5076 !------------------------------------------------------------------- 5077 ! 5078 ! NH3rc (see depac v3.6 is based on Avero workshop Marc Sutton. p. 173. 5079 ! Sutton 1998 AE 473-480) 5080 ! 5081 ! Documentation by Ferd Sauter, 2008; see also documentation block in header of depac subroutine. 5082 ! FS 2009-01-29: variable names made consistent with DEPAC 5083 ! FS 2009-03-04: use total compensation point 5084 ! 5085 ! C: with total compensation point ! D: approximation of C 5086 ! ! with classical approach 5087 ! zr --------- Catm ! zr --------- Catm 5088 ! | ! | 5089 ! Ra ! Ra 5090 ! | ! | 5091 ! Rb ! Rb 5092 ! | ! | 5093 ! z0 --------- Cc ! z0 --------- Cc 5094 ! | ! | 5095 ! Rc ! Rc_eff 5096 ! | ! | 5097 ! --------- Ccomp_tot ! --------- C=0 5098 ! 5099 ! 5100 ! The effective Rc is defined such that instead of using 5101 ! 5102 ! F = -vd*[Catm - Ccomp_tot] (1) 5103 ! 5104 ! we can use the 'normal' flux formula 5105 ! 5106 ! F = -vd'*Catm, (2) 5107 ! 5108 ! with vd' = 1/(Ra + Rb + Rc') (3) 5109 ! 5110 ! and Rc' the effective Rc (rc_eff). 5111 ! (Catm - Ccomp_tot) 5112 ! vd'*Catm = vd*(Catm - Ccomp_tot) <=> vd' = vd* ------------------ 5113 ! Catm 5114 ! 5115 ! (Catm - Ccomp_tot) 5116 ! 1/(Ra + Rb + Rc') = (1/Ra + Rb + Rc) * ------------------ 5117 ! Catm 5118 ! 5119 ! Catm 5120 ! (Ra + Rb + Rc') = (Ra + Rb + Rc) * ------------------ 5121 ! (Catm - Ccomp_tot) 5122 ! 5123 ! Catm 5124 ! Rc' = (Ra + Rb + Rc) * ------------------ - Ra - Rb 5125 ! (Catm - Ccomp_tot) 5126 ! 5127 ! Catm Catm 5128 ! Rc' = (Ra + Rb) [------------------ - 1 ] + Rc * ------------------ 5129 ! (Catm - Ccomp_tot) (Catm - Ccomp_tot) 5130 ! 5131 ! Rc' = [(Ra + Rb)*Ccomp_tot + Rc*Catm ] / (Catm - Ccomp_tot) 5132 ! 5133 ! This is not the most efficient way to DO this; 5134 ! in the current LE version, (1) is used directly in the CALLing routine 5135 ! and this routine is not used anymore. 5136 ! 5137 ! ------------------------------------------------------------------------------------------- 5138 ! In fact it is the question IF this correct; note the difference in differential equations 5139 ! 5140 ! dCatm ! dCatm 5141 ! H ----- = F = -vd'*Catm, ! H ----- = F = -vd*[Catm - Ccomp_tot], 5142 ! dt ! dt 5143 ! 5144 ! where in the left colum vd' is a function of Catm and not a constant anymore. 5145 ! 5146 ! Another problem is that IF Catm -> 0, we end up with an infinite value of the exchange 5147 ! velocity vd. 5148 ! ------------------------------------------------------------------------------------------- 5149 5150 SUBROUTINE rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff) 5151 5152 IMPLICIT NONE 5153 5154 ! Input/output variables: 5155 REAL(kind=wp), INTENT(in) :: ccomp_tot ! total compensation point (weighed average of separate compensation points) (ug/m3) 5156 REAL(kind=wp), INTENT(in) :: catm ! atmospheric concentration (ug/m3) 5157 REAL(kind=wp), INTENT(in) :: ra ! aerodynamic resistance (s/m) 5158 REAL(kind=wp), INTENT(in) :: rb ! boundary layer resistance (s/m) 5159 REAL(kind=wp), INTENT(in) :: rc_tot ! total canopy resistance (s/m) 5160 REAL(kind=wp), INTENT(out) :: rc_eff ! effective total canopy resistance (s/m) 5161 5162 ! Compute effective resistance: 5163 IF ( ccomp_tot == 0.0 ) then 5164 ! trace with no compensiation point ( or compensation point equal to zero) 5165 rc_eff = rc_tot 5166 5167 ELSE IF ( ccomp_tot > 0 .and. ( abs( catm - ccomp_tot ) < 1.e-8 ) ) then 5168 ! surface concentration (almoast) equal to atmospheric concentration 5169 ! no exchange between surface and atmosphere, infinite RC --> vd=0 5170 rc_eff = 9999999999. 5171 5172 ELSE IF ( ccomp_tot > 0 ) then 5173 ! compensation point available, calculate effective restisance 5174 rc_eff = ((ra + rb)*ccomp_tot + rc_tot*catm)/(catm-ccomp_tot) 5175 5176 ELSE 5177 rc_eff = -999. 5178 message_string = 'This should not be possible, check ccomp_tot' 5179 CALL message( 'rc_comp_point_rc_eff', 'CM0461', 1, 2, 0, 6, 0 ) 5180 ENDIF 5181 5182 return 5183 END SUBROUTINE rc_comp_point_rc_eff 5184 5185 5186 5187 !------------------------------------------------------------------- 5188 ! missing: check for data that correspond with a missing deposition path 5189 ! this data is represented by -999 5190 !------------------------------------------------------------------- 5191 5192 LOGICAL function missing(x) 5193 5194 REAL(kind=wp), INTENT(in) :: x 5195 5196 ! bandwidth for checking (in)equalities of floats 5197 REAL(kind=wp), PARAMETER :: EPS = 1.0e-5 5198 5199 missing = (abs(x + 999.) .le. EPS) 5200 5201 END function missing 5202 5203 5204 5205 ELEMENTAL FUNCTION sedimentation_velocity(rhopart, partsize, slipcor, visc) RESULT (vs) 5206 5207 5208 IMPLICIT NONE 5209 5210 ! --- in/out --------------------------------- 5211 5212 REAL(kind=wp), INTENT(in) :: rhopart ! kg/m3 5213 REAL(kind=wp), INTENT(in) :: partsize ! m 5214 REAL(kind=wp), INTENT(in) :: slipcor ! m 5215 REAL(kind=wp), INTENT(in) :: visc ! viscosity 5216 REAL(kind=wp) :: vs 5217 5218 ! acceleration of gravity: 5219 REAL(kind=wp), PARAMETER :: grav = 9.80665 ! m/s2 5220 5221 ! --- begin ---------------------------------- 5222 5223 !viscosity & sedimentation velocity 5224 vs = rhopart * (partsize**2) * grav * slipcor / (18*visc) 5225 5226 END FUNCTION sedimentation_velocity 5227 5228 !------------------------------------------------------------------------ 5229 ! Boundary-layer deposition resistance following Zhang (2001) 5230 !------------------------------------------------------------------------ 5231 SUBROUTINE drydepo_aero_zhang_vd( vd, Rs, vs1, partsize, slipcor, & 5232 nwet, tsurf, dens1, viscos1, & 5233 luc, ftop_lu, ustar_lu) 5234 5235 5236 IMPLICIT NONE 5237 5238 ! --- in/out --------------------------------- 5239 5240 REAL(kind=wp), INTENT(out) :: vd ! deposition velocity (m/s) 5241 REAL(kind=wp), INTENT(out) :: Rs ! Sedimentaion resistance (s/m) 5242 REAL(kind=wp), INTENT(in) :: vs1 ! sedimentation velocity in lowest layer 5243 REAL(kind=wp), INTENT(in) :: partsize ! particle diameter (m) 5244 REAL(kind=wp), INTENT(in) :: slipcor ! slip correction factor 5245 REAL(kind=wp), INTENT(in) :: tsurf ! surface temperature (K) 5246 REAL(kind=wp), INTENT(in) :: dens1 ! air density (kg/m3) in lowest layer 5247 REAL(kind=wp), INTENT(in) :: viscos1 ! air viscosity in lowest layer 5248 REAL(kind=wp), INTENT(in) :: ftop_lu ! atmospheric resistnace Ra 5249 REAL(kind=wp), INTENT(in) :: ustar_lu ! friction velocity u* 5250 5251 INTEGER(iwp), INTENT(in) :: nwet ! 1=rain, 9=snowcover 5252 INTEGER(iwp), INTENT(in) :: luc ! DEPAC LU 5253 5254 5255 ! --- const ---------------------------------- 5256 5257 5258 5259 ! acceleration of gravity: 5260 REAL(kind=wp), PARAMETER :: grav = 9.80665 ! m/s2 5261 5262 REAL(kind=wp), PARAMETER :: beta = 2.0 5263 REAL(kind=wp), PARAMETER :: epsilon0 = 3.0 5264 REAL(kind=wp), PARAMETER :: kb = 1.38066e-23 5265 REAL(kind=wp), PARAMETER :: pi = 3.141592654_wp !< PI 5266 5267 REAL, PARAMETER :: alfa_lu(nlu_dep) = & 5268 (/1.2, 1.2, 1.2, 1.0, 1.0, 100.0, 1.5, 1.2, 50.0, 100.0, 1.2, 1.0, 100.0, 1.2, 50.0/) 5269 ! grass arabl crops conif decid water urba othr desr ice sav trf wai med sem 5270 REAL, PARAMETER :: gamma_lu(nlu_dep) = & 5271 (/0.54, 0.54, 0.54, 0.56, 0.56, 0.50, 0.56, 0.54, 0.58, 0.50, 0.54, 0.56, 0.50, 0.54, 0.54/) 5272 ! grass arabl crops conif decid water urba othr desr ice sav trf wai med sem 5273 REAL, PARAMETER ::A_lu(nlu_dep) = & 5274 (/3.0, 3.0, 2.0, 2.0, 7.0, -99., 10.0, 3.0, -99., -99., 3.0, 7.0, -99., 2.0, -99./) 5275 ! grass arabl crops conif decid water urba othr desr ice sav trf wai med sem 5276 5277 5278 ! --- local ---------------------------------- 5279 5280 REAL(kind=wp) :: kinvisc, diff_part 5281 REAL(kind=wp) :: schmidt,stokes, Ebrown, Eimpac, Einterc, Reffic 5282 5283 ! --- begin ---------------------------------- 5284 5285 ! kinetic viscosity & diffusivity 5286 kinvisc = viscos1 / dens1 ! only needed at surface 5287 5288 diff_part = kb * tsurf * slipcor / (3*pi*viscos1*partsize) 5289 5290 ! Schmidt number 5291 schmidt = kinvisc / diff_part 5292 5293 ! calculate collection efficiencie E 5294 Ebrown = Schmidt**(-gamma_lu(luc)) ! Brownian diffusion 5295 ! determine Stokes number, interception efficiency 5296 ! and sticking efficiency R (1 = no rebound) 5297 IF ( luc == ilu_ice .or. nwet.eq.9 .or. luc == ilu_water_sea .or. luc == ilu_water_inland ) THEN 5298 stokes=vs1*ustar_lu**2/(grav*kinvisc) 5299 Einterc=0.0 5300 Reffic=1.0 5301 ELSE IF ( luc == ilu_other .or. luc == ilu_desert ) THEN !tundra of desert 5302 stokes=vs1*ustar_lu**2/(grav*kinvisc) 5303 Einterc=0.0 5304 Reffic=exp(-Stokes**0.5) 5305 ELSE 5306 stokes=vs1*ustar_lu/(grav*A_lu(luc)*1.e-3) 5307 Einterc=0.5*(partsize/(A_lu(luc)*1e-3))**2 5308 Reffic=exp(-Stokes**0.5) 5309 END IF 5310 ! when surface is wet all particles DO not rebound: 5311 IF(nwet.eq.1) Reffic = 1.0 5312 ! determine impaction efficiency: 5313 Eimpac = ( stokes / (alfa_lu(luc)+stokes) )**beta 5314 5315 ! sedimentation resistance: 5316 Rs = 1.0 / ( epsilon0 * ustar_lu * (Ebrown+Eimpac+Einterc) * Reffic ) 5317 5318 ! deposition velocity according to Seinfeld and Pandis (= SP, 2006; eq 19.7): 5319 ! 5320 ! 1 5321 ! vd = ------------------ + vs 5322 ! Ra + Rs + Ra*Rs*vs 5323 ! 5324 ! where: Rs = Rb (in SP, 2006) 5325 ! 5326 ! 5327 vd = 1.0 / ( ftop_lu + Rs + ftop_lu*Rs*vs1) + vs1 5328 5329 5330 5331 END SUBROUTINE drydepo_aero_zhang_vd 5332 5333 5334 5335 SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffc, Rb ) 5336 5337 ! Compute quasi-laminar boundary layer resistance as a function of landuse and tracer. 5338 ! 5339 ! Original EMEP formulation by (Simpson et al, 2003) is used. 5340 ! 5341 5342 5343 IMPLICIT NONE 5344 5345 ! --- in/out --------------------------------- 5346 5347 LOGICAL , INTENT(in) :: is_water 5348 REAL(kind=wp), INTENT(in) :: z0h ! roughness length for heat 5349 REAL(kind=wp), INTENT(in) :: ustar ! friction velocity 5350 REAL(kind=wp), INTENT(in) :: diffc ! coefficient of diffusivity 5351 REAL(kind=wp), INTENT(out) :: Rb ! boundary layer resistance 5352 5353 ! --- const ---------------------------------- 5354 5355 ! thermal diffusivity of dry air 20 C 5356 REAL(kind=wp), PARAMETER :: thk = 0.19e-4 ! http://en.wikipedia.org/wiki/Thermal_diffusivity (@ T=300K) 5357 REAL(kind=wp), PARAMETER :: kappa_stab = 0.35 ! von Karman constant 5358 5359 ! --- begin --------------------------------- 5360 5361 5362 ! Use Simpson et al. (2003) 5363 !IF ( is_water ) THEN 5364 !Rb = 1.0_wp / (kappa_stab*max(0.01_wp,ustar)) * log(z0h/diffc*kappa_stab*max(0.01_wp,ustar)) 5365 !TODO: Check Rb over water calculation!!! 5366 ! Rb = 1.0_wp / (kappa_stab*max(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*max(0.1_wp,ustar)) 5367 !ELSE 5368 Rb = 5.0_wp / max(0.01_wp,ustar) * (thk/diffc)**0.67_wp 5369 !END IF 5370 5371 END SUBROUTINE get_rb_cell 5372 5373 5374 !----------------------------------------------------------------- 5375 ! Compute water vapor partial pressure (e_w) 5376 ! given specific humidity Q [(kg water)/(kg air)]. 5377 ! 5378 ! Use that gas law for volume V with temperature T 5379 ! holds for the total mixture as well as the water part: 5380 ! 5381 ! R T / V = p_air / n_air = p_water / n_water 5382 ! 5383 ! thus: 5384 ! 5385 ! p_water = p_air n_water / n_air 5386 ! 5387 ! Use: 5388 ! n_air = m_air / xm_air 5389 ! [kg air] / [(kg air)/(mole air)] 5390 ! and: 5391 ! n_water = m_air * Q / xm_water 5392 ! [kg water] / [(kg water)/(mole water)] 5393 ! thus: 5394 ! p_water = p_air Q / (xm_water/xm_air) 5395 ! 5396 5397 ELEMENTAL FUNCTION watervaporpartialpressure( Q, p ) result ( p_w ) 5398 5399 5400 ! --- in/out --------------------------------- 5401 5402 REAL(kind=wp), INTENT(in) :: Q ! specific humidity [(kg water)/(kg air)] 5403 REAL(kind=wp), INTENT(in) :: p ! air pressure [Pa] 5404 REAL(kind=wp) :: p_w ! water vapor partial pressure [Pa] 5405 5406 ! --- const ---------------------------------- 5407 5408 ! mole mass ratio: 5409 REAL(kind=wp), PARAMETER :: eps = xm_H2O / xm_air ! ~ 0.622 5410 5411 ! --- begin ---------------------------------- 5412 5413 ! partial pressure of water vapor: 5414 p_w = p * Q / eps 5415 5416 END function watervaporpartialpressure 5417 5418 5419 5420 !------------------------------------------------------------------ 5421 ! Saturation vapor pressure. 5422 ! From (Stull 1988, eq. 7.5.2d): 5423 ! 5424 ! e_sat = p0 exp( 17.67 * (T-273.16) / (T-29.66) ) [Pa] 5425 ! 5426 ! where: 5427 ! p0 = 611.2 [Pa] : reference pressure 5428 ! 5429 ! Arguments: 5430 ! T [K] : air temperature 5431 ! Result: 5432 ! e_sat_w [Pa] : saturation vapor pressure 5433 ! 5434 ! References: 5435 ! Roland B. Stull, 1988 5436 ! An introduction to boundary layer meteorology. 5437 ! 5438 5439 ELEMENTAL FUNCTION saturationvaporpressure( T ) result( e_sat_w ) 5440 5441 5442 ! --- in/out --------------------------------- 5443 5444 REAL(kind=wp), INTENT(in) :: T ! temperature [K] 5445 REAL(kind=wp) :: e_sat_w ! saturation vapor pressure [Pa] 5446 5447 ! --- const ---------------------------------- 5448 5449 ! base pressure: 5450 REAL(kind=wp), PARAMETER :: p0 = 611.2 ! [Pa] 5451 5452 ! --- begin ---------------------------------- 5453 5454 ! saturation vapor pressure: 5455 e_sat_w = p0 * exp( 17.67 * (T-273.16) / (T-29.66) ) ! [Pa] 5456 5457 END FUNCTION saturationvaporpressure 5458 5459 5460 5461 !------------------------------------------------------------------------ 5462 ! Relative humidity RH [%] is by definition: 5463 ! 5464 ! e_w # water vapor partial pressure 5465 ! Rh = -------- * 100 5466 ! e_sat_w # saturation vapor pressure 5467 ! 5468 ! Compute here from: 5469 ! Q specific humidity [(kg water)/(kg air)] 5470 ! T temperature [K] 5471 ! P pressure [Pa] 5472 ! 5473 5474 ELEMENTAL FUNCTION relativehumidity_from_specifichumidity( Q, T, p ) result( Rh ) 5475 5476 5477 ! --- in/out --------------------------------- 5478 5479 REAL(kind=wp), INTENT(in) :: Q ! specific humidity [(kg water)/(kg air)] 5480 REAL(kind=wp), INTENT(in) :: T ! temperature [K] 5481 REAL(kind=wp), INTENT(in) :: p ! air pressure [Pa] 5482 REAL(kind=wp) :: Rh ! relative humidity [%] 5483 5484 ! --- begin ---------------------------------- 5485 5486 ! relative humidity: 5487 Rh = watervaporpartialpressure( Q, p ) / saturationvaporpressure( T ) * 100.0 5488 5489 END FUNCTION relativehumidity_from_specifichumidity 5490 5491 5492 2434 5493 END MODULE chemistry_model_mod 2435 5494 -
palm/trunk/SOURCE/date_and_time_mod.f90
r3298 r3458 19 19 ! 20 20 ! Current revisions: 21 ! ----------------- 21 ! ------------------ 22 22 ! 23 23 ! … … 25 25 ! ----------------- 26 26 ! $Id$ 27 ! from chemistry branch r3443, banzhafs: 28 ! Added initial hour_of_day, hour_of_year, day_of_year and month_of_year to 29 ! init_date_and_time 30 ! 31 ! 3298 2018-10-02 12:21:11Z kanani 27 32 ! - Minor formatting (kanani) 28 33 ! - Added Routines for DEFAULT mode of chemistry emissions (Russo) … … 61 66 !> already implemented changes for calculating it from date_init in 62 67 !> calc_date_and_time 68 !> @todo time_utc during spin-up 63 69 !------------------------------------------------------------------------------! 64 70 MODULE date_and_time_mod … … 91 97 INTEGER(iwp) :: index_mm !< index months of the default emission mode 92 98 INTEGER(iwp) :: index_dd !< index days of the default emission mode 93 INTEGER(iwp) :: index_hh !< index hours of the default emission mode 99 INTEGER(iwp) :: index_hh !< index hours of the emission mode 100 INTEGER(iwp) :: hours_since_reference_point !< hours of current simulation 94 101 95 102 REAL(wp) :: time_utc !< current model time in UTC … … 102 109 REAL(wp), PARAMETER :: d_seconds_year = 1.0_wp / 31536000.0_wp !< inverse of the seconds per year (1/(365*86400)) 103 110 104 CHARACTER(len=8) :: date_init = " 31122018" !< Starting date of simulation: We selected this because it was a monday111 CHARACTER(len=8) :: date_init = "21062017" !< Starting date of simulation: We selected this because it was a monday 105 112 106 113 !> --- Parameters … … 121 128 END INTERFACE time_default_indices 122 129 130 !-- Get hour index in the PRE-PROCESSED case of chemistry emissions : 131 INTERFACE time_preprocessed_indices 132 MODULE PROCEDURE time_preprocessed_indices 133 END INTERFACE time_preprocessed_indices 134 135 123 136 !-- Calculate current date and time 124 137 INTERFACE calc_date_and_time … … 128 141 129 142 !-- Public Interfaces 130 PUBLIC calc_date_and_time, time_default_indices, init_date_and_time 143 PUBLIC calc_date_and_time, time_default_indices, init_date_and_time, time_preprocessed_indices 131 144 132 145 !-- Public Variables … … 147 160 IMPLICIT NONE 148 161 162 !-- Variables Definition 163 INTEGER :: i_mon !< Index for going through the different months 164 149 165 IF (day_of_year_init == 0) THEN 150 166 ! Day of the month at starting time … … 159 175 ENDIF 160 176 161 END SUBROUTINE init_date_and_time 162 163 !------------------------------------------------------------------------------! 164 ! Description: 165 ! ------------ 166 !> Calculate current date and time of the simulation 167 !------------------------------------------------------------------------------! 168 169 SUBROUTINE calc_date_and_time 170 171 IMPLICIT NONE 172 173 !-- Variables Definition 174 INTEGER :: i_mon !< Index for going through the different months 175 176 !> Update simulation time in seconds 177 time_update = simulated_time-coupling_start_time 178 179 !-- Calculate current day of the simulated time 180 days_since_reference_point=INT(FLOOR( (time_utc_init + time_update) & 181 / 86400.0_wp ) ) 182 183 !-- Calculate actual UTC time 184 time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp) 185 186 !sB PRILIMINARY workaround for time_utc bug concerning time_since_reference_point: 187 time_utc_emis = MOD((time_utc_init + time_update), 86400.0_wp) 188 189 !-- Calculate initial day of the year: it is calculated only once. In fact, day_of_year_init is initialized to 0 and then a positive value is passed. This condition is also called only when day_of_year_init is not given in the namelist. 190 177 178 !-- Calculate initial hour of the day: the first hour of the day is from 00:00:00 to 00:59:59. 179 180 hour_of_day = INT( FLOOR( time_utc_init/3600.0_wp ) ) + 1 181 182 !-- Calculate initial day day_of_year_init in case date_init is given or day_of_year_init is given 191 183 IF ( day_of_year_init == 0 ) THEN 192 184 … … 214 206 ENDIF 215 207 208 209 !-- Initial day of the year 210 day_of_year = day_of_year_init 211 212 !-- Initial hour of the year 213 hour_of_year = ( (day_of_year-1) * 24 ) + hour_of_day 214 215 !--Initial day of the month and month of the year 216 !> -------------------------------------------------------------------------------- 217 !> The first case is when date_init is not provided: we only know day_of_year_init 218 IF ( month_of_year == 0 .AND. day_of_month == 0) THEN 219 220 221 IF ( day_of_year .LE. 31 ) THEN 222 223 month_of_year=1 224 day_of_month=day_of_year 225 226 ELSE 227 228 DO i_mon=2,12 !january is considered in the first case 229 IF ( day_of_year .LE. SUM(days(1:i_mon)) .AND. day_of_year .GT. SUM(days(1:(i_mon-1))) ) THEN 230 231 month_of_year=i_mon 232 233 day_of_month=INT(MOD(day_of_year, SUM(days(1:(i_mon-1))))) 234 235 GOTO 38 236 237 ENDIF 238 239 38 ENDDO 240 ENDIF 241 !> -------------------------------------------------------------------------------- 242 !> in the second condition both day of month and month_of_year are either given in input (passed to date_init) or we are in some day successive to the initial one, so that day_of_month has already be computed in previous step 243 !>TBD: something to calculate the current year is missing 244 ELSEIF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN 245 246 !> calculate month_of_year. TBD: test the condition when day_of_year==31 247 248 IF (day_of_year==1) THEN !> this allows to turn from december to January when passing from a year to another 249 250 month_of_year = 1 251 252 ELSE IF (day_of_year .GT. 1 .AND. day_of_year .GT. SUM(days(1:month_of_year))) THEN 253 254 month_of_year = month_of_year + 1 255 256 ENDIF 257 258 !> calculate day_of_month 259 IF ( month_of_year == 1 ) THEN 260 261 day_of_month=day_of_year 262 263 ELSE 264 265 day_of_month=INT(MOD(day_of_year, SUM(days(1:(month_of_year-1))))) 266 267 ENDIF 268 269 270 ELSE 271 272 !> Condition when date_init is provided but it is given in the wrong format 273 message_string = 'date_init not provided in the namelist or' // & 274 ' given in the wrong format: MUST BE DDMMYYYY' 275 CALL message( 'init_date_and_time', 'DT0102', 2, 2, 0, 6, 0 ) 276 277 ENDIF 278 279 280 END SUBROUTINE init_date_and_time 281 282 !------------------------------------------------------------------------------! 283 ! Description: 284 ! ------------ 285 !> Calculate current date and time of the simulation 286 !------------------------------------------------------------------------------! 287 288 SUBROUTINE calc_date_and_time 289 290 IMPLICIT NONE 291 292 !-- Variables Definition 293 INTEGER :: i_mon !< Index for going through the different months 294 295 !> Update simulation time in seconds 296 time_update = simulated_time-coupling_start_time 297 298 !-- Calculate current day of the simulated time 299 days_since_reference_point=INT(FLOOR( (time_utc_init + time_update) & 300 / 86400.0_wp ) ) 301 302 !-- Calculate actual UTC time 303 time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp) 304 305 !sB PRILIMINARY workaround for time_utc changes due to changes in time_since_reference_point in 306 !sB radiation_model_mod during runtime: 307 time_utc_emis = MOD((time_utc_init + time_update), 86400.0_wp) 308 309 !-- Calculate initial day of the year: it is calculated only once. In fact, day_of_year_init is initialized to 0 and then a positive value is passed. This condition is also called only when day_of_year_init is not given in the namelist. 310 311 IF ( day_of_year_init == 0 ) THEN 312 313 !> Condition for printing an error when date_init is not provided when day_of_year_init is not given in the namelist or when the format of the date is not the one required by PALM. 314 IF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN 315 316 IF ( month_of_year == 1 ) THEN !!month of year is read in input 317 318 day_of_year_init = day_of_month 319 320 ELSE 321 322 day_of_year_init= SUM(days( 1:(month_of_year-1) )) + day_of_month !day_of_month is read in input in this case 323 324 ENDIF 325 !kanani: Revise, we cannot force users to provide date_init, maybe set a default value? 326 ! ELSE 327 ! 328 ! message_string = 'date_init not provided in the namelist or' // & 329 ! ' given in the wrong format: MUST BE DDMMYYYY' 330 ! CALL message( 'calc_date_and_time', 'DT0100', 2, 2, 0, 6, 0 ) 331 332 ENDIF 333 334 ENDIF 335 216 336 !-- Calculate actual hour of the day: the first hour of the day is from 00:00:00 to 00:59:59. 217 337 … … 302 422 303 423 ENDIF 304 424 305 425 END SUBROUTINE calc_date_and_time 306 426 307 !TBD: check the routines used for update of emission indices in the DEFAULT MODE 427 428 !------------------------------------------------------------------------------! 429 ! Description: 430 ! ------------ 431 !> This routine determines the time factor index in the PRE-PROCESSED emissions mode. 432 !------------------------------------------------------------------------------! 433 434 SUBROUTINE time_preprocessed_indices(index_hh) 435 436 USE indices 437 438 IMPLICIT NONE 439 440 ! 441 !-- In/output 442 INTEGER, INTENT(INOUT) :: index_hh !> Index Hour 443 ! 444 !-- Additional Variables for calculateing indices 445 !-- Constants 446 INTEGER, PARAMETER :: nhour = 24 447 448 IF (days_since_reference_point == 0) THEN 449 450 index_hh=hour_of_day 451 452 ELSE 453 454 index_hh=(days_since_reference_point*nhour)+(hour_of_day) 455 456 ENDIF 457 458 459 END SUBROUTINE time_preprocessed_indices 460 461 308 462 !------------------------------------------------------------------------------! 309 463 ! Description: -
palm/trunk/SOURCE/flow_statistics.f90
r3452 r3458 25 25 ! ----------------- 26 26 ! $Id$ 27 ! from chemistry branch r3443, basit: 28 ! bug fixed in chemistry profiles 29 ! 30 ! 3452 2018-10-30 13:13:34Z schwenkel 27 31 ! Bugfix for profiles output 28 32 ! … … 1843 1847 ENDIF 1844 1848 1845 IF ( air_chemistry .AND.max_pr_cs > 0 ) THEN1849 IF ( air_chemistry .AND. max_pr_cs > 0 ) THEN 1846 1850 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1847 CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+1,0), sums(nzb,pr_palm+1), & 1848 nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr ) 1851 DO i = 1, max_pr_cs 1852 CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+i,0), & 1853 sums(nzb,pr_palm+max_pr_user+i), & 1854 nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr ) 1855 ENDDO 1849 1856 ENDIF 1850 1857 -
palm/trunk/SOURCE/init_3d_model.f90
r3448 r3458 25 25 ! ----------------- 26 26 ! $Id$ 27 ! from chemistry branch r3443, basit: 28 ! bug fixed in sums and sums_l for chemistry profile output 29 ! 30 ! 3448 2018-10-29 18:14:31Z kanani 27 31 ! Add biometeorology 28 32 ! … … 711 715 ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions), & 712 716 rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions), & 713 sums(nzb:nzt+1,pr_palm+max_pr_user ),&714 sums_l(nzb:nzt+1,pr_palm+max_pr_user ,0:threads_per_task-1), &717 sums(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs), & 718 sums_l(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs,0:threads_per_task-1), & 715 719 sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), & 716 720 sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions), & -
palm/trunk/SOURCE/netcdf_data_input_mod.f90
r3429 r3458 25 25 ! ----------------- 26 26 ! $Id$ 27 ! from chemistry branch r3443, banzhafs, Russo: 28 ! Uncommented lines on dimension of surface_fractions 29 ! Removed par_emis_time_factor variable, moved to chem_emissions_mod 30 ! Initialized nspec and other emission variables at time of declaration 31 ! Modified EXPERT mode to PRE-PROCESSED mode 32 ! Introduced Chemistry static netcdf file 33 ! Added the routine for reading-in netcdf data for chemistry 34 ! Added routines to get_variable interface specific for chemistry files 35 ! 36 ! 3429 2018-10-25 13:04:23Z knoop 27 37 ! add default values of origin_x/y/z 28 38 ! … … 368 378 369 379 !-DIMENSIONS 370 INTEGER(iwp) :: nspec 371 INTEGER(iwp) :: ncat 372 INTEGER(iwp) :: nvoc 373 INTEGER(iwp) :: npm 380 INTEGER(iwp) :: nspec=0 !< number of chem species for which emission values are provided 381 INTEGER(iwp) :: ncat=0 !< number of emission categories 382 INTEGER(iwp) :: nvoc=0 !< number of VOCs components 383 INTEGER(iwp) :: npm=0 !< number of PMs components 374 384 INTEGER(iwp) :: nnox=2 !< number of NOx components: NO and NO2 375 385 INTEGER(iwp) :: nsox=2 !< number of SOx components: SO and SO4 … … 391 401 INTEGER(iwp),ALLOCATABLE, DIMENSION(:) :: species_index !< Index of emission chem species 392 402 393 REAL(wp), DIMENSION(24) :: par_emis_time_factor !< time factors394 ! for the parameterized mode: these are fixed for each hour395 ! of a single day.396 403 REAL(wp),ALLOCATABLE, DIMENSION(:) :: xm !< Molecular masses of emission chem species 397 404 … … 928 935 ENDDO 929 936 930 !Assign Constant Values of time factors:931 emt_att%par_emis_time_factor( : ) = (/ 0.01, 0.01, 0.01, 0.02, 0.03, 0.07, 0.09, 0.09, 0.05, 0.03, 0.03, 0.03, &932 0.03, 0.03, 0.03, 0.04, 0.05, 0.09, 0.09, 0.09, 0.04, 0.02, 0.01, 0.01 /)933 937 934 938 !> DEFAULT AND PRE-PROCESSED MODE -
palm/trunk/SOURCE/palm.f90
r3337 r3458 25 25 ! ----------------- 26 26 ! $Id$ 27 ! from chemistry branch r3443, forkel: 28 ! removed double do_emis check around CALL chem_init 29 ! replaced call to calc_date_and_time to init_date_and_time 30 ! 31 ! 3337 2018-10-12 15:17:09Z kanani 27 32 ! (from branch resler) 28 33 ! Fix chemistry call … … 436 441 437 442 CALL chem_init 443 438 444 ! CALL photolysis_init ! probably also required for restart 439 440 CALL init_date_and_time !initialize the time of chemistry emissions441 445 442 446 ENDIF … … 448 452 ! 449 453 !-- Initialize all necessary variables 450 CALL calc_date_and_time !this is required for chemistry emissions 454 ! 455 !-- Initial time for chem_emissions_mod 456 CALL init_date_and_time 451 457 452 458 CALL init_3d_model -
palm/trunk/SOURCE/prognostic_equations.f90
r3386 r3458 25 25 ! ----------------- 26 26 ! $Id$ 27 ! remove duplicate USE chem_modules 28 ! from chemistry branch r3443, banzhafs, basit: 29 ! chem_depo call introduced 30 ! code added for decycling chemistry 31 ! 32 ! 3386 2018-10-19 16:28:22Z gronemeier 27 33 ! Renamed tcm_prognostic to tcm_prognostic_equations 28 34 ! … … 340 346 USE buoyancy_mod, & 341 347 ONLY: buoyancy 342 343 USE chemistry_model_mod, &344 ONLY: chem_integrate, chem_species, chem_prognostic_equations, &345 nspec, nvar, spc_names, chem_boundary_conds346 348 347 349 USE chem_modules, & 348 ONLY: call_chem_at_all_substeps, chem_gasphase_on 350 ONLY: call_chem_at_all_substeps, chem_gasphase_on, cs_name, do_depo 349 351 350 352 USE chem_photolysis_mod, & 351 353 ONLY: photolysis_control 352 354 353 USE chem_modules, & 354 ONLY: call_chem_at_all_substeps, chem_gasphase_on, cs_name 355 USE chemistry_model_mod, & 356 ONLY: chem_boundary_conds, chem_depo, chem_integrate, & 357 chem_prognostic_equations, chem_species, & 358 nspec, nvar, spc_names 355 359 356 360 USE control_parameters, & … … 495 499 496 500 IF ( intermediate_timestep_count == 1 .OR. & 497 call_chem_at_all_substeps ) THEN 498 CALL chem_integrate (i,j) 501 call_chem_at_all_substeps ) THEN 502 CALL chem_integrate (i,j) 503 IF ( do_depo ) THEN 504 CALL chem_depo(i,j) 505 ENDIF 499 506 ENDIF 500 507 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
