Changeset 4012 for palm/trunk/SOURCE/salsa_mod.f90
- Timestamp:
- May 31, 2019 3:19:05 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/salsa_mod.f90 (modified) (80 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/salsa_mod.f90
r3956 r4012 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Merge salsa branch to trunk. List of changes: 29 ! - Error corrected in distr_update that resulted in the aerosol number size 30 ! distribution not converging if the concentration was nclim. 31 ! - Added a separate output for aerosol liquid water (s_H2O) 32 ! - aerosol processes for a size bin are now calculated only if the aerosol 33 ! number of concentration of that bin is > 2*nclim 34 ! - An initialisation error in the subroutine "deposition" corrected and the 35 ! subroutine reformatted. 36 ! - stuff from salsa_util_mod.f90 moved into salsa_mod.f90 37 ! - calls for closing the netcdf input files added 38 ! 39 ! 3956 2019-05-07 12:32:52Z monakurppa 28 40 ! - Conceptual bug in depo_surf correct for urban and land surface model 29 41 ! - Subroutine salsa_tendency_ij optimized. … … 147 159 !> @todo emission mode "parameterized", i.e. based on street type 148 160 !> @todo Allow insoluble emissions 149 !> @todo two-way nesting is not working properly161 !> @todo Apply flux limiter in prognostic equations 150 162 !------------------------------------------------------------------------------! 151 163 MODULE salsa_mod 152 164 153 USE basic_constants_and_equations_mod, &165 USE basic_constants_and_equations_mod, & 154 166 ONLY: c_p, g, p_0, pi, r_d 155 167 156 USE chem_gasphase_mod, &168 USE chem_gasphase_mod, & 157 169 ONLY: nspec, nvar, spc_names 158 170 159 USE chem_modules, &171 USE chem_modules, & 160 172 ONLY: call_chem_at_all_substeps, chem_gasphase_on, chem_species 161 173 162 174 USE control_parameters 163 175 164 USE indices, & 165 ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, & 166 nzb_s_inner, nz, nzt, wall_flags_0 176 USE indices, & 177 ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nz, nzt, wall_flags_0 167 178 168 179 USE kinds 169 180 181 USE netcdf_data_input_mod, & 182 ONLY: chem_emis_att_type, chem_emis_val_type 183 170 184 USE pegrid 171 185 172 USE salsa_util_mod 173 174 USE statistics, & 186 USE statistics, & 175 187 ONLY: sums_salsa_ws_l 176 188 … … 452 464 (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) 453 465 REAL(wp), DIMENSION(nmod) :: surface_aerosol_flux = & 454 (/1.04e+11_wp, 3.23E+10_wp, 5.4E+6_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)466 (/1.0E+8_wp, 1.0E+9_wp, 1.0E+5_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) 455 467 456 468 REAL(wp), DIMENSION(:), ALLOCATABLE :: bin_low_limits !< to deliver information about … … 464 476 ! 465 477 !-- SALSA derived datatypes: 478 ! 479 !-- Component index 480 TYPE component_index 481 CHARACTER(len=3), ALLOCATABLE :: comp(:) !< Component name 482 INTEGER(iwp) :: ncomp !< Number of components 483 INTEGER(iwp), ALLOCATABLE :: ind(:) !< Component index 484 END TYPE component_index 466 485 ! 467 486 !-- For matching LSM and USM surface types and the deposition module surface types … … 607 626 TYPE(salsa_emission_value_type) :: aero_emission !< emission values 608 627 TYPE(salsa_emission_mode_type) :: def_modes !< default emission modes 628 629 TYPE(chem_emis_att_type) :: chem_emission_att !< chemistry emission attributes 630 TYPE(chem_emis_val_type) :: chem_emission !< chemistry emission values 609 631 610 632 TYPE(t_section), DIMENSION(:), ALLOCATABLE :: aero !< local aerosol properties … … 979 1001 SUBROUTINE salsa_header( io ) 980 1002 1003 USE indices, & 1004 ONLY: nx, ny, nz 1005 981 1006 IMPLICIT NONE 982 1007 … … 987 1012 WRITE( io, 2 ) skip_time_do_salsa 988 1013 WRITE( io, 3 ) dt_salsa 989 WRITE( io, 4 ) SHAPE( aerosol_number(1)%conc ), nbins_aerosol1014 WRITE( io, 4 ) nz, ny, nx, nbins_aerosol 990 1015 IF ( advect_particle_water ) THEN 991 WRITE( io, 5 ) SHAPE( aerosol_mass(1)%conc ), ncomponents_mass*nbins_aerosol, &1016 WRITE( io, 5 ) SHAPE( aerosol_mass(1)%conc ), ncomponents_mass*nbins_aerosol, & 992 1017 advect_particle_water 993 1018 ELSE … … 1154 1179 aerosol_number(i)%init(nzb:nzt+1), & 1155 1180 aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 1181 aerosol_number(i)%init = nclim 1156 1182 IF ( include_emission .OR. ( nldepo .AND. nldepo_surf ) ) THEN 1157 1183 ALLOCATE( aerosol_number(i)%source(nys:nyn,nxl:nxr) ) … … 1179 1205 aerosol_mass(i)%init(nzb:nzt+1), & 1180 1206 aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 1207 aerosol_mass(i)%init = mclim 1181 1208 IF ( include_emission .OR. ( nldepo .AND. nldepo_surf ) ) THEN 1182 1209 ALLOCATE( aerosol_mass(i)%source(nys:nyn,nxl:nxr) ) … … 1283 1310 salsa_gas(i)%init(nzb:nzt+1), & 1284 1311 salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 1312 salsa_gas(i)%init = nclim 1285 1313 IF ( include_emission ) ALLOCATE( salsa_gas(i)%source(nys:nys,nxl:nxr) ) 1286 1314 ENDDO … … 1402 1430 IF ( nldepo ) sedim_vd = 0.0_wp 1403 1431 1404 DO ib = 1, nbins_aerosol1405 IF ( .NOT. read_restart_data_salsa ) aerosol_number(ib)%conc = nclim1406 aerosol_number(ib)%conc_p = 0.0_wp1407 aerosol_number(ib)%tconc_m = 0.0_wp1408 aerosol_number(ib)%flux_s = 0.0_wp1409 aerosol_number(ib)%diss_s = 0.0_wp1410 aerosol_number(ib)%flux_l = 0.0_wp1411 aerosol_number(ib)%diss_l = 0.0_wp1412 aerosol_number(ib)%init = nclim1413 aerosol_number(ib)%sums_ws_l = 0.0_wp1414 ENDDO1415 DO ic = 1, ncomponents_mass*nbins_aerosol1416 IF ( .NOT. read_restart_data_salsa ) aerosol_mass(ic)%conc = mclim1417 aerosol_mass(ic)%conc_p = 0.0_wp1418 aerosol_mass(ic)%tconc_m = 0.0_wp1419 aerosol_mass(ic)%flux_s = 0.0_wp1420 aerosol_mass(ic)%diss_s = 0.0_wp1421 aerosol_mass(ic)%flux_l = 0.0_wp1422 aerosol_mass(ic)%diss_l = 0.0_wp1423 aerosol_mass(ic)%init = mclim1424 aerosol_mass(ic)%sums_ws_l = 0.0_wp1425 ENDDO1426 1427 1432 IF ( .NOT. salsa_gases_from_chem ) THEN 1433 IF ( .NOT. read_restart_data_salsa ) THEN 1434 salsa_gas(1)%conc = h2so4_init 1435 salsa_gas(2)%conc = hno3_init 1436 salsa_gas(3)%conc = nh3_init 1437 salsa_gas(4)%conc = ocnv_init 1438 salsa_gas(5)%conc = ocsv_init 1439 ENDIF 1428 1440 DO ig = 1, ngases_salsa 1429 1441 salsa_gas(ig)%conc_p = 0.0_wp … … 1434 1446 salsa_gas(ig)%diss_l = 0.0_wp 1435 1447 salsa_gas(ig)%sums_ws_l = 0.0_wp 1448 salsa_gas(ig)%conc_p = salsa_gas(ig)%conc 1436 1449 ENDDO 1437 IF ( .NOT. read_restart_data_salsa ) THEN 1438 salsa_gas(1)%conc = h2so4_init 1439 salsa_gas(2)%conc = hno3_init 1440 salsa_gas(3)%conc = nh3_init 1441 salsa_gas(4)%conc = ocnv_init 1442 salsa_gas(5)%conc = ocsv_init 1443 ENDIF 1444 ! 1445 !-- Set initial value for gas compound tracers and initial values 1450 ! 1451 !-- Set initial value for gas compound tracer 1446 1452 salsa_gas(1)%init = h2so4_init 1447 1453 salsa_gas(2)%init = hno3_init … … 1466 1472 !-- Initialise location-dependent aerosol size distributions and chemical compositions: 1467 1473 CALL aerosol_init 1468 ! 1474 1469 1475 !-- Initalisation run of SALSA + calculate the vertical top index of the topography 1470 1476 DO i = nxl, nxr … … 1477 1483 ENDDO 1478 1484 ENDDO 1485 1486 DO ib = 1, nbins_aerosol 1487 aerosol_number(ib)%conc_p = aerosol_number(ib)%conc 1488 aerosol_number(ib)%tconc_m = 0.0_wp 1489 aerosol_number(ib)%flux_s = 0.0_wp 1490 aerosol_number(ib)%diss_s = 0.0_wp 1491 aerosol_number(ib)%flux_l = 0.0_wp 1492 aerosol_number(ib)%diss_l = 0.0_wp 1493 aerosol_number(ib)%sums_ws_l = 0.0_wp 1494 ENDDO 1495 DO ic = 1, ncomponents_mass*nbins_aerosol 1496 aerosol_mass(ic)%conc_p = aerosol_mass(ic)%conc 1497 aerosol_mass(ic)%tconc_m = 0.0_wp 1498 aerosol_mass(ic)%flux_s = 0.0_wp 1499 aerosol_mass(ic)%diss_s = 0.0_wp 1500 aerosol_mass(ic)%flux_l = 0.0_wp 1501 aerosol_mass(ic)%diss_l = 0.0_wp 1502 aerosol_mass(ic)%sums_ws_l = 0.0_wp 1503 ENDDO 1504 ! 1479 1505 ! 1480 1506 !-- Initialise the deposition scheme and surface types … … 1587 1613 1588 1614 USE netcdf_data_input_mod, & 1589 ONLY: get_attribute, get_variable, netcdf_data_input_get_dimension_length, open_read_file 1615 ONLY: close_input_file, get_attribute, get_variable, & 1616 netcdf_data_input_get_dimension_length, open_read_file 1590 1617 1591 1618 IMPLICIT NONE … … 1615 1642 1616 1643 REAL(wp), DIMENSION(nbins_aerosol) :: core !< size of the bin mid aerosol particle 1617 REAL(wp), DIMENSION(nbins_aerosol) :: nsect !< size distribution (#/m3)1618 1644 1619 1645 REAL(wp), DIMENSION(0:nz+1) :: pnf2a !< number fraction in 2a … … 1638 1664 ! 1639 1665 !-- Set concentrations to zero 1640 nsect(:) = 0.0_wp1641 1666 pndist(:,:) = 0.0_wp 1642 1667 pnf2a(:) = nf2a … … 1667 1692 ! 1668 1693 !-- Allocate memory 1669 ALLOCATE( pr_z(1:pr_nz), pr_mass_fracs_a(nzb:nzt+1,pr_ncc), &1694 ALLOCATE( pr_z(1:pr_nz), pr_mass_fracs_a(nzb:nzt+1,pr_ncc), & 1670 1695 pr_mass_fracs_b(nzb:nzt+1,pr_ncc) ) 1671 1696 pr_mass_fracs_a = 0.0_wp … … 1806 1831 ' for SALSA missing!' 1807 1832 CALL message( 'salsa_mod: aerosol_init', 'PA0607', 1, 2, 0, 6, 0 ) 1808 1833 ! 1834 !-- Close input file 1835 CALL close_input_file( id_dyn ) 1809 1836 ENDIF ! netcdf_extend 1810 1837 … … 1892 1919 salsa_gas(ig)%init(nzb) = salsa_gas(ig)%init(nzb+1) 1893 1920 salsa_gas(ig)%init(nzt+1) = salsa_gas(ig)%init(nzt) 1894 DO k = nzb, nzt+1 1895 salsa_gas(ig)%conc(k,:,:) = salsa_gas(ig)%init(k) 1896 ENDDO 1921 IF ( .NOT. read_restart_data_salsa ) THEN 1922 DO k = nzb, nzt+1 1923 salsa_gas(ig)%conc(k,:,:) = salsa_gas(ig)%init(k) 1924 ENDDO 1925 ENDIF 1897 1926 ENDDO 1898 1927 … … 1901 1930 ' for SALSA missing!' 1902 1931 CALL message( 'salsa_mod: aerosol_init', 'PA0610', 1, 2, 0, 6, 0 ) 1932 ! 1933 !-- Close input file 1934 CALL close_input_file( id_dyn ) 1903 1935 ENDIF ! netcdf_extend 1904 1936 #else … … 1940 1972 !-- Region 1: 1941 1973 DO ib = start_subrange_1a, end_subrange_1a 1942 aerosol_number(ib)%conc(k,j,i) = pndist(k,ib) * flag 1974 IF ( .NOT. read_restart_data_salsa ) THEN 1975 aerosol_number(ib)%conc(k,j,i) = pndist(k,ib) * flag 1976 ENDIF 1943 1977 IF ( prunmode == 1 ) THEN 1944 1978 aerosol_number(ib)%init = pndist(:,ib) … … 1949 1983 IF ( nreg > 1 ) THEN 1950 1984 DO ib = start_subrange_2a, end_subrange_2a 1951 aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) * pndist(k,ib) * flag 1985 IF ( .NOT. read_restart_data_salsa ) THEN 1986 aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) * pndist(k,ib) * flag 1987 ENDIF 1952 1988 IF ( prunmode == 1 ) THEN 1953 1989 aerosol_number(ib)%init = MAX( 0.0_wp, nf2a ) * pndist(:,ib) … … 1957 1993 DO ib = start_subrange_2b, end_subrange_2b 1958 1994 IF ( pnf2a(k) < 1.0_wp ) THEN 1959 aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pnf2a(k) ) * & 1960 pndist(k,ib) * flag 1995 IF ( .NOT. read_restart_data_salsa ) THEN 1996 aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pnf2a(k) ) * & 1997 pndist(k,ib) * flag 1998 ENDIF 1961 1999 IF ( prunmode == 1 ) THEN 1962 2000 aerosol_number(ib)%init = MAX( 0.0_wp, 1.0_wp - nf2a ) * pndist(:,ib) … … 1976 2014 ib = start_subrange_1a 1977 2015 DO ic = ss, ee 1978 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) * pndist(k,ib)& 1979 * core(ib) * arhoh2so4 * flag 2016 IF ( .NOT. read_restart_data_salsa ) THEN 2017 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) * & 2018 pndist(k,ib) * core(ib) * arhoh2so4 * flag 2019 ENDIF 1980 2020 IF ( prunmode == 1 ) THEN 1981 2021 aerosol_mass(ic)%init(k) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) * pndist(k,ib) & … … 1991 2031 ee = ( index_oc - 1 ) * nbins_aerosol + end_subrange_1a !< end 1992 2032 ib = start_subrange_1a 1993 DO ic = ss, ee 1994 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) * & 1995 core(ib) * arhooc * flag 2033 DO ic = ss, ee 2034 IF ( .NOT. read_restart_data_salsa ) THEN 2035 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) * & 2036 core(ib) * arhooc * flag 2037 ENDIF 1996 2038 IF ( prunmode == 1 ) THEN 1997 2039 aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) * & … … 2143 2185 ib = start_subrange_2a 2144 2186 DO ic = ss, ee 2145 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib) * & 2146 pcore(ib) * prho * flag 2187 IF ( .NOT. read_restart_data_salsa ) THEN 2188 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib)& 2189 * pcore(ib) * prho * flag 2190 ENDIF 2147 2191 IF ( prunmode == 1 ) THEN 2148 2192 aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib) * & … … 2158 2202 ib = start_subrange_2a 2159 2203 DO ic = ss, ee 2160 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k) ) *& 2161 pndist(k,ib) * pcore(ib) * prho * flag 2204 IF ( .NOT. read_restart_data_salsa ) THEN 2205 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k))& 2206 * pndist(k,ib) * pcore(ib) * prho * flag 2207 ENDIF 2162 2208 IF ( prunmode == 1 ) THEN 2163 2209 aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k) ) * & … … 2810 2856 ENDDO 2811 2857 ! 2812 !-- On EACH call of salsa_driver, calculate the ambient sizes of 2813 !-- particles by equilibrating soluble fraction of particles with water 2814 !-- using the ZSR method. 2858 !-- Calculate the ambient sizes of particles by equilibrating soluble fraction of particles with 2859 !-- water using the ZSR method. 2815 2860 in_rh = in_cw(k) / in_cs(k) 2816 2861 IF ( prunmode==1 .OR. .NOT. advect_particle_water ) THEN … … 2865 2910 !-- Calculate changes in concentrations 2866 2911 DO ib = 1, nbins_aerosol 2867 aerosol_number(ib)%conc(k,j,i) = aerosol_number(ib)%conc(k,j,i) + ( lo_aero(ib)%numc - 2912 aerosol_number(ib)%conc(k,j,i) = aerosol_number(ib)%conc(k,j,i) + ( lo_aero(ib)%numc - & 2868 2913 aero_old(ib)%numc ) * flag 2869 2914 ENDDO … … 2875 2920 ic = 1 2876 2921 DO ss = str, endi 2877 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) - 2922 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -& 2878 2923 aero_old(ic)%volc(vc) ) * arhoh2so4 * flag 2879 2924 ic = ic+1 … … 2887 2932 ic = 1 2888 2933 DO ss = str, endi 2889 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) - 2934 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -& 2890 2935 aero_old(ic)%volc(vc) ) * arhooc * flag 2891 2936 ic = ic+1 … … 2899 2944 ic = 1 + end_subrange_1a 2900 2945 DO ss = str, endi 2901 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) - 2946 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -& 2902 2947 aero_old(ic)%volc(vc) ) * arhobc * flag 2903 2948 ic = ic+1 … … 2911 2956 ic = 1 + end_subrange_1a 2912 2957 DO ss = str, endi 2913 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) - 2958 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -& 2914 2959 aero_old(ic)%volc(vc) ) * arhodu * flag 2915 2960 ic = ic+1 … … 2923 2968 ic = 1 + end_subrange_1a 2924 2969 DO ss = str, endi 2925 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) - 2970 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -& 2926 2971 aero_old(ic)%volc(vc) ) * arhoss * flag 2927 2972 ic = ic+1 … … 2935 2980 ic = 1 2936 2981 DO ss = str, endi 2937 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) - 2982 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -& 2938 2983 aero_old(ic)%volc(vc) ) * arhohno3 * flag 2939 2984 ic = ic+1 … … 2947 2992 ic = 1 2948 2993 DO ss = str, endi 2949 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) - 2994 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -& 2950 2995 aero_old(ic)%volc(vc) ) * arhonh3 * flag 2951 2996 ic = ic+1 … … 2960 3005 ic = 1 2961 3006 DO ss = str, endi 2962 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) - 3007 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -& 2963 3008 aero_old(ic)%volc(vc) ) * arhoh2o * flag 2964 IF ( prunmode == 1 ) THEN 2965 aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), & 2966 aerosol_mass(ss)%conc(k,j,i) ) 2967 IF ( k == nzb+1 ) THEN 2968 aerosol_mass(ss)%init(k-1) = 0.0_wp 2969 ELSEIF ( k == nzt ) THEN 2970 aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k) 2971 ENDIF 3009 ic = ic+1 3010 ENDDO 3011 ENDIF 3012 IF ( prunmode == 1 ) THEN 3013 nc_h2o = get_index( prtcl,'H2O' ) 3014 vc = 8 3015 str = ( nc_h2o-1 ) * nbins_aerosol + 1 3016 endi = nc_h2o * nbins_aerosol 3017 ic = 1 3018 DO ss = str, endi 3019 aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), ( lo_aero(ic)%volc(vc) - & 3020 aero_old(ic)%volc(vc) ) * arhoh2o ) 3021 IF ( k == nzb+1 ) THEN 3022 aerosol_mass(ss)%init(k-1) = aerosol_mass(ss)%init(k) 3023 ELSEIF ( k == nzt ) THEN 3024 aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k) 3025 aerosol_mass(ss)%conc(k+1,j,i) = aerosol_mass(ss)%init(k) 2972 3026 ENDIF 2973 3027 ic = ic+1 … … 3446 3500 REAL(wp) :: avis !< molecular viscocity of air (kg/(m*s)) 3447 3501 REAL(wp) :: beta_im !< parameter for turbulent impaction 3448 REAL(wp) :: beta !< Cunningham slip-flow correction factor3449 3502 REAL(wp) :: c_brownian_diff !< coefficient for Brownian diffusion 3450 3503 REAL(wp) :: c_impaction !< coefficient for inertial impaction … … 3453 3506 REAL(wp) :: depo !< deposition velocity (m/s) 3454 3507 REAL(wp) :: gamma !< parameter, Table 3 in Z01 3455 REAL(wp) :: Kn !< Knudsen number3456 3508 REAL(wp) :: lambda !< molecular mean free path (m) 3457 3509 REAL(wp) :: mdiff !< particle diffusivity coefficient … … 3461 3513 REAL(wp) :: pdn !< particle density (kg/m3) 3462 3514 REAL(wp) :: ustar !< friction velocity (m/s) 3463 REAL(wp) :: va !< thermal speed of an air molecule (m/s) 3464 REAL(wp) :: zdwet !< wet diameter (m) 3515 REAL(wp) :: va !< thermal speed of an air molecule (m/s) 3465 3516 3466 3517 REAL(wp), INTENT(in) :: adn !< air density (kg/m3) … … 3471 3522 REAL(wp), INTENT(inout) :: kvis !< kinematic viscosity of air (m2/s) 3472 3523 3524 REAL(wp), DIMENSION(nbins_aerosol) :: beta !< Cunningham slip-flow correction factor 3525 REAL(wp), DIMENSION(nbins_aerosol) :: Kn !< Knudsen number 3526 REAL(wp), DIMENSION(nbins_aerosol) :: zdwet !< wet diameter (m) 3527 3473 3528 REAL(wp), DIMENSION(:), INTENT(inout) :: schmidt_num !< particle Schmidt number 3474 3529 REAL(wp), DIMENSION(:), INTENT(inout) :: vc !< critical fall speed i.e. settling velocity of … … 3494 3549 lambda = 2.0_wp * avis / ( adn * va ) 3495 3550 ! 3496 !-- Parameters for the land use category 'deciduous broadleaf trees'(Table 3) 3497 alpha = alpha_z01(depo_pcm_type_num) 3498 gamma = gamma_z01(depo_pcm_type_num) 3499 par_a = A_z01(depo_pcm_type_num, season) * 1.0E-3_wp 3500 ! 3501 !-- Deposition efficiencies from Table 1. Constants from Table 2. 3502 par_l = l_p10(depo_pcm_type_num) * 0.01_wp 3503 c_brownian_diff = c_b_p10(depo_pcm_type_num) 3504 c_interception = c_in_p10(depo_pcm_type_num) 3505 c_impaction = c_im_p10(depo_pcm_type_num) 3506 beta_im = beta_im_p10(depo_pcm_type_num) 3507 c_turb_impaction = c_it_p10(depo_pcm_type_num) 3508 3509 DO ib = 1, nbins_aerosol 3510 3511 IF ( paero(ib)%numc < nclim ) CYCLE 3512 zdwet = paero(ib)%dwet 3513 ! 3514 !-- Knudsen number (Eq. 15.23) 3515 Kn = MAX( 1.0E-2_wp, lambda / ( zdwet * 0.5_wp ) ) ! To avoid underflow 3516 ! 3517 !-- Cunningham slip-flow correction (Eq. 15.30) 3518 beta = 1.0_wp + Kn * ( 1.249_wp + 0.42_wp * EXP( -0.87_wp / Kn ) ) 3519 3520 !-- Particle diffusivity coefficient (Eq. 15.29) 3521 mdiff = ( abo * tk * beta ) / ( 3.0_wp * pi * avis * zdwet ) 3522 ! 3523 !-- Particle Schmidt number (Eq. 15.36) 3524 schmidt_num(ib) = kvis / mdiff 3525 ! 3526 !-- Critical fall speed i.e. settling velocity (Eq. 20.4) 3527 vc(ib) = MIN( 1.0_wp, terminal_vel( 0.5_wp * zdwet, pdn, adn, avis, beta) ) 3528 ! 3529 !-- Friction velocity for deposition on vegetation. Calculated following Prandtl (1925): 3530 IF ( lsdepo_pcm .AND. plant_canopy .AND. lad > 0.0_wp ) THEN 3551 !-- Particle wet diameter (m) 3552 zdwet = paero(:)%dwet 3553 ! 3554 !-- Knudsen number (Eq. 15.23) 3555 Kn = MAX( 1.0E-2_wp, lambda / ( zdwet * 0.5_wp ) ) ! To avoid underflow 3556 ! 3557 !-- Cunningham slip-flow correction (Eq. 15.30) 3558 beta = 1.0_wp + Kn * ( 1.249_wp + 0.42_wp * EXP( -0.87_wp / Kn ) ) 3559 ! 3560 !-- Critical fall speed i.e. settling velocity (Eq. 20.4) 3561 vc = MIN( 1.0_wp, zdwet**2 * ( pdn - adn ) * g * beta / ( 18.0_wp * avis ) ) 3562 ! 3563 !-- Deposition on vegetation 3564 IF ( lsdepo_pcm .AND. plant_canopy .AND. lad > 0.0_wp ) THEN 3565 ! 3566 !-- Parameters for the land use category 'deciduous broadleaf trees'(Table 3) 3567 alpha = alpha_z01(depo_pcm_type_num) 3568 gamma = gamma_z01(depo_pcm_type_num) 3569 par_a = A_z01(depo_pcm_type_num, season) * 1.0E-3_wp 3570 ! 3571 !-- Deposition efficiencies from Table 1. Constants from Table 2. 3572 par_l = l_p10(depo_pcm_type_num) * 0.01_wp 3573 c_brownian_diff = c_b_p10(depo_pcm_type_num) 3574 c_interception = c_in_p10(depo_pcm_type_num) 3575 c_impaction = c_im_p10(depo_pcm_type_num) 3576 beta_im = beta_im_p10(depo_pcm_type_num) 3577 c_turb_impaction = c_it_p10(depo_pcm_type_num) 3578 3579 DO ib = 1, nbins_aerosol 3580 3581 IF ( paero(ib)%numc < ( 2.0_wp * nclim ) ) CYCLE 3582 3583 !-- Particle diffusivity coefficient (Eq. 15.29) 3584 mdiff = ( abo * tk * beta(ib) ) / ( 3.0_wp * pi * avis * zdwet(ib) ) 3585 ! 3586 !-- Particle Schmidt number (Eq. 15.36) 3587 schmidt_num(ib) = kvis / mdiff 3588 ! 3589 !-- Friction velocity for deposition on vegetation. Calculated following Prandtl (1925): 3531 3590 ustar = SQRT( cdc ) * mag_u 3532 3591 SELECT CASE ( depo_pcm_par_num ) … … 3546 3605 paero(ib)%volc(ic) = paero(ib)%volc(ic) - depo * lad * paero(ib)%volc(ic) * dt_salsa 3547 3606 ENDDO 3548 ENDIF 3549 ENDDO 3607 ENDDO 3608 3609 ENDIF 3550 3610 3551 3611 END SUBROUTINE deposition … … 3779 3839 3780 3840 DO ib = 1, nbins_aerosol 3781 IF ( aerosol_number(ib)%conc(k,j,i) < = nclim .OR. schmidt_num(k+1,ib) < 1.0_wp )&3782 CYCLE3841 IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim ) .OR. & 3842 schmidt_num(k+1,ib) < 1.0_wp ) CYCLE 3783 3843 3784 3844 SELECT CASE ( depo_surf_par_num ) … … 3810 3870 3811 3871 DO ib = 1, nbins_aerosol 3812 IF ( aerosol_number(ib)%conc(k,j,i) < = nclim .OR. schmidt_num(k+1,ib) < 1.0_wp )&3813 CYCLE3872 IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim ) .OR. & 3873 schmidt_num(k+1,ib) < 1.0_wp ) CYCLE 3814 3874 3815 3875 SELECT CASE ( depo_surf_par_num ) … … 3841 3901 3842 3902 DO ib = 1, nbins_aerosol 3843 IF ( aerosol_number(ib)%conc(k,j,i) < = nclim .OR. schmidt_num(k+1,ib) < 1.0_wp )&3844 CYCLE3903 IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim ) .OR. & 3904 schmidt_num(k+1,ib) < 1.0_wp ) CYCLE 3845 3905 3846 3906 SELECT CASE ( depo_surf_par_num ) … … 3860 3920 3861 3921 DO ib = 1, nbins_aerosol 3922 IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim ) ) CYCLE 3862 3923 ! 3863 3924 !-- Calculate changes in surface fluxes due to dry deposition … … 3891 3952 3892 3953 DO ib = 1, nbins_aerosol 3893 IF ( aerosol_number(ib)%conc(k,j,i) < = nclim .OR. schmidt_num(k+1,ib) < 1.0_wp )&3894 CYCLE3954 IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim ) .OR. & 3955 schmidt_num(k+1,ib) < 1.0_wp ) CYCLE 3895 3956 3896 3957 SELECT CASE ( depo_surf_par_num ) … … 3932 3993 ! Description: 3933 3994 ! ------------ 3934 ! Function for calculating terminal velocities for different particles sizes.3935 !------------------------------------------------------------------------------!3936 REAL(wp) FUNCTION terminal_vel( radius, rhop, rhoa, visc, beta )3937 3938 IMPLICIT NONE3939 3940 REAL(wp), INTENT(in) :: beta !< Cunningham correction factor3941 REAL(wp), INTENT(in) :: radius !< particle radius (m)3942 REAL(wp), INTENT(in) :: rhop !< particle density (kg/m3)3943 REAL(wp), INTENT(in) :: rhoa !< air density (kg/m3)3944 REAL(wp), INTENT(in) :: visc !< molecular viscosity of air (kg/(m*s))3945 3946 !3947 !-- Stokes law with Cunningham slip correction factor3948 terminal_vel = 4.0_wp * radius**2 * ( rhop - rhoa ) * g * beta / ( 18.0_wp * visc ) ! (m/s)3949 3950 END FUNCTION terminal_vel3951 3952 !------------------------------------------------------------------------------!3953 ! Description:3954 ! ------------3955 3995 !> Calculates particle loss and change in size distribution due to (Brownian) 3956 3996 !> coagulation. Only for particles with dwet < 30 micrometres. … … 4034 4074 !-- Aero-aero coagulation 4035 4075 DO mm = 1, end_subrange_2b ! smaller colliding particle 4036 IF ( paero(mm)%numc < nclim) CYCLE4076 IF ( paero(mm)%numc < ( 2.0_wp * nclim ) ) CYCLE 4037 4077 DO nn = mm, end_subrange_2b ! larger colliding particle 4038 IF ( paero(nn)%numc < nclim) CYCLE4078 IF ( paero(nn)%numc < ( 2.0_wp * nclim ) ) CYCLE 4039 4079 4040 4080 zdpart_mm = MIN( paero(mm)%dwet, 30.0E-6_wp ) ! Limit to 30 um … … 4053 4093 !-- Aerosols in subrange 1a: 4054 4094 DO ib = start_subrange_1a, end_subrange_1a 4055 IF ( paero(ib)%numc < nclim) CYCLE4095 IF ( paero(ib)%numc < ( 2.0_wp * nclim ) ) CYCLE 4056 4096 zminusterm = 0.0_wp 4057 4097 zplusterm(:) = 0.0_wp … … 4080 4120 !-- Aerosols in subrange 2a: 4081 4121 DO ib = start_subrange_2a, end_subrange_2a 4082 IF ( paero(ib)%numc < nclim) CYCLE4122 IF ( paero(ib)%numc < ( 2.0_wp * nclim ) ) CYCLE 4083 4123 zminusterm = 0.0_wp 4084 4124 zplusterm(:) = 0.0_wp … … 4128 4168 IF ( .NOT. no_insoluble ) THEN 4129 4169 DO ib = start_subrange_2b, end_subrange_2b 4130 IF ( paero(ib)%numc < nclim) CYCLE4170 IF ( paero(ib)%numc < ( 2.0_wp * nclim ) ) CYCLE 4131 4171 zminusterm = 0.0_wp 4132 4172 zplusterm(:) = 0.0_wp … … 4337 4377 REAL(wp), INTENT(inout) :: pcw !< Water vapor concentration (kg/m3) 4338 4378 4339 REAL(wp), DIMENSION(nbins_aerosol) 4340 REAL(wp), DIMENSION(nbins_aerosol) 4341 REAL(wp), DIMENSION(nbins_aerosol) :: zcolrate_ocnv !< collision rate of non-vol. OC(1/s)4342 REAL(wp), DIMENSION(start_subrange_1a+1) :: zdfpart !< particle diffusion coef ficient(m2/s)4343 REAL(wp), DIMENSION(nbins_aerosol) 4344 REAL(wp), DIMENSION(nbins_aerosol) 4379 REAL(wp), DIMENSION(nbins_aerosol) :: zbeta !< transitional correction factor 4380 REAL(wp), DIMENSION(nbins_aerosol) :: zcolrate !< collision rate (1/s) 4381 REAL(wp), DIMENSION(nbins_aerosol) :: zcolrate_ocnv !< collision rate of OCNV (1/s) 4382 REAL(wp), DIMENSION(start_subrange_1a+1) :: zdfpart !< particle diffusion coef. (m2/s) 4383 REAL(wp), DIMENSION(nbins_aerosol) :: zdvoloc !< change of organics volume 4384 REAL(wp), DIMENSION(nbins_aerosol) :: zdvolsa !< change of sulphate volume 4345 4385 REAL(wp), DIMENSION(2) :: zj3n3 !< Formation massrate of molecules 4346 4386 !< in nucleation, (molec/m3s), 4347 4387 !< 1: H2SO4 and 2: organic vapor 4348 REAL(wp), DIMENSION(nbins_aerosol) :: zknud !< particle Knudsen number4388 REAL(wp), DIMENSION(nbins_aerosol) :: zknud !< particle Knudsen number 4349 4389 4350 4390 TYPE(component_index), INTENT(in) :: prtcl !< Keeps track which substances are used … … 7095 7135 REAL(wp) :: znfrac !< number fraction to be moved to the larger bin 7096 7136 REAL(wp) :: zvfrac !< volume fraction to be moved to the larger bin 7097 REAL(wp) :: z Vexc !< Volume in the grown bin which exceeds the bin upper limit7098 REAL(wp) :: z Vihi !< particle volume at the high end of the bin7099 REAL(wp) :: z Vilo !< particle volume at the low end of the bin7137 REAL(wp) :: zvexc !< Volume in the grown bin which exceeds the bin upper limit 7138 REAL(wp) :: zvihi !< particle volume at the high end of the bin 7139 REAL(wp) :: zvilo !< particle volume at the low end of the bin 7100 7140 REAL(wp) :: zvpart !< particle volume (m3) 7101 REAL(wp) :: zVrat !< volume ratio of a size bin 7141 REAL(wp) :: zvrat !< volume ratio of a size bin 7142 7143 real(wp), dimension(nbins_aerosol) :: dummy 7102 7144 7103 7145 TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) :: paero !< aerosol properties … … 7106 7148 zvfrac = 0.0_wp 7107 7149 within_bins = .FALSE. 7150 7151 dummy = paero(:)%numc 7108 7152 ! 7109 7153 !-- Check if the volume of the bin is within bin limits after update … … 7150 7194 ! 7151 7195 !-- Volume ratio of the size bin 7152 z Vrat = paero(ib)%vhilim / paero(ib)%vlolim7196 zvrat = paero(ib)%vhilim / paero(ib)%vlolim 7153 7197 ! 7154 7198 !-- Particle volume at the low end of the bin 7155 z Vilo = 2.0_wp * zvpart / ( 1.0_wp + zVrat )7199 zvilo = 2.0_wp * zvpart / ( 1.0_wp + zvrat ) 7156 7200 ! 7157 7201 !-- Particle volume at the high end of the bin 7158 z Vihi = zVrat * zVilo7202 zvihi = zvrat * zvilo 7159 7203 ! 7160 7204 !-- Volume in the grown bin which exceeds the bin upper limit 7161 z Vexc = 0.5_wp * ( zVihi + paero(ib)%vhilim )7205 zvexc = 0.5_wp * ( zvihi + paero(ib)%vhilim ) 7162 7206 ! 7163 7207 !-- Number fraction to be moved to the larger bin 7164 znfrac = MIN( 1.0_wp, ( z Vihi - paero(ib)%vhilim) / ( zVihi - zVilo ) )7208 znfrac = MIN( 1.0_wp, ( zvihi - paero(ib)%vhilim) / ( zvihi - zvilo ) ) 7165 7209 ! 7166 7210 !-- Volume fraction to be moved to the larger bin 7167 zvfrac = MIN( 0.99_wp, znfrac * z Vexc / zvpart )7211 zvfrac = MIN( 0.99_wp, znfrac * zvexc / zvpart ) 7168 7212 IF ( zvfrac < 0.0_wp ) THEN 7169 7213 message_string = 'Error: zvfrac < 0' … … 7175 7219 ! 7176 7220 !-- Volume (cm3/cm3) 7177 paero(mm)%volc(:) = paero(mm)%volc(:) + znfrac * paero(ib)%numc * z Vexc * &7178 paero(ib)%volc(:) / SUM( paero(ib)%volc( :) )7179 paero(ib)%volc(:) = paero(ib)%volc(:) - znfrac * paero(ib)%numc * z Vexc * &7180 paero(ib)%volc(:) / SUM( paero(ib)%volc( :) )7221 paero(mm)%volc(:) = paero(mm)%volc(:) + znfrac * paero(ib)%numc * zvexc * & 7222 paero(ib)%volc(:) / SUM( paero(ib)%volc(1:7) ) 7223 paero(ib)%volc(:) = paero(ib)%volc(:) - znfrac * paero(ib)%numc * zvexc * & 7224 paero(ib)%volc(:) / SUM( paero(ib)%volc(1:7) ) 7181 7225 7182 7226 !-- Number concentration (#/m3) … … 7266 7310 !-- Calculate total mass concentration per bin 7267 7311 mcsum = 0.0_wp 7268 DO ic = 1, nc c7312 DO ic = 1, ncomponents_mass 7269 7313 icc = ( ic - 1 ) * nbins_aerosol + ib 7270 7314 mcsum = mcsum + aerosol_mass(icc)%conc(:,j,i) * flag 7315 aerosol_mass(icc)%conc(:,j,i) = MAX( mclim, aerosol_mass(icc)%conc(:,j,i) ) * flag 7271 7316 ENDDO 7272 7317 ! 7273 7318 !-- Check that number and mass concentration match qualitatively 7274 IF ( ANY ( aerosol_number(ib)%conc(:,j,i) >=nclim .AND. mcsum <= 0.0_wp ) ) THEN7319 IF ( ANY( aerosol_number(ib)%conc(:,j,i) > nclim .AND. mcsum <= 0.0_wp ) ) THEN 7275 7320 DO k = nzb+1, nzt 7276 7321 IF ( aerosol_number(ib)%conc(k,j,i) >= nclim .AND. mcsum(k) <= 0.0_wp ) THEN … … 7349 7394 flag_zddry > 0.0_wp ) 7350 7395 ENDDO 7351 aerosol_number(ib)%conc(:,j,i) = MERGE( nclim , aerosol_number(ib)%conc(:,j,i),&7396 aerosol_number(ib)%conc(:,j,i) = MERGE( nclim * flag, aerosol_number(ib)%conc(:,j,i), & 7352 7397 flag_zddry > 0.0_wp ) 7353 7398 ra_dry(:,j,i,ib) = MAX( 1.0E-10_wp, 0.5_wp * zddry ) … … 8387 8432 8388 8433 USE netcdf_data_input_mod, & 8389 ONLY: check_existence, get_attribute, get_variable, inquire_num_variables, & 8390 inquire_variable_names, netcdf_data_input_get_dimension_length, open_read_file 8434 ONLY: check_existence, close_input_file, get_attribute, get_variable, & 8435 inquire_num_variables, inquire_variable_names, & 8436 netcdf_data_input_get_dimension_length, open_read_file 8391 8437 8392 8438 USE surface_mod, & … … 8482 8528 8483 8529 #if defined( __netcdf ) 8530 ! 8531 !-- Check existence of PIDS_SALSA file 8532 INQUIRE( FILE = TRIM( input_file_salsa ) // TRIM( coupling_char ), EXIST = netcdf_extend ) 8533 IF ( .NOT. netcdf_extend ) THEN 8534 message_string = 'Input file '// TRIM( input_file_salsa ) // TRIM( coupling_char )& 8535 // ' missing!' 8536 CALL message( 'salsa_emission_setup', 'PA0629', 1, 2, 0, 6, 0 ) 8537 ENDIF 8538 ! 8539 !-- Open file in read-only mode 8540 CALL open_read_file( TRIM( input_file_salsa ) // TRIM( coupling_char ), id_salsa ) 8541 8484 8542 IF ( init ) THEN 8485 !8486 !-- Check existence of PIDS_SALSA file8487 INQUIRE( FILE = TRIM( input_file_salsa ) // TRIM( coupling_char ), &8488 EXIST = netcdf_extend )8489 IF ( .NOT. netcdf_extend ) THEN8490 message_string = 'Input file '// TRIM( input_file_salsa ) // TRIM( coupling_char )&8491 // ' missing!'8492 CALL message( 'salsa_emission_setup', 'PA0629', 1, 2, 0, 6, 0 )8493 ENDIF8494 !8495 !-- Open file in read-only mode8496 CALL open_read_file( TRIM( input_file_salsa ) // TRIM( coupling_char ), id_salsa )8497 8543 ! 8498 8544 !-- Read the index and name of chemical components … … 8822 8868 8823 8869 8824 DEALLOCATE( source_array )8870 DEALLOCATE( nsect_emission, source_array ) 8825 8871 ! 8826 8872 !-- Pre-processed: … … 8857 8903 8858 8904 ENDIF 8859 8905 ! 8906 !-- Close input file 8907 CALL close_input_file( id_salsa ) 8860 8908 #else 8861 8909 message_string = 'salsa_emission_mode = "read_from_file", but preprocessor directive ' //& … … 9073 9121 9074 9122 USE netcdf_data_input_mod, & 9075 ONLY: check_existence, c hem_emis_att_type, chem_emis_val_type, get_attribute,&9076 get_variable, inquire_num_variables, inquire_variable_names,&9123 ONLY: check_existence, close_input_file, get_attribute, get_variable, & 9124 inquire_num_variables, inquire_variable_names, & 9077 9125 netcdf_data_input_get_dimension_length, open_read_file 9078 9126 … … 9098 9146 REAL(wp), DIMENSION(:), ALLOCATABLE :: time_factor !< emission time factor 9099 9147 9100 TYPE(chem_emis_att_type) :: chem_emission_att !< chemistry emission attributes9101 TYPE(chem_emis_val_type) :: chem_emission !< chemistry emission values9102 9103 9148 ! 9104 9149 !-- Reset surface fluxes … … 9108 9153 9109 9154 #if defined( __netcdf ) 9155 ! 9156 !-- Check existence of PIDS_CHEM file 9157 INQUIRE( FILE = 'PIDS_CHEM' // TRIM( coupling_char ), EXIST = netcdf_extend ) 9158 IF ( .NOT. netcdf_extend ) THEN 9159 message_string = 'Input file PIDS_CHEM' // TRIM( coupling_char ) // ' missing!' 9160 CALL message( 'salsa_gas_emission_setup', 'PA0640', 1, 2, 0, 6, 0 ) 9161 ENDIF 9162 ! 9163 !-- Open file in read-only mode 9164 CALL open_read_file( 'PIDS_CHEM' // TRIM( coupling_char ), id_chem ) 9165 9110 9166 IF ( init ) THEN 9111 !9112 !-- Check existence of PIDS_CHEM file9113 INQUIRE( FILE = 'PIDS_CHEM' // TRIM( coupling_char ), EXIST = netcdf_extend )9114 IF ( .NOT. netcdf_extend ) THEN9115 message_string = 'Input file PIDS_CHEM' // TRIM( coupling_char ) // ' missing!'9116 CALL message( 'salsa_gas_emission_setup', 'PA0640', 1, 2, 0, 6, 0 )9117 ENDIF9118 !9119 !-- Open file in read-only mode9120 CALL open_read_file( 'PIDS_CHEM' // TRIM( coupling_char ), id_chem )9121 9167 ! 9122 9168 !-- Read the index and name of chemical components … … 9301 9347 DEALLOCATE( chem_emission%default_emission_data ) 9302 9348 ENDIF 9349 ! 9350 !-- Close input file 9351 CALL close_input_file( id_chem ) 9352 9303 9353 #else 9304 9354 message_string = 'salsa_emission_mode = "read_from_file", but preprocessor directive ' // & … … 9347 9397 TYPE(surf_type), INTENT(inout) :: surface !< respective surface type 9348 9398 9399 conv = 1.0_wp 9349 9400 use_time_fac = PRESENT( time_fac ) 9350 9401 … … 9765 9816 ENDDO 9766 9817 9767 CASE ( 's_BC', 's_DU', 's_ H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )9818 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 9768 9819 IF ( is_used( prtcl, TRIM( variable(3:) ) ) ) THEN 9769 9820 found_index = get_index( prtcl, TRIM( variable(3:) ) ) … … 9786 9837 ENDDO 9787 9838 ENDIF 9839 9840 CASE ( 's_H2O' ) 9841 found_index = get_index( prtcl,'H2O' ) 9842 DO i = nxlg, nxrg 9843 DO j = nysg, nyng 9844 DO k = nzb, nzt+1 9845 DO ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol 9846 to_be_resorted(k,j,i) = to_be_resorted(k,j,i) + & 9847 aerosol_mass(ic)%conc(k,j,i) 9848 ENDDO 9849 ENDDO 9850 ENDDO 9851 ENDDO 9788 9852 9789 9853 CASE DEFAULT … … 9885 9949 ENDDO 9886 9950 9887 CASE ( 's_BC', 's_DU', 's_ H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )9951 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 9888 9952 IF ( is_used( prtcl, TRIM( variable(3:) ) ) ) THEN 9889 found_index = get_index( prtcl, TRIM( variable(3:) ) )9890 9953 IF ( TRIM( variable(3:) ) == 'BC' ) to_be_resorted => s_bc_av 9891 9954 IF ( TRIM( variable(3:) ) == 'DU' ) to_be_resorted => s_du_av … … 9904 9967 ENDDO 9905 9968 ENDIF 9969 9970 CASE ( 's_H2O' ) 9971 to_be_resorted => s_h2o_av 9972 DO i = nxlg, nxrg 9973 DO j = nysg, nyng 9974 DO k = nzb, nzt+1 9975 to_be_resorted(k,j,i) = to_be_resorted(k,j,i) / & 9976 REAL( average_count_3d, KIND=wp ) 9977 ENDDO 9978 ENDDO 9979 ENDDO 9906 9980 9907 9981 END SELECT … … 10212 10286 IF ( mode == 'xy' ) grid = 'zu' 10213 10287 10214 CASE ( 's_BC', 's_DU', 's_ H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )10288 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 10215 10289 vari = TRIM( variable( 3:LEN( TRIM( variable ) ) - 3 ) ) 10216 10290 IF ( is_used( prtcl, vari ) ) THEN … … 10252 10326 IF ( mode == 'xy' ) grid = 'zu' 10253 10327 10328 CASE ( 's_H2O' ) 10329 found_index = get_index( prtcl, 'H2O' ) 10330 IF ( av == 0 ) THEN 10331 DO i = nxl, nxr 10332 DO j = nys, nyn 10333 DO k = nzb_do, nzt_do 10334 temp_bin = 0.0_wp 10335 DO ic = ( found_index-1 ) * nbins_aerosol+1, found_index * nbins_aerosol 10336 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10337 ENDDO 10338 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10339 BTEST( wall_flags_0(k,j,i), 0 ) ) 10340 ENDDO 10341 ENDDO 10342 ENDDO 10343 ELSE 10344 to_be_resorted => s_h2o_av 10345 DO i = nxl, nxr 10346 DO j = nys, nyn 10347 DO k = nzb_do, nzt_do 10348 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, & 10349 KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 10350 ENDDO 10351 ENDDO 10352 ENDDO 10353 ENDIF 10354 10355 IF ( mode == 'xy' ) grid = 'zu' 10356 10254 10357 CASE DEFAULT 10255 10358 found = .FALSE. … … 10534 10637 ENDIF 10535 10638 10536 CASE ( 's_BC', 's_DU', 's_ H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )10639 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 10537 10640 IF ( is_used( prtcl, TRIM( variable(3:) ) ) ) THEN 10538 10641 found_index = get_index( prtcl, TRIM( variable(3:) ) ) … … 10567 10670 ENDDO 10568 10671 ENDIF 10672 ENDIF 10673 10674 CASE ( 's_H2O' ) 10675 found_index = get_index( prtcl, 'H2O' ) 10676 IF ( av == 0 ) THEN 10677 DO i = nxl, nxr 10678 DO j = nys, nyn 10679 DO k = nzb_do, nzt_do 10680 temp_bin = 0.0_wp 10681 DO ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol 10682 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10683 ENDDO 10684 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10685 BTEST( wall_flags_0(k,j,i), 0 ) ) 10686 ENDDO 10687 ENDDO 10688 ENDDO 10689 ELSE 10690 to_be_resorted => s_h2o_av 10691 DO i = nxl, nxr 10692 DO j = nys, nyn 10693 DO k = nzb_do, nzt_do 10694 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, & 10695 KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 10696 ENDDO 10697 ENDDO 10698 ENDDO 10569 10699 ENDIF 10570 10700 … … 10962 11092 ENDIF 10963 11093 11094 CASE ( 's_H2O' ) 11095 IF ( av == 0 ) THEN 11096 found_index = get_index( prtcl, 'H2O' ) 11097 DO i = nxl, nxr 11098 DO j = nys, nyn 11099 DO k = nzb, nz_do3d 11100 temp_bin = 0.0_wp 11101 DO ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol 11102 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 11103 ENDDO 11104 tend(k,j,i) = temp_bin 11105 ENDDO 11106 ENDDO 11107 ENDDO 11108 IF ( .NOT. mask_surface(mid) ) THEN 11109 DO i = 1, mask_size_l(mid,1) 11110 DO j = 1, mask_size_l(mid,2) 11111 DO k = 1, mask_size_l(mid,3) 11112 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11113 ENDDO 11114 ENDDO 11115 ENDDO 11116 ELSE 11117 DO i = 1, mask_size_l(mid,1) 11118 DO j = 1, mask_size_l(mid,2) 11119 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11120 grid ) 11121 DO k = 1, mask_size_l(mid,3) 11122 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11123 mask_j(mid,j), mask_i(mid,i) ) 11124 ENDDO 11125 ENDDO 11126 ENDDO 11127 ENDIF 11128 resorted = .TRUE. 11129 ELSE 11130 to_be_resorted => s_h2o_av 11131 ENDIF 11132 10964 11133 CASE DEFAULT 10965 11134 found = .FALSE. … … 10999 11168 END SUBROUTINE salsa_data_output_mask 11000 11169 11170 !------------------------------------------------------------------------------! 11171 ! Description: 11172 ! ------------ 11173 !> Creates index tables for different (aerosol) components 11174 !------------------------------------------------------------------------------! 11175 SUBROUTINE component_index_constructor( self, ncomp, nlist, listcomp ) 11176 11177 IMPLICIT NONE 11178 11179 INTEGER(iwp) :: ii !< 11180 INTEGER(iwp) :: jj !< 11181 11182 INTEGER(iwp), INTENT(in) :: nlist ! < Maximum number of components 11183 11184 INTEGER(iwp), INTENT(inout) :: ncomp !< Number of components 11185 11186 CHARACTER(len=3), INTENT(in) :: listcomp(nlist) !< List cof component names 11187 11188 TYPE(component_index), INTENT(inout) :: self !< Object containing the indices of different 11189 !< aerosol components 11190 11191 ncomp = 0 11192 11193 DO WHILE ( listcomp(ncomp+1) /= ' ' .AND. ncomp < nlist ) 11194 ncomp = ncomp + 1 11195 ENDDO 11196 11197 self%ncomp = ncomp 11198 ALLOCATE( self%ind(ncomp), self%comp(ncomp) ) 11199 11200 DO ii = 1, ncomp 11201 self%ind(ii) = ii 11202 ENDDO 11203 11204 jj = 1 11205 DO ii = 1, nlist 11206 IF ( listcomp(ii) == '') CYCLE 11207 self%comp(jj) = listcomp(ii) 11208 jj = jj + 1 11209 ENDDO 11210 11211 END SUBROUTINE component_index_constructor 11212 11213 !------------------------------------------------------------------------------! 11214 ! Description: 11215 ! ------------ 11216 !> Gives the index of a component in the component list 11217 !------------------------------------------------------------------------------! 11218 INTEGER FUNCTION get_index( self, incomp ) 11219 11220 IMPLICIT NONE 11221 11222 CHARACTER(len=*), INTENT(in) :: incomp !< Component name 11223 11224 INTEGER(iwp) :: ii !< index 11225 11226 TYPE(component_index), INTENT(in) :: self !< Object containing the indices of different 11227 !< aerosol components 11228 IF ( ANY( self%comp == incomp ) ) THEN 11229 ii = 1 11230 DO WHILE ( (self%comp(ii) /= incomp) ) 11231 ii = ii + 1 11232 ENDDO 11233 get_index = ii 11234 ELSEIF ( incomp == 'H2O' ) THEN 11235 get_index = self%ncomp + 1 11236 ELSE 11237 WRITE( message_string, * ) 'Incorrect component name given!' 11238 CALL message( 'get_index', 'PA0591', 1, 2, 0, 6, 0 ) 11239 ENDIF 11240 11241 END FUNCTION get_index 11242 11243 !------------------------------------------------------------------------------! 11244 ! Description: 11245 ! ------------ 11246 !> Tells if the (aerosol) component is being used in the simulation 11247 !------------------------------------------------------------------------------! 11248 LOGICAL FUNCTION is_used( self, icomp ) 11249 11250 IMPLICIT NONE 11251 11252 CHARACTER(len=*), INTENT(in) :: icomp !< Component name 11253 11254 TYPE(component_index), INTENT(in) :: self !< Object containing the indices of different 11255 !< aerosol components 11256 11257 IF ( ANY(self%comp == icomp) ) THEN 11258 is_used = .TRUE. 11259 ELSE 11260 is_used = .FALSE. 11261 ENDIF 11262 11263 END FUNCTION 11264 11001 11265 END MODULE salsa_mod
Note: See TracChangeset
for help on using the changeset viewer.
