Changeset 4131 for palm/trunk/SOURCE/salsa_mod.f90
- Timestamp:
- Aug 2, 2019 11:06:18 AM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/salsa_mod.f90 (modified) (100 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/salsa_mod.f90
r4118 r4131 26 26 ! ----------------- 27 27 ! $Id$ 28 ! - Add "salsa_" before each salsa output variable 29 ! - Add a possibility to output the number (salsa_N_UFP) and mass concentration 30 ! (salsa_PM0.1) of ultrafine particles, i.e. particles with a diameter smaller 31 ! than 100 nm 32 ! - Implement aerosol emission mode "parameterized" which is based on the street 33 ! type (similar to the chemistry module). 34 ! - Remove unnecessary nucleation subroutines. 35 ! - Add the z-dimension for gaseous emissions to correspond the implementation 36 ! in the chemistry module 37 ! 38 ! 4118 2019-07-25 16:11:45Z suehring 28 39 ! - When Dirichlet condition is applied in decycling, the boundary conditions are 29 40 ! only set at the ghost points and not at the prognostic grid points as done … … 211 222 bc_lr, bc_lr_cyc, bc_ns, bc_ns_cyc, bc_radiation_l, bc_radiation_n, bc_radiation_r, & 212 223 bc_radiation_s, coupling_char, debug_output, dt_3d, intermediate_timestep_count, & 213 intermediate_timestep_count_max, land_surface, message_string, monotonic_limiter_z, & 214 plant_canopy, pt_surface, salsa, scalar_advec, surface_pressure, & 215 time_since_reference_point, timestep_scheme, tsc, urban_surface, ws_scheme_sca 224 intermediate_timestep_count_max, land_surface, max_pr_salsa, message_string, & 225 monotonic_limiter_z, plant_canopy, pt_surface, salsa, scalar_advec, & 226 surface_pressure, time_since_reference_point, timestep_scheme, tsc, urban_surface, & 227 ws_scheme_sca 216 228 217 229 USE indices, & … … 392 404 !< 1 = read vertical profile from an input file 393 405 INTEGER(iwp) :: lod_gas_emissions = 0 !< level of detail of the gaseous emission data 406 INTEGER(iwp) :: main_street_id = 0 !< lower bound of main street IDs (OpenStreetMaps) for parameterized mode 407 INTEGER(iwp) :: max_street_id = 0 !< upper bound of main street IDs (OpenStreetMaps) for parameterized mode 394 408 INTEGER(iwp) :: nbins_aerosol = 1 !< total number of size bins 395 409 INTEGER(iwp) :: ncc = 1 !< number of chemical components used … … 413 427 !< 9 = homomolecular nucleation of H2SO4 and ORG 414 428 !< + heteromolecular nucleation with H2SO4*ORG 429 INTEGER(iwp) :: salsa_pr_count = 0 !< counter for salsa variable profiles 430 INTEGER(iwp) :: side_street_id = 0 !< lower bound of side street IDs (OpenStreetMaps) for parameterized mode 415 431 INTEGER(iwp) :: start_subrange_1a = 1 !< start index for bin subranges: subrange 1a 416 432 INTEGER(iwp) :: start_subrange_2a = 1 !< subrange 2a … … 419 435 INTEGER(iwp), DIMENSION(nreg) :: nbin = (/ 3, 7/) !< Number of size bins per subrange: 1 & 2 420 436 421 INTEGER(iwp), DIMENSION(ngases_salsa) :: gas_index_chem = &422 (/ 1, 1, 1, 1, 1/) !< gas indices in chemistry_model_mod423 !< 1 = H2SO4, 2 = HNO3, 3 = NH3,4 = OCNV, 5 = OCSV437 INTEGER(iwp), DIMENSION(ngases_salsa) :: gas_index_chem = (/ 1, 1, 1, 1, 1/) !< gas indices in chemistry_model_mod 438 !< 1 = H2SO4, 2 = HNO3, 439 !< 3 = NH3, 4 = OCNV, 5 = OCSV 424 440 INTEGER(iwp), DIMENSION(ngases_salsa) :: emission_index_chem !< gas indices in the gas emission file 441 INTEGER(iwp), DIMENSION(99) :: salsa_pr_index = 0 !< index for salsa profiles 425 442 426 443 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: k_topo_top !< vertical index of the topography top … … 466 483 LOGICAL :: lspartition = .FALSE. !< Partition of HNO3 and NH3 467 484 468 REAL(wp) :: act_coeff = 1.0E-7_wp !< Activation coefficient 485 REAL(wp) :: act_coeff = 1.0E-7_wp !< Activation coefficient (1/s) 469 486 REAL(wp) :: dt_salsa = 0.00001_wp !< Time step of SALSA 487 REAL(wp) :: emiss_factor_main = 0.0_wp !< relative emission factor for main streets 488 REAL(wp) :: emiss_factor_side = 0.0_wp !< relative emission factor for side streets 470 489 REAL(wp) :: h2so4_init = nclim !< Init value for sulphuric acid gas 471 490 REAL(wp) :: hno3_init = nclim !< Init value for nitric acid gas … … 557 576 INTEGER(iwp) :: tind !< time index for reference time in salsa emission data 558 577 559 INTEGER(iwp), DIMENSION(maxspec) :: cc_in put_to_model!<578 INTEGER(iwp), DIMENSION(maxspec) :: cc_in2mod = 0 !< 560 579 561 580 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: cat_index !< Index of emission categories … … 630 649 TYPE salsa_variable 631 650 632 REAL(wp), ALLOCATABLE, DIMENSION(:):: init !<633 634 REAL(wp), ALLOCATABLE, DIMENSION(:,:):: diss_s !<635 REAL(wp), ALLOCATABLE, DIMENSION(:,:):: flux_s !<636 REAL(wp), ALLOCATABLE, DIMENSION(:,:):: source !<637 REAL(wp), ALLOCATABLE, DIMENSION(:,:):: sums_ws_l !<638 639 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:):: diss_l !<640 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:):: flux_l !<641 642 REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS :: conc !<643 REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS :: conc_p !<644 REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS :: tconc_m !<651 REAL(wp), DIMENSION(:), ALLOCATABLE :: init !< 652 653 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: diss_s !< 654 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: flux_s !< 655 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: source !< 656 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sums_ws_l !< 657 658 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: diss_l !< 659 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: flux_l !< 660 661 REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS :: conc !< 662 REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS :: conc_p !< 663 REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS :: tconc_m !< 645 664 646 665 END TYPE salsa_variable … … 673 692 674 693 TYPE(chem_emis_att_type) :: chem_emission_att !< chemistry emission attributes 675 TYPE(chem_emis_val_type) :: chem_emission !< chemistry emission values 694 695 TYPE(chem_emis_val_type), DIMENSION(:), ALLOCATABLE :: chem_emission !< chemistry emissions 676 696 677 697 TYPE(t_section), DIMENSION(:), ALLOCATABLE :: aero !< local aerosol properties … … 683 703 TYPE(match_surface), DIMENSION(0:3) :: usm_to_depo_v !< to match the deposition mod. and vertical USM surfaces 684 704 ! 685 !-- SALSA variables: as x = x(k,j,i,bin). 705 !-- SALSA variables: as x = x(k,j,i,bin). 686 706 !-- The 4th dimension contains all the size bins sequentially for each aerosol species + water. 687 707 ! … … 689 709 ! 690 710 !-- Number concentration (#/m3) 691 TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET :: aerosol_number !<692 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: nconc_1 !<693 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: nconc_2 !<694 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: nconc_3 !<711 TYPE(salsa_variable), DIMENSION(:), ALLOCATABLE, TARGET :: aerosol_number !< 712 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: nconc_1 !< 713 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: nconc_2 !< 714 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: nconc_3 !< 695 715 ! 696 716 !-- Mass concentration (kg/m3) 697 TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET :: aerosol_mass !<698 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mconc_1 !<699 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mconc_2 !<700 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mconc_3 !<717 TYPE(salsa_variable), DIMENSION(:), ALLOCATABLE, TARGET :: aerosol_mass !< 718 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: mconc_1 !< 719 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: mconc_2 !< 720 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: mconc_3 !< 701 721 ! 702 722 !-- Gaseous concentrations (#/m3) 703 TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET :: salsa_gas !<704 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: gconc_1 !<705 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: gconc_2 !<706 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: gconc_3 !<723 TYPE(salsa_variable), DIMENSION(:), ALLOCATABLE, TARGET :: salsa_gas !< 724 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: gconc_1 !< 725 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: gconc_2 !< 726 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: gconc_3 !< 707 727 ! 708 728 !-- Diagnostic tracers 709 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:):: sedim_vd !< sedimentation velocity per bin (m/s)710 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:):: ra_dry !< aerosol dry radius (m)729 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: sedim_vd !< sedimentation velocity per bin (m/s) 730 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ra_dry !< aerosol dry radius (m) 711 731 712 732 !-- Particle component index tables … … 717 737 ! 718 738 !-- Gases: 719 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_h2so4_av !< H2SO4720 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_hno3_av !< HNO3721 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_nh3_av !< NH3722 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_ocnv_av !< non-volatile OC723 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_ocsv_av !< semi-volatile OC739 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: g_h2so4_av !< H2SO4 740 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: g_hno3_av !< HNO3 741 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: g_nh3_av !< NH3 742 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: g_ocnv_av !< non-volatile OC 743 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: g_ocsv_av !< semi-volatile OC 724 744 ! 725 745 !-- Integrated: 726 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: ldsa_av !< lung-deposited surface area 727 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: ntot_av !< total number concentration 728 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: pm25_av !< PM2.5 729 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: pm10_av !< PM10 746 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ldsa_av !< lung-deposited surface area 747 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ntot_av !< total number concentration 748 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: nufp_av !< ultrafine particles (UFP) 749 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: pm01_av !< PM0.1 750 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: pm25_av !< PM2.5 751 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: pm10_av !< PM10 730 752 ! 731 753 !-- In the particle phase: 732 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_bc_av !< black carbon733 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_du_av !< dust734 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_h2o_av !< liquid water735 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_nh_av !< ammonia736 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_no_av !< nitrates737 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_oc_av !< org. carbon738 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_so4_av !< sulphates739 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_ss_av !< sea salt754 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: s_bc_av !< black carbon 755 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: s_du_av !< dust 756 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: s_h2o_av !< liquid water 757 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: s_nh_av !< ammonia 758 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: s_no_av !< nitrates 759 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: s_oc_av !< org. carbon 760 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: s_so4_av !< sulphates 761 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: s_ss_av !< sea salt 740 762 ! 741 763 !-- Bin specific mass and number concentrations: 742 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mbins_av !< bin mas743 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: nbins_av !< bin number764 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: mbins_av !< bin mas 765 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: nbins_av !< bin number 744 766 745 767 ! 746 768 !-- PALM interfaces: 747 ! 748 !-- Boundary conditions: 769 770 INTERFACE salsa_actions 771 MODULE PROCEDURE salsa_actions 772 MODULE PROCEDURE salsa_actions_ij 773 END INTERFACE salsa_actions 774 775 INTERFACE salsa_3d_data_averaging 776 MODULE PROCEDURE salsa_3d_data_averaging 777 END INTERFACE salsa_3d_data_averaging 778 749 779 INTERFACE salsa_boundary_conds 750 780 MODULE PROCEDURE salsa_boundary_conds 751 781 MODULE PROCEDURE salsa_boundary_conds_decycle 752 782 END INTERFACE salsa_boundary_conds 753 ! 754 !-- Data output checks for 2D/3D data to be done in check_parameters 783 755 784 INTERFACE salsa_check_data_output 756 785 MODULE PROCEDURE salsa_check_data_output 757 786 END INTERFACE salsa_check_data_output 758 ! 759 !-- Input parameter checks to be done in check_parameters 787 788 INTERFACE salsa_check_data_output_pr 789 MODULE PROCEDURE salsa_check_data_output_pr 790 END INTERFACE salsa_check_data_output_pr 791 760 792 INTERFACE salsa_check_parameters 761 793 MODULE PROCEDURE salsa_check_parameters 762 794 END INTERFACE salsa_check_parameters 763 ! 764 !-- Averaging of 3D data for output 765 INTERFACE salsa_3d_data_averaging 766 MODULE PROCEDURE salsa_3d_data_averaging 767 END INTERFACE salsa_3d_data_averaging 768 ! 769 !-- Data output of 2D quantities 795 770 796 INTERFACE salsa_data_output_2d 771 797 MODULE PROCEDURE salsa_data_output_2d 772 798 END INTERFACE salsa_data_output_2d 773 ! 774 !-- Data output of 3D data 799 775 800 INTERFACE salsa_data_output_3d 776 801 MODULE PROCEDURE salsa_data_output_3d 777 802 END INTERFACE salsa_data_output_3d 778 ! 779 !-- Data output of 3D data 803 780 804 INTERFACE salsa_data_output_mask 781 805 MODULE PROCEDURE salsa_data_output_mask 782 806 END INTERFACE salsa_data_output_mask 783 ! 784 !-- Definition of data output quantities 807 785 808 INTERFACE salsa_define_netcdf_grid 786 809 MODULE PROCEDURE salsa_define_netcdf_grid 787 810 END INTERFACE salsa_define_netcdf_grid 788 ! 789 !-- Output of information to the header file 811 812 INTERFACE salsa_driver 813 MODULE PROCEDURE salsa_driver 814 END INTERFACE salsa_driver 815 816 INTERFACE salsa_emission_update 817 MODULE PROCEDURE salsa_emission_update 818 END INTERFACE salsa_emission_update 819 820 INTERFACE salsa_exchange_horiz_bounds 821 MODULE PROCEDURE salsa_exchange_horiz_bounds 822 END INTERFACE salsa_exchange_horiz_bounds 823 790 824 INTERFACE salsa_header 791 825 MODULE PROCEDURE salsa_header 792 826 END INTERFACE salsa_header 793 ! 794 !-- Initialization actions 827 795 828 INTERFACE salsa_init 796 829 MODULE PROCEDURE salsa_init 797 830 END INTERFACE salsa_init 798 ! 799 !-- Initialization of arrays 831 800 832 INTERFACE salsa_init_arrays 801 833 MODULE PROCEDURE salsa_init_arrays 802 834 END INTERFACE salsa_init_arrays 803 ! 804 !-- Writing of binary output for restart runs !!! renaming?! 805 INTERFACE salsa_wrd_local 806 MODULE PROCEDURE salsa_wrd_local 807 END INTERFACE salsa_wrd_local 808 ! 809 !-- Reading of NAMELIST parameters 810 INTERFACE salsa_parin 811 MODULE PROCEDURE salsa_parin 812 END INTERFACE salsa_parin 813 ! 814 !-- Reading of parameters for restart runs 815 INTERFACE salsa_rrd_local 816 MODULE PROCEDURE salsa_rrd_local 817 END INTERFACE salsa_rrd_local 818 ! 819 !-- Swapping of time levels (required for prognostic variables) 820 INTERFACE salsa_swap_timelevel 821 MODULE PROCEDURE salsa_swap_timelevel 822 END INTERFACE salsa_swap_timelevel 823 ! 824 !-- Interface between PALM and salsa 825 INTERFACE salsa_driver 826 MODULE PROCEDURE salsa_driver 827 END INTERFACE salsa_driver 828 829 !-- Actions salsa variables 830 INTERFACE salsa_actions 831 MODULE PROCEDURE salsa_actions 832 MODULE PROCEDURE salsa_actions_ij 833 END INTERFACE salsa_actions 834 ! 835 !-- Non-advective processes (i.e. aerosol dynamic reactions) for salsa variables 835 836 836 INTERFACE salsa_non_advective_processes 837 837 MODULE PROCEDURE salsa_non_advective_processes 838 838 MODULE PROCEDURE salsa_non_advective_processes_ij 839 839 END INTERFACE salsa_non_advective_processes 840 ! 841 !-- Exchange horiz for salsa variables 842 INTERFACE salsa_exchange_horiz_bounds 843 MODULE PROCEDURE salsa_exchange_horiz_bounds 844 END INTERFACE salsa_exchange_horiz_bounds 845 ! 846 !-- Prognostics equations for salsa variables 840 841 INTERFACE salsa_parin 842 MODULE PROCEDURE salsa_parin 843 END INTERFACE salsa_parin 844 847 845 INTERFACE salsa_prognostic_equations 848 846 MODULE PROCEDURE salsa_prognostic_equations 849 847 MODULE PROCEDURE salsa_prognostic_equations_ij 850 848 END INTERFACE salsa_prognostic_equations 851 ! 852 !-- Tendency salsa variables 849 850 INTERFACE salsa_rrd_local 851 MODULE PROCEDURE salsa_rrd_local 852 END INTERFACE salsa_rrd_local 853 854 INTERFACE salsa_statistics 855 MODULE PROCEDURE salsa_statistics 856 END INTERFACE salsa_statistics 857 858 INTERFACE salsa_swap_timelevel 859 MODULE PROCEDURE salsa_swap_timelevel 860 END INTERFACE salsa_swap_timelevel 861 853 862 INTERFACE salsa_tendency 854 863 MODULE PROCEDURE salsa_tendency 855 864 MODULE PROCEDURE salsa_tendency_ij 856 865 END INTERFACE salsa_tendency 866 867 INTERFACE salsa_wrd_local 868 MODULE PROCEDURE salsa_wrd_local 869 END INTERFACE salsa_wrd_local 857 870 858 871 … … 867 880 salsa_emission_update, salsa_header, salsa_init, salsa_init_arrays, salsa_parin, & 868 881 salsa_rrd_local, salsa_swap_timelevel, salsa_prognostic_equations, salsa_wrd_local, & 869 salsa_actions, salsa_non_advective_processes, salsa_exchange_horiz_bounds 882 salsa_actions, salsa_non_advective_processes, salsa_exchange_horiz_bounds, & 883 salsa_check_data_output_pr, salsa_statistics 870 884 ! 871 885 !-- Public parameters, constants and initial values … … 888 902 SUBROUTINE salsa_parin 889 903 904 USE control_parameters, & 905 ONLY: data_output_pr 906 890 907 IMPLICIT NONE 891 908 892 CHARACTER(LEN=80) :: line !< dummy string that contains the current line 893 !< of the parameter file 894 895 NAMELIST /salsa_parameters/ aerosol_flux_dpg, aerosol_flux_mass_fracs_a, & 896 aerosol_flux_mass_fracs_b, aerosol_flux_sigmag, & 897 advect_particle_water, bc_salsa_b, bc_salsa_t, & 898 decycle_salsa_lr, decycle_method_salsa, decycle_salsa_ns, & 899 depo_pcm_par, depo_pcm_type, & 900 depo_surf_par, dpg, dt_salsa, feedback_to_palm, h2so4_init, & 901 hno3_init, init_gases_type, init_aerosol_type, listspec, & 902 mass_fracs_a, mass_fracs_b, n_lognorm, nbin, nest_salsa, nf2a,& 903 nh3_init, nj3, nlcnd, nlcndgas, nlcndh2oae, nlcoag, nldepo, & 904 nldepo_pcm, nldepo_surf, nldistupdate, nsnucl, ocnv_init, & 905 ocsv_init, read_restart_data_salsa, reglim, salsa, & 906 salsa_emission_mode, sigmag, skip_time_do_salsa, & 907 surface_aerosol_flux, van_der_waals_coagc, write_binary_salsa 909 CHARACTER(LEN=80) :: line !< dummy string that contains the current line of parameter file 910 911 INTEGER(iwp) :: i !< loop index 912 INTEGER(iwp) :: max_pr_salsa_tmp !< dummy variable 913 914 NAMELIST /salsa_parameters/ aerosol_flux_dpg, & 915 aerosol_flux_mass_fracs_a, & 916 aerosol_flux_mass_fracs_b, & 917 aerosol_flux_sigmag, & 918 advect_particle_water, & 919 bc_salsa_b, & 920 bc_salsa_t, & 921 decycle_salsa_lr, & 922 decycle_method_salsa, & 923 decycle_salsa_ns, & 924 depo_pcm_par, & 925 depo_pcm_type, & 926 depo_surf_par, & 927 dpg, & 928 dt_salsa, & 929 emiss_factor_main, & 930 emiss_factor_side, & 931 feedback_to_palm, & 932 h2so4_init, & 933 hno3_init, & 934 init_gases_type, & 935 init_aerosol_type, & 936 listspec, & 937 main_street_id, & 938 mass_fracs_a, & 939 mass_fracs_b, & 940 max_street_id, & 941 n_lognorm, & 942 nbin, & 943 nest_salsa, & 944 nf2a, & 945 nh3_init, & 946 nj3, & 947 nlcnd, & 948 nlcndgas, & 949 nlcndh2oae, & 950 nlcoag, & 951 nldepo, & 952 nldepo_pcm, & 953 nldepo_surf, & 954 nldistupdate, & 955 nsnucl, & 956 ocnv_init, & 957 ocsv_init, & 958 read_restart_data_salsa, & 959 reglim, & 960 salsa, & 961 salsa_emission_mode, & 962 sigmag, & 963 side_street_id, & 964 skip_time_do_salsa, & 965 surface_aerosol_flux, & 966 van_der_waals_coagc, & 967 write_binary_salsa 908 968 909 969 line = ' ' … … 924 984 925 985 10 CONTINUE 986 ! 987 !-- Update the number of output profiles 988 max_pr_salsa_tmp = 0 989 i = 1 990 DO WHILE ( data_output_pr(i) /= ' ' .AND. i <= 100 ) 991 IF ( TRIM( data_output_pr(i)(1:6) ) == 'salsa_' ) max_pr_salsa_tmp = max_pr_salsa_tmp + 1 992 i = i + 1 993 ENDDO 994 IF ( max_pr_salsa_tmp > 0 ) max_pr_salsa = max_pr_salsa_tmp 926 995 927 996 END SUBROUTINE salsa_parin … … 1003 1072 !-- Check for the grid 1004 1073 1005 IF ( var(1:2) == 'g_' ) THEN 1006 grid_x = 'x' 1007 grid_y = 'y' 1008 grid_z = 'zu' 1009 ELSEIF ( var(1:4) == 'LDSA' ) THEN 1010 grid_x = 'x' 1011 grid_y = 'y' 1012 grid_z = 'zu' 1013 ELSEIF ( var(1:5) == 'm_bin' ) THEN 1014 grid_x = 'x' 1015 grid_y = 'y' 1016 grid_z = 'zu' 1017 ELSEIF ( var(1:5) == 'N_bin' ) THEN 1018 grid_x = 'x' 1019 grid_y = 'y' 1020 grid_z = 'zu' 1021 ELSEIF ( var(1:4) == 'Ntot' ) THEN 1022 grid_x = 'x' 1023 grid_y = 'y' 1024 grid_z = 'zu' 1025 ELSEIF ( var(1:2) == 'PM' ) THEN 1026 grid_x = 'x' 1027 grid_y = 'y' 1028 grid_z = 'zu' 1029 ELSEIF ( var(1:2) == 's_' ) THEN 1074 IF ( var(1:6) == 'salsa_' ) THEN ! same grid for all salsa output variables 1030 1075 grid_x = 'x' 1031 1076 grid_y = 'y' … … 1163 1208 INTEGER(iwp) :: gases_available !< Number of available gas components in the chemistry model 1164 1209 INTEGER(iwp) :: i !< loop index for allocating 1210 INTEGER(iwp) :: ii !< index for indexing chemical components 1165 1211 INTEGER(iwp) :: l !< loop index for allocating: surfaces 1166 1212 INTEGER(iwp) :: lsp !< loop index for chem species in the chemistry model … … 1194 1240 1195 1241 ncomponents_mass = ncc 1196 IF ( advect_particle_water ) ncomponents_mass = ncc + 1 ! Add water 1197 1198 ! 1199 !-- Allocate: 1200 ALLOCATE( aero(nbins_aerosol), bc_am_t_val(nbins_aerosol*ncomponents_mass), & 1201 bc_an_t_val(nbins_aerosol), bc_gt_t_val(ngases_salsa), bin_low_limits(nbins_aerosol),& 1202 nsect(nbins_aerosol), massacc(nbins_aerosol) ) 1203 ALLOCATE( k_topo_top(nysg:nyng,nxlg:nxrg) ) 1204 IF ( nldepo ) ALLOCATE( sedim_vd(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 1205 ALLOCATE( ra_dry(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 1206 1207 ! 1208 !-- Aerosol number concentration 1209 ALLOCATE( aerosol_number(nbins_aerosol) ) 1210 ALLOCATE( nconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol), & 1211 nconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol), & 1212 nconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 1213 nconc_1 = 0.0_wp 1214 nconc_2 = 0.0_wp 1215 nconc_3 = 0.0_wp 1216 1217 DO i = 1, nbins_aerosol 1218 aerosol_number(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_1(:,:,:,i) 1219 aerosol_number(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_2(:,:,:,i) 1220 aerosol_number(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_3(:,:,:,i) 1221 ALLOCATE( aerosol_number(i)%flux_s(nzb+1:nzt,0:threads_per_task-1), & 1222 aerosol_number(i)%diss_s(nzb+1:nzt,0:threads_per_task-1), & 1223 aerosol_number(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 1224 aerosol_number(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 1225 aerosol_number(i)%init(nzb:nzt+1), & 1226 aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 1227 aerosol_number(i)%init = nclim 1228 IF ( include_emission .OR. ( nldepo .AND. nldepo_surf ) ) THEN 1229 ALLOCATE( aerosol_number(i)%source(nys:nyn,nxl:nxr) ) 1230 ENDIF 1231 ENDDO 1232 1233 ! 1234 !-- Aerosol mass concentration 1235 ALLOCATE( aerosol_mass(ncomponents_mass*nbins_aerosol) ) 1236 ALLOCATE( mconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol), & 1237 mconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol), & 1238 mconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol) ) 1239 mconc_1 = 0.0_wp 1240 mconc_2 = 0.0_wp 1241 mconc_3 = 0.0_wp 1242 1243 DO i = 1, ncomponents_mass*nbins_aerosol 1244 aerosol_mass(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_1(:,:,:,i) 1245 aerosol_mass(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_2(:,:,:,i) 1246 aerosol_mass(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_3(:,:,:,i) 1247 ALLOCATE( aerosol_mass(i)%flux_s(nzb+1:nzt,0:threads_per_task-1), & 1248 aerosol_mass(i)%diss_s(nzb+1:nzt,0:threads_per_task-1), & 1249 aerosol_mass(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 1250 aerosol_mass(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 1251 aerosol_mass(i)%init(nzb:nzt+1), & 1252 aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 1253 aerosol_mass(i)%init = mclim 1254 IF ( include_emission .OR. ( nldepo .AND. nldepo_surf ) ) THEN 1255 ALLOCATE( aerosol_mass(i)%source(nys:nyn,nxl:nxr) ) 1256 ENDIF 1257 ENDDO 1258 1259 ! 1260 !-- Surface fluxes: answs = aerosol number, amsws = aerosol mass 1261 ! 1262 !-- Horizontal surfaces: default type 1263 DO l = 0, 2 ! upward (l=0), downward (l=1) and model top (l=2) 1264 ALLOCATE( surf_def_h(l)%answs( 1:surf_def_h(l)%ns, nbins_aerosol ) ) 1265 ALLOCATE( surf_def_h(l)%amsws( 1:surf_def_h(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1266 surf_def_h(l)%answs = 0.0_wp 1267 surf_def_h(l)%amsws = 0.0_wp 1268 ENDDO 1269 ! 1270 !-- Horizontal surfaces: natural type 1271 ALLOCATE( surf_lsm_h%answs( 1:surf_lsm_h%ns, nbins_aerosol ) ) 1272 ALLOCATE( surf_lsm_h%amsws( 1:surf_lsm_h%ns, nbins_aerosol*ncomponents_mass ) ) 1273 surf_lsm_h%answs = 0.0_wp 1274 surf_lsm_h%amsws = 0.0_wp 1275 ! 1276 !-- Horizontal surfaces: urban type 1277 ALLOCATE( surf_usm_h%answs( 1:surf_usm_h%ns, nbins_aerosol ) ) 1278 ALLOCATE( surf_usm_h%amsws( 1:surf_usm_h%ns, nbins_aerosol*ncomponents_mass ) ) 1279 surf_usm_h%answs = 0.0_wp 1280 surf_usm_h%amsws = 0.0_wp 1281 1282 ! 1283 !-- Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and westward (l=3) facing 1284 DO l = 0, 3 1285 ALLOCATE( surf_def_v(l)%answs( 1:surf_def_v(l)%ns, nbins_aerosol ) ) 1286 surf_def_v(l)%answs = 0.0_wp 1287 ALLOCATE( surf_def_v(l)%amsws( 1:surf_def_v(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1288 surf_def_v(l)%amsws = 0.0_wp 1289 1290 ALLOCATE( surf_lsm_v(l)%answs( 1:surf_lsm_v(l)%ns, nbins_aerosol ) ) 1291 surf_lsm_v(l)%answs = 0.0_wp 1292 ALLOCATE( surf_lsm_v(l)%amsws( 1:surf_lsm_v(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1293 surf_lsm_v(l)%amsws = 0.0_wp 1294 1295 ALLOCATE( surf_usm_v(l)%answs( 1:surf_usm_v(l)%ns, nbins_aerosol ) ) 1296 surf_usm_v(l)%answs = 0.0_wp 1297 ALLOCATE( surf_usm_v(l)%amsws( 1:surf_usm_v(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1298 surf_usm_v(l)%amsws = 0.0_wp 1299 1300 ENDDO 1301 1302 ! 1303 !-- Concentration of gaseous tracers (1. SO4, 2. HNO3, 3. NH3, 4. OCNV, 5. OCSV) 1304 !-- (number concentration (#/m3) ) 1305 ! 1306 !-- If chemistry is on, read gas phase concentrations from there. Otherwise, 1307 !-- allocate salsa_gas array. 1308 1309 IF ( air_chemistry ) THEN 1310 DO lsp = 1, nvar 1311 SELECT CASE ( TRIM( chem_species(lsp)%name ) ) 1312 CASE ( 'H2SO4', 'h2so4' ) 1313 gases_available = gases_available + 1 1314 gas_index_chem(1) = lsp 1315 CASE ( 'HNO3', 'hno3' ) 1316 gases_available = gases_available + 1 1317 gas_index_chem(2) = lsp 1318 CASE ( 'NH3', 'nh3' ) 1319 gases_available = gases_available + 1 1320 gas_index_chem(3) = lsp 1321 CASE ( 'OCNV', 'ocnv' ) 1322 gases_available = gases_available + 1 1323 gas_index_chem(4) = lsp 1324 CASE ( 'OCSV', 'ocsv' ) 1325 gases_available = gases_available + 1 1326 gas_index_chem(5) = lsp 1327 END SELECT 1328 ENDDO 1329 1330 IF ( gases_available == ngases_salsa ) THEN 1331 salsa_gases_from_chem = .TRUE. 1332 ELSE 1333 WRITE( message_string, * ) 'SALSA is run together with chemistry but not all gaseous '// & 1334 'components are provided by kpp (H2SO4, HNO3, NH3, OCNV, OCSV)' 1335 CALL message( 'check_parameters', 'PA0599', 1, 2, 0, 6, 0 ) 1336 ENDIF 1337 1338 ELSE 1339 1340 ALLOCATE( salsa_gas(ngases_salsa) ) 1341 ALLOCATE( gconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa), & 1342 gconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa), & 1343 gconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) ) 1344 gconc_1 = 0.0_wp 1345 gconc_2 = 0.0_wp 1346 gconc_3 = 0.0_wp 1347 1348 DO i = 1, ngases_salsa 1349 salsa_gas(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_1(:,:,:,i) 1350 salsa_gas(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_2(:,:,:,i) 1351 salsa_gas(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_3(:,:,:,i) 1352 ALLOCATE( salsa_gas(i)%flux_s(nzb+1:nzt,0:threads_per_task-1), & 1353 salsa_gas(i)%diss_s(nzb+1:nzt,0:threads_per_task-1), & 1354 salsa_gas(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),& 1355 salsa_gas(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),& 1356 salsa_gas(i)%init(nzb:nzt+1), & 1357 salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 1358 salsa_gas(i)%init = nclim 1359 IF ( include_emission ) ALLOCATE( salsa_gas(i)%source(nys:nys,nxl:nxr) ) 1360 ENDDO 1361 ! 1362 !-- Surface fluxes: gtsws = gaseous tracer flux 1363 ! 1364 !-- Horizontal surfaces: default type 1365 DO l = 0, 2 ! upward (l=0), downward (l=1) and model top (l=2) 1366 ALLOCATE( surf_def_h(l)%gtsws( 1:surf_def_h(l)%ns, ngases_salsa ) ) 1367 surf_def_h(l)%gtsws = 0.0_wp 1368 ENDDO 1369 !-- Horizontal surfaces: natural type 1370 ALLOCATE( surf_lsm_h%gtsws( 1:surf_lsm_h%ns, ngases_salsa ) ) 1371 surf_lsm_h%gtsws = 0.0_wp 1372 !-- Horizontal surfaces: urban type 1373 ALLOCATE( surf_usm_h%gtsws( 1:surf_usm_h%ns, ngases_salsa ) ) 1374 surf_usm_h%gtsws = 0.0_wp 1375 ! 1376 !-- Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and 1377 !-- westward (l=3) facing 1378 DO l = 0, 3 1379 ALLOCATE( surf_def_v(l)%gtsws( 1:surf_def_v(l)%ns, ngases_salsa ) ) 1380 surf_def_v(l)%gtsws = 0.0_wp 1381 ALLOCATE( surf_lsm_v(l)%gtsws( 1:surf_lsm_v(l)%ns, ngases_salsa ) ) 1382 surf_lsm_v(l)%gtsws = 0.0_wp 1383 ALLOCATE( surf_usm_v(l)%gtsws( 1:surf_usm_v(l)%ns, ngases_salsa ) ) 1384 surf_usm_v(l)%gtsws = 0.0_wp 1385 ENDDO 1386 ENDIF 1387 1388 IF ( ws_scheme_sca ) THEN 1389 1390 IF ( salsa ) THEN 1391 ALLOCATE( sums_salsa_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 1392 sums_salsa_ws_l = 0.0_wp 1393 ENDIF 1394 1395 ENDIF 1396 ! 1397 !-- Set control flags for decycling only at lateral boundary cores. Within the inner cores the 1398 !-- decycle flag is set to .FALSE.. Even though it does not affect the setting of chemistry boundary 1399 !-- conditions, this flag is used to set advection control flags appropriately. 1400 decycle_salsa_lr = MERGE( decycle_salsa_lr, .FALSE., nxl == 0 .OR. nxr == nx ) 1401 decycle_salsa_ns = MERGE( decycle_salsa_ns, .FALSE., nys == 0 .OR. nyn == ny ) 1402 ! 1403 !-- Decycling can be applied separately for aerosol variables, while wind and other scalars may have 1404 !-- cyclic or nested boundary conditions. However, large gradients near the boundaries may produce 1405 !-- stationary numerical oscillations near the lateral boundaries when a higher-order scheme is 1406 !-- applied near these boundaries. To get rid-off this, set-up additional flags that control the 1407 !-- order of the scalar advection scheme near the lateral boundaries for passive scalars with 1408 !-- decycling. 1409 IF ( scalar_advec == 'ws-scheme' ) THEN 1410 ALLOCATE( salsa_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1411 ! 1412 !-- In case of decycling, set Neuman boundary conditions for wall_flags_0 bit 31 instead of 1413 !-- cyclic boundary conditions. Bit 31 is used to identify extended degradation zones (please see 1414 !-- the following comment). Note, since several also other modules may access this bit but may 1415 !-- have other boundary conditions, the original value of wall_flags_0 bit 31 must not be 1416 !-- modified. Hence, store the boundary conditions directly on salsa_advc_flags_s. 1417 !-- salsa_advc_flags_s will be later overwritten in ws_init_flags_scalar and bit 31 won't be used 1418 !-- to control the numerical order. 1419 !-- Initialize with flag 31 only. 1420 salsa_advc_flags_s = 0 1421 salsa_advc_flags_s = MERGE( IBSET( salsa_advc_flags_s, 31 ), 0, BTEST( wall_flags_0, 31 ) ) 1422 1423 IF ( decycle_salsa_ns ) THEN 1424 IF ( nys == 0 ) THEN 1425 DO i = 1, nbgp 1426 salsa_advc_flags_s(:,nys-i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nys,:), 31 ), & 1427 IBCLR( salsa_advc_flags_s(:,nys,:), 31 ), & 1428 BTEST( salsa_advc_flags_s(:,nys,:), 31 ) ) 1429 ENDDO 1430 ENDIF 1431 IF ( nyn == ny ) THEN 1432 DO i = 1, nbgp 1433 salsa_advc_flags_s(:,nyn+i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nyn,:), 31 ), & 1434 IBCLR( salsa_advc_flags_s(:,nyn,:), 31 ), & 1435 BTEST( salsa_advc_flags_s(:,nyn,:), 31 ) ) 1436 ENDDO 1437 ENDIF 1438 ENDIF 1439 IF ( decycle_salsa_lr ) THEN 1440 IF ( nxl == 0 ) THEN 1441 DO i = 1, nbgp 1442 salsa_advc_flags_s(:,:,nxl-i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxl), 31 ), & 1443 IBCLR( salsa_advc_flags_s(:,:,nxl), 31 ), & 1444 BTEST( salsa_advc_flags_s(:,:,nxl), 31 ) ) 1445 ENDDO 1446 ENDIF 1447 IF ( nxr == nx ) THEN 1448 DO i = 1, nbgp 1449 salsa_advc_flags_s(:,:,nxr+i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxr), 31 ), & 1450 IBCLR( salsa_advc_flags_s(:,:,nxr), 31 ), & 1451 BTEST( salsa_advc_flags_s(:,:,nxr), 31 ) ) 1452 ENDDO 1453 ENDIF 1454 ENDIF 1455 ! 1456 !-- To initialise the advection flags appropriately, pass the boundary flags to 1457 !-- ws_init_flags_scalar. The last argument in ws_init_flags_scalar indicates that a passive 1458 !-- scalar is being treated and the horizontal advection terms are degraded already 2 grid points 1459 !-- before the lateral boundary. Also, extended degradation zones are applied, where 1460 !-- horizontal advection of scalars is discretised by the first-order scheme at all grid points 1461 !-- in the vicinity of buildings (<= 3 grid points). Even though no building is within the 1462 !-- numerical stencil, the first-order scheme is used. At fourth and fifth grid points, the order 1463 !-- of the horizontal advection scheme is successively upgraded. 1464 !-- These degradations of the advection scheme are done to avoid stationary numerical 1465 !-- oscillations, which are responsible for high concentration maxima that may appear e.g. under 1466 !-- shear-free stable conditions. 1467 CALL ws_init_flags_scalar( bc_dirichlet_l .OR. bc_radiation_l .OR. decycle_salsa_lr, & 1468 bc_dirichlet_n .OR. bc_radiation_n .OR. decycle_salsa_ns, & 1469 bc_dirichlet_r .OR. bc_radiation_r .OR. decycle_salsa_lr, & 1470 bc_dirichlet_s .OR. bc_radiation_s .OR. decycle_salsa_ns, & 1471 salsa_advc_flags_s, .TRUE. ) 1472 ENDIF 1473 1474 1475 END SUBROUTINE salsa_init_arrays 1476 1477 !------------------------------------------------------------------------------! 1478 ! Description: 1479 ! ------------ 1480 !> Initialization of SALSA. Based on salsa_initialize in UCLALES-SALSA. 1481 !> Subroutines salsa_initialize, SALSAinit and DiagInitAero in UCLALES-SALSA are 1482 !> also merged here. 1483 !------------------------------------------------------------------------------! 1484 SUBROUTINE salsa_init 1485 1486 IMPLICIT NONE 1487 1488 INTEGER(iwp) :: i !< 1489 INTEGER(iwp) :: ib !< loop index for aerosol number bins 1490 INTEGER(iwp) :: ic !< loop index for aerosol mass bins 1491 INTEGER(iwp) :: ig !< loop index for gases 1492 INTEGER(iwp) :: ii !< index for indexing 1493 INTEGER(iwp) :: j !< 1494 1495 IF ( debug_output ) CALL debug_message( 'salsa_init', 'start' ) 1496 1497 bin_low_limits = 0.0_wp 1498 k_topo_top = 0 1499 nsect = 0.0_wp 1500 massacc = 1.0_wp 1501 1242 IF ( advect_particle_water ) ncomponents_mass = ncc + 1 ! Add water 1502 1243 ! 1503 1244 !-- Indices for chemical components used (-1 = not used) … … 1538 1279 ENDIF 1539 1280 ! 1540 !-- Partition and dissolutional growth by gaseous HNO3 and NH3 1541 IF ( index_no > 0 .AND. index_nh > 0 .AND. index_so4 > 0 ) lspartition = .TRUE. 1281 !-- Allocate: 1282 ALLOCATE( aero(nbins_aerosol), bc_am_t_val(nbins_aerosol*ncomponents_mass), & 1283 bc_an_t_val(nbins_aerosol), bc_gt_t_val(ngases_salsa), bin_low_limits(nbins_aerosol),& 1284 nsect(nbins_aerosol), massacc(nbins_aerosol) ) 1285 ALLOCATE( k_topo_top(nysg:nyng,nxlg:nxrg) ) 1286 IF ( nldepo ) ALLOCATE( sedim_vd(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 1287 ALLOCATE( ra_dry(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 1288 ! 1289 !-- Initialise the sectional particle size distribution 1290 CALL set_sizebins 1291 ! 1292 !-- Aerosol number concentration 1293 ALLOCATE( aerosol_number(nbins_aerosol) ) 1294 ALLOCATE( nconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol), & 1295 nconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol), & 1296 nconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 1297 nconc_1 = 0.0_wp 1298 nconc_2 = 0.0_wp 1299 nconc_3 = 0.0_wp 1300 1301 DO i = 1, nbins_aerosol 1302 aerosol_number(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_1(:,:,:,i) 1303 aerosol_number(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_2(:,:,:,i) 1304 aerosol_number(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_3(:,:,:,i) 1305 ALLOCATE( aerosol_number(i)%flux_s(nzb+1:nzt,0:threads_per_task-1), & 1306 aerosol_number(i)%diss_s(nzb+1:nzt,0:threads_per_task-1), & 1307 aerosol_number(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 1308 aerosol_number(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 1309 aerosol_number(i)%init(nzb:nzt+1), & 1310 aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 1311 aerosol_number(i)%init = nclim 1312 IF ( include_emission .OR. ( nldepo .AND. nldepo_surf ) ) THEN 1313 ALLOCATE( aerosol_number(i)%source(nys:nyn,nxl:nxr) ) 1314 aerosol_number(i)%source = 0.0_wp 1315 ENDIF 1316 ENDDO 1317 1318 ! 1319 !-- Aerosol mass concentration 1320 ALLOCATE( aerosol_mass(ncomponents_mass*nbins_aerosol) ) 1321 ALLOCATE( mconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol), & 1322 mconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol), & 1323 mconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol) ) 1324 mconc_1 = 0.0_wp 1325 mconc_2 = 0.0_wp 1326 mconc_3 = 0.0_wp 1327 1328 DO i = 1, ncomponents_mass*nbins_aerosol 1329 aerosol_mass(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_1(:,:,:,i) 1330 aerosol_mass(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_2(:,:,:,i) 1331 aerosol_mass(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_3(:,:,:,i) 1332 ALLOCATE( aerosol_mass(i)%flux_s(nzb+1:nzt,0:threads_per_task-1), & 1333 aerosol_mass(i)%diss_s(nzb+1:nzt,0:threads_per_task-1), & 1334 aerosol_mass(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 1335 aerosol_mass(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 1336 aerosol_mass(i)%init(nzb:nzt+1), & 1337 aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 1338 aerosol_mass(i)%init = mclim 1339 IF ( include_emission .OR. ( nldepo .AND. nldepo_surf ) ) THEN 1340 ALLOCATE( aerosol_mass(i)%source(nys:nyn,nxl:nxr) ) 1341 aerosol_mass(i)%source = 0.0_wp 1342 ENDIF 1343 ENDDO 1344 1345 ! 1346 !-- Surface fluxes: answs = aerosol number, amsws = aerosol mass 1347 ! 1348 !-- Horizontal surfaces: default type 1349 DO l = 0, 2 ! upward (l=0), downward (l=1) and model top (l=2) 1350 ALLOCATE( surf_def_h(l)%answs( 1:surf_def_h(l)%ns, nbins_aerosol ) ) 1351 ALLOCATE( surf_def_h(l)%amsws( 1:surf_def_h(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1352 surf_def_h(l)%answs = 0.0_wp 1353 surf_def_h(l)%amsws = 0.0_wp 1354 ENDDO 1355 ! 1356 !-- Horizontal surfaces: natural type 1357 ALLOCATE( surf_lsm_h%answs( 1:surf_lsm_h%ns, nbins_aerosol ) ) 1358 ALLOCATE( surf_lsm_h%amsws( 1:surf_lsm_h%ns, nbins_aerosol*ncomponents_mass ) ) 1359 surf_lsm_h%answs = 0.0_wp 1360 surf_lsm_h%amsws = 0.0_wp 1361 ! 1362 !-- Horizontal surfaces: urban type 1363 ALLOCATE( surf_usm_h%answs( 1:surf_usm_h%ns, nbins_aerosol ) ) 1364 ALLOCATE( surf_usm_h%amsws( 1:surf_usm_h%ns, nbins_aerosol*ncomponents_mass ) ) 1365 surf_usm_h%answs = 0.0_wp 1366 surf_usm_h%amsws = 0.0_wp 1367 1368 ! 1369 !-- Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and westward (l=3) facing 1370 DO l = 0, 3 1371 ALLOCATE( surf_def_v(l)%answs( 1:surf_def_v(l)%ns, nbins_aerosol ) ) 1372 surf_def_v(l)%answs = 0.0_wp 1373 ALLOCATE( surf_def_v(l)%amsws( 1:surf_def_v(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1374 surf_def_v(l)%amsws = 0.0_wp 1375 1376 ALLOCATE( surf_lsm_v(l)%answs( 1:surf_lsm_v(l)%ns, nbins_aerosol ) ) 1377 surf_lsm_v(l)%answs = 0.0_wp 1378 ALLOCATE( surf_lsm_v(l)%amsws( 1:surf_lsm_v(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1379 surf_lsm_v(l)%amsws = 0.0_wp 1380 1381 ALLOCATE( surf_usm_v(l)%answs( 1:surf_usm_v(l)%ns, nbins_aerosol ) ) 1382 surf_usm_v(l)%answs = 0.0_wp 1383 ALLOCATE( surf_usm_v(l)%amsws( 1:surf_usm_v(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1384 surf_usm_v(l)%amsws = 0.0_wp 1385 1386 ENDDO 1387 1388 ! 1389 !-- Concentration of gaseous tracers (1. SO4, 2. HNO3, 3. NH3, 4. OCNV, 5. OCSV) 1390 !-- (number concentration (#/m3) ) 1391 ! 1392 !-- If chemistry is on, read gas phase concentrations from there. Otherwise, 1393 !-- allocate salsa_gas array. 1394 1395 IF ( air_chemistry ) THEN 1396 DO lsp = 1, nvar 1397 SELECT CASE ( TRIM( chem_species(lsp)%name ) ) 1398 CASE ( 'H2SO4', 'h2so4' ) 1399 gases_available = gases_available + 1 1400 gas_index_chem(1) = lsp 1401 CASE ( 'HNO3', 'hno3' ) 1402 gases_available = gases_available + 1 1403 gas_index_chem(2) = lsp 1404 CASE ( 'NH3', 'nh3' ) 1405 gases_available = gases_available + 1 1406 gas_index_chem(3) = lsp 1407 CASE ( 'OCNV', 'ocnv' ) 1408 gases_available = gases_available + 1 1409 gas_index_chem(4) = lsp 1410 CASE ( 'OCSV', 'ocsv' ) 1411 gases_available = gases_available + 1 1412 gas_index_chem(5) = lsp 1413 END SELECT 1414 ENDDO 1415 1416 IF ( gases_available == ngases_salsa ) THEN 1417 salsa_gases_from_chem = .TRUE. 1418 ELSE 1419 WRITE( message_string, * ) 'SALSA is run together with chemistry but not all gaseous '// & 1420 'components are provided by kpp (H2SO4, HNO3, NH3, OCNV, OCSV)' 1421 CALL message( 'check_parameters', 'PA0599', 1, 2, 0, 6, 0 ) 1422 ENDIF 1423 1424 ELSE 1425 1426 ALLOCATE( salsa_gas(ngases_salsa) ) 1427 ALLOCATE( gconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa), & 1428 gconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa), & 1429 gconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) ) 1430 gconc_1 = 0.0_wp 1431 gconc_2 = 0.0_wp 1432 gconc_3 = 0.0_wp 1433 1434 DO i = 1, ngases_salsa 1435 salsa_gas(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_1(:,:,:,i) 1436 salsa_gas(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_2(:,:,:,i) 1437 salsa_gas(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_3(:,:,:,i) 1438 ALLOCATE( salsa_gas(i)%flux_s(nzb+1:nzt,0:threads_per_task-1), & 1439 salsa_gas(i)%diss_s(nzb+1:nzt,0:threads_per_task-1), & 1440 salsa_gas(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),& 1441 salsa_gas(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),& 1442 salsa_gas(i)%init(nzb:nzt+1), & 1443 salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 1444 salsa_gas(i)%init = nclim 1445 IF ( include_emission ) THEN 1446 ALLOCATE( salsa_gas(i)%source(nys:nys,nxl:nxr) ) 1447 salsa_gas(i)%source = 0.0_wp 1448 ENDIF 1449 ENDDO 1450 ! 1451 !-- Surface fluxes: gtsws = gaseous tracer flux 1452 ! 1453 !-- Horizontal surfaces: default type 1454 DO l = 0, 2 ! upward (l=0), downward (l=1) and model top (l=2) 1455 ALLOCATE( surf_def_h(l)%gtsws( 1:surf_def_h(l)%ns, ngases_salsa ) ) 1456 surf_def_h(l)%gtsws = 0.0_wp 1457 ENDDO 1458 !-- Horizontal surfaces: natural type 1459 ALLOCATE( surf_lsm_h%gtsws( 1:surf_lsm_h%ns, ngases_salsa ) ) 1460 surf_lsm_h%gtsws = 0.0_wp 1461 !-- Horizontal surfaces: urban type 1462 ALLOCATE( surf_usm_h%gtsws( 1:surf_usm_h%ns, ngases_salsa ) ) 1463 surf_usm_h%gtsws = 0.0_wp 1464 ! 1465 !-- Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and 1466 !-- westward (l=3) facing 1467 DO l = 0, 3 1468 ALLOCATE( surf_def_v(l)%gtsws( 1:surf_def_v(l)%ns, ngases_salsa ) ) 1469 surf_def_v(l)%gtsws = 0.0_wp 1470 ALLOCATE( surf_lsm_v(l)%gtsws( 1:surf_lsm_v(l)%ns, ngases_salsa ) ) 1471 surf_lsm_v(l)%gtsws = 0.0_wp 1472 ALLOCATE( surf_usm_v(l)%gtsws( 1:surf_usm_v(l)%ns, ngases_salsa ) ) 1473 surf_usm_v(l)%gtsws = 0.0_wp 1474 ENDDO 1475 ENDIF 1476 1477 IF ( ws_scheme_sca ) THEN 1478 1479 IF ( salsa ) THEN 1480 ALLOCATE( sums_salsa_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 1481 sums_salsa_ws_l = 0.0_wp 1482 ENDIF 1483 1484 ENDIF 1485 ! 1486 !-- Set control flags for decycling only at lateral boundary cores. Within the inner cores the 1487 !-- decycle flag is set to .FALSE.. Even though it does not affect the setting of chemistry boundary 1488 !-- conditions, this flag is used to set advection control flags appropriately. 1489 decycle_salsa_lr = MERGE( decycle_salsa_lr, .FALSE., nxl == 0 .OR. nxr == nx ) 1490 decycle_salsa_ns = MERGE( decycle_salsa_ns, .FALSE., nys == 0 .OR. nyn == ny ) 1491 ! 1492 !-- Decycling can be applied separately for aerosol variables, while wind and other scalars may have 1493 !-- cyclic or nested boundary conditions. However, large gradients near the boundaries may produce 1494 !-- stationary numerical oscillations near the lateral boundaries when a higher-order scheme is 1495 !-- applied near these boundaries. To get rid-off this, set-up additional flags that control the 1496 !-- order of the scalar advection scheme near the lateral boundaries for passive scalars with 1497 !-- decycling. 1498 IF ( scalar_advec == 'ws-scheme' ) THEN 1499 ALLOCATE( salsa_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1500 ! 1501 !-- In case of decycling, set Neuman boundary conditions for wall_flags_0 bit 31 instead of 1502 !-- cyclic boundary conditions. Bit 31 is used to identify extended degradation zones (please see 1503 !-- the following comment). Note, since several also other modules may access this bit but may 1504 !-- have other boundary conditions, the original value of wall_flags_0 bit 31 must not be 1505 !-- modified. Hence, store the boundary conditions directly on salsa_advc_flags_s. 1506 !-- salsa_advc_flags_s will be later overwritten in ws_init_flags_scalar and bit 31 won't be used 1507 !-- to control the numerical order. 1508 !-- Initialize with flag 31 only. 1509 salsa_advc_flags_s = 0 1510 salsa_advc_flags_s = MERGE( IBSET( salsa_advc_flags_s, 31 ), 0, BTEST( wall_flags_0, 31 ) ) 1511 1512 IF ( decycle_salsa_ns ) THEN 1513 IF ( nys == 0 ) THEN 1514 DO i = 1, nbgp 1515 salsa_advc_flags_s(:,nys-i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nys,:), 31 ), & 1516 IBCLR( salsa_advc_flags_s(:,nys,:), 31 ), & 1517 BTEST( salsa_advc_flags_s(:,nys,:), 31 ) ) 1518 ENDDO 1519 ENDIF 1520 IF ( nyn == ny ) THEN 1521 DO i = 1, nbgp 1522 salsa_advc_flags_s(:,nyn+i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nyn,:), 31 ), & 1523 IBCLR( salsa_advc_flags_s(:,nyn,:), 31 ), & 1524 BTEST( salsa_advc_flags_s(:,nyn,:), 31 ) ) 1525 ENDDO 1526 ENDIF 1527 ENDIF 1528 IF ( decycle_salsa_lr ) THEN 1529 IF ( nxl == 0 ) THEN 1530 DO i = 1, nbgp 1531 salsa_advc_flags_s(:,:,nxl-i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxl), 31 ), & 1532 IBCLR( salsa_advc_flags_s(:,:,nxl), 31 ), & 1533 BTEST( salsa_advc_flags_s(:,:,nxl), 31 ) ) 1534 ENDDO 1535 ENDIF 1536 IF ( nxr == nx ) THEN 1537 DO i = 1, nbgp 1538 salsa_advc_flags_s(:,:,nxr+i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxr), 31 ), & 1539 IBCLR( salsa_advc_flags_s(:,:,nxr), 31 ), & 1540 BTEST( salsa_advc_flags_s(:,:,nxr), 31 ) ) 1541 ENDDO 1542 ENDIF 1543 ENDIF 1544 ! 1545 !-- To initialise the advection flags appropriately, pass the boundary flags to 1546 !-- ws_init_flags_scalar. The last argument in ws_init_flags_scalar indicates that a passive 1547 !-- scalar is being treated and the horizontal advection terms are degraded already 2 grid points 1548 !-- before the lateral boundary. Also, extended degradation zones are applied, where 1549 !-- horizontal advection of scalars is discretised by the first-order scheme at all grid points 1550 !-- in the vicinity of buildings (<= 3 grid points). Even though no building is within the 1551 !-- numerical stencil, the first-order scheme is used. At fourth and fifth grid points, the order 1552 !-- of the horizontal advection scheme is successively upgraded. 1553 !-- These degradations of the advection scheme are done to avoid stationary numerical 1554 !-- oscillations, which are responsible for high concentration maxima that may appear e.g. under 1555 !-- shear-free stable conditions. 1556 CALL ws_init_flags_scalar( bc_dirichlet_l .OR. bc_radiation_l .OR. decycle_salsa_lr, & 1557 bc_dirichlet_n .OR. bc_radiation_n .OR. decycle_salsa_ns, & 1558 bc_dirichlet_r .OR. bc_radiation_r .OR. decycle_salsa_lr, & 1559 bc_dirichlet_s .OR. bc_radiation_s .OR. decycle_salsa_ns, & 1560 salsa_advc_flags_s, .TRUE. ) 1561 ENDIF 1562 1563 1564 END SUBROUTINE salsa_init_arrays 1565 1566 !------------------------------------------------------------------------------! 1567 ! Description: 1568 ! ------------ 1569 !> Initialization of SALSA. Based on salsa_initialize in UCLALES-SALSA. 1570 !> Subroutines salsa_initialize, SALSAinit and DiagInitAero in UCLALES-SALSA are 1571 !> also merged here. 1572 !------------------------------------------------------------------------------! 1573 SUBROUTINE salsa_init 1574 1575 IMPLICIT NONE 1576 1577 INTEGER(iwp) :: i !< 1578 INTEGER(iwp) :: ib !< loop index for aerosol number bins 1579 INTEGER(iwp) :: ic !< loop index for aerosol mass bins 1580 INTEGER(iwp) :: ig !< loop index for gases 1581 INTEGER(iwp) :: j !< 1582 1583 IF ( debug_output ) CALL debug_message( 'salsa_init', 'start' ) 1584 1585 bin_low_limits = 0.0_wp 1586 k_topo_top = 0 1587 nsect = 0.0_wp 1588 massacc = 1.0_wp 1542 1589 ! 1543 1590 !-- Initialise 1544 !1545 !-- Aerosol size distribution (TYPE t_section)1546 aero(:)%dwet = 1.0E-10_wp1547 aero(:)%veqh2o = 1.0E-10_wp1548 aero(:)%numc = nclim1549 aero(:)%core = 1.0E-10_wp1550 DO ic = 1, maxspec+1 ! 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O1551 aero(:)%volc(ic) = 0.0_wp1552 ENDDO1553 1554 1591 IF ( nldepo ) sedim_vd = 0.0_wp 1555 1592 … … 1583 1620 !-- Aerosol radius in each bin: dry and wet (m) 1584 1621 ra_dry = 1.0E-10_wp 1585 !1586 !-- Initialise aerosol tracers1587 aero(:)%vhilim = 0.0_wp1588 aero(:)%vlolim = 0.0_wp1589 aero(:)%vratiohi = 0.0_wp1590 aero(:)%vratiolo = 0.0_wp1591 aero(:)%dmid = 0.0_wp1592 !1593 !-- Initialise the sectional particle size distribution1594 CALL set_sizebins1595 1622 ! 1596 1623 !-- Initialise location-dependent aerosol size distributions and chemical compositions: … … 1639 1666 ENDIF 1640 1667 ENDIF 1668 ! 1669 !-- Partition and dissolutional growth by gaseous HNO3 and NH3 1670 IF ( index_no > 0 .AND. index_nh > 0 .AND. index_so4 > 0 ) lspartition = .TRUE. 1641 1671 1642 1672 IF ( debug_output ) CALL debug_message( 'salsa_init', 'end' ) … … 1676 1706 1677 1707 REAL(wp) :: ratio_d !< ratio of the upper and lower diameter of subranges 1708 1709 aero(:)%dwet = 1.0E-10_wp 1710 aero(:)%veqh2o = 1.0E-10_wp 1711 aero(:)%numc = nclim 1712 aero(:)%core = 1.0E-10_wp 1713 DO cc = 1, maxspec+1 ! 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O 1714 aero(:)%volc(cc) = 0.0_wp 1715 ENDDO 1678 1716 ! 1679 1717 !-- vlolim&vhilim: min & max *dry* volumes [fxm] … … 1737 1775 1738 1776 USE netcdf_data_input_mod, & 1739 ONLY: close_input_file, get_attribute, get_variable, & 1777 ONLY: check_existence, close_input_file, get_attribute, get_variable, & 1778 inquire_num_variables, inquire_variable_names, & 1740 1779 netcdf_data_input_get_dimension_length, open_read_file 1741 1780 1742 1781 IMPLICIT NONE 1743 1782 1744 CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE :: cc_name !< chemical component name 1783 CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE :: cc_name !< chemical component name 1784 CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE :: var_names !< variable names 1745 1785 1746 1786 INTEGER(iwp) :: ee !< index: end … … 1753 1793 INTEGER(iwp) :: k !< loop index: z-direction 1754 1794 INTEGER(iwp) :: lod_aero !< level of detail of inital aerosol concentrations 1755 INTEGER(iwp) :: pr_nbins !< Number of aerosol size bins in file 1756 INTEGER(iwp) :: pr_ncc !< Number of aerosol chemical components in file 1757 INTEGER(iwp) :: pr_nz !< Number of vertical grid-points in file 1795 INTEGER(iwp) :: num_vars !< number of variables 1796 INTEGER(iwp) :: pr_nbins !< number of aerosol size bins in file 1797 INTEGER(iwp) :: pr_ncc !< number of aerosol chemical components in file 1798 INTEGER(iwp) :: pr_nz !< number of vertical grid-points in file 1758 1799 INTEGER(iwp) :: prunmode !< running mode of SALSA 1759 1800 INTEGER(iwp) :: ss !< index: start 1760 1801 1761 INTEGER(iwp), DIMENSION(maxspec) :: cc_in put_to_model1802 INTEGER(iwp), DIMENSION(maxspec) :: cc_in2mod 1762 1803 1763 1804 LOGICAL :: netcdf_extend = .FALSE. !< Flag: netcdf file exists … … 1770 1811 REAL(wp), DIMENSION(0:nz+1) :: pmfoc1a !< mass fraction of OC in 1a 1771 1812 1772 REAL(wp), DIMENSION(0:nz+1,nbins_aerosol) :: pndist !< size dist as a function of height(#/m3)1813 REAL(wp), DIMENSION(0:nz+1,nbins_aerosol) :: pndist !< vertical profile of size dist. (#/m3) 1773 1814 REAL(wp), DIMENSION(0:nz+1,maxspec) :: pmf2a !< mass distributions in subrange 2a 1774 1815 REAL(wp), DIMENSION(0:nz+1,maxspec) :: pmf2b !< mass distributions in subrange 2b … … 1780 1821 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pr_mass_fracs_b !< and b 1781 1822 1782 cc_in put_to_model= 01823 cc_in2mod = 0 1783 1824 prunmode = 1 1784 1825 ! 1785 1826 !-- Bin mean aerosol particle volume (m3) 1786 core(:) = 0.0_wp1787 1827 core(1:nbins_aerosol) = api6 * aero(1:nbins_aerosol)%dmid**3 1788 1828 ! … … 1806 1846 CALL open_read_file( TRIM( input_file_dynamic ) // TRIM( coupling_char ), id_dyn ) 1807 1847 ! 1808 !-- Inquire dimensions: 1848 !-- At first, inquire all variable names 1849 CALL inquire_num_variables( id_dyn, num_vars ) 1850 ! 1851 !-- Allocate memory to store variable names 1852 ALLOCATE( var_names(1:num_vars) ) 1853 CALL inquire_variable_names( id_dyn, var_names ) 1854 ! 1855 !-- Inquire vertical dimension and number of aerosol chemical components 1809 1856 CALL netcdf_data_input_get_dimension_length( id_dyn, pr_nz, 'z' ) 1810 1857 IF ( pr_nz /= nz ) THEN … … 1824 1871 CALL get_variable( id_dyn, 'z', pr_z ) 1825 1872 ! 1826 !-- Read name and index of chemical components 1827 CALL get_variable( id_dyn, 'composition_name', cc_name, pr_ncc ) 1873 !-- Read the names of chemical components 1874 IF ( check_existence( var_names, 'composition_name' ) ) THEN 1875 CALL get_variable( id_dyn, 'composition_name', cc_name, pr_ncc ) 1876 ELSE 1877 WRITE( message_string, * ) 'Missing composition_name in ' // TRIM( input_file_dynamic ) 1878 CALL message( 'aerosol_init', 'PA0655', 1, 2, 0, 6, 0 ) 1879 ENDIF 1880 ! 1881 !-- Define the index of each chemical component in the model 1828 1882 DO ic = 1, pr_ncc 1829 1883 SELECT CASE ( TRIM( cc_name(ic) ) ) 1830 1884 CASE ( 'H2SO4', 'SO4', 'h2so4', 'so4' ) 1831 cc_in put_to_model(1) = ic1885 cc_in2mod(1) = ic 1832 1886 CASE ( 'OC', 'oc' ) 1833 cc_in put_to_model(2) = ic1887 cc_in2mod(2) = ic 1834 1888 CASE ( 'BC', 'bc' ) 1835 cc_in put_to_model(3) = ic1889 cc_in2mod(3) = ic 1836 1890 CASE ( 'DU', 'du' ) 1837 cc_in put_to_model(4) = ic1891 cc_in2mod(4) = ic 1838 1892 CASE ( 'SS', 'ss' ) 1839 cc_in put_to_model(5) = ic1893 cc_in2mod(5) = ic 1840 1894 CASE ( 'HNO3', 'hno3', 'NO', 'no' ) 1841 cc_in put_to_model(6) = ic1895 cc_in2mod(6) = ic 1842 1896 CASE ( 'NH3', 'nh3', 'NH', 'nh' ) 1843 cc_in put_to_model(7) = ic1897 cc_in2mod(7) = ic 1844 1898 END SELECT 1845 1899 ENDDO 1846 1900 1847 IF ( SUM( cc_in put_to_model) == 0 ) THEN1901 IF ( SUM( cc_in2mod ) == 0 ) THEN 1848 1902 message_string = 'None of the aerosol chemical components in ' // TRIM( & 1849 1903 input_file_dynamic ) // ' correspond to ones applied in SALSA.' … … 1852 1906 ! 1853 1907 !-- Vertical profiles of mass fractions of different chemical components: 1854 CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_a', pr_mass_fracs_a, & 1855 0, pr_ncc-1, 0, pr_nz-1 ) 1908 IF ( check_existence( var_names, 'init_atmosphere_mass_fracs_a' ) ) THEN 1909 CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_a', pr_mass_fracs_a, & 1910 0, pr_ncc-1, 0, pr_nz-1 ) 1911 ELSE 1912 WRITE( message_string, * ) 'Missing init_atmosphere_mass_fracs_a in ' // & 1913 TRIM( input_file_dynamic ) 1914 CALL message( 'aerosol_init', 'PA0656', 1, 2, 0, 6, 0 ) 1915 ENDIF 1856 1916 CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_b', pr_mass_fracs_b, & 1857 1917 0, pr_ncc-1, 0, pr_nz-1 ) … … 1859 1919 !-- Match the input data with the chemical composition applied in the model 1860 1920 DO ic = 1, maxspec 1861 ss = cc_in put_to_model(ic)1921 ss = cc_in2mod(ic) 1862 1922 IF ( ss == 0 ) CYCLE 1863 1923 pmf2a(nzb+1:nzt+1,ic) = pr_mass_fracs_a(nzb:nzt,ss) … … 1865 1925 ENDDO 1866 1926 ! 1867 !-- Aerosol concentrations: lod=1 ( total PM) or lod=2 (sectional number size distribution)1927 !-- Aerosol concentrations: lod=1 (vertical profile of sectional number size distribution) 1868 1928 CALL get_attribute( id_dyn, 'lod', lod_aero, .FALSE., 'init_atmosphere_aerosol' ) 1869 1929 IF ( lod_aero /= 1 ) THEN … … 3020 3080 !-- Condensation 3021 3081 IF ( lscnd ) THEN 3022 CALL condensation( lo_aero, zgso4, zgocnv, zgocsv, zghno3, zgnh3, in_cw(k), in_cs(k), 3082 CALL condensation( lo_aero, zgso4, zgocnv, zgocsv, zghno3, zgnh3, in_cw(k), in_cs(k), & 3023 3083 in_t(k), in_p(k), dt_salsa, prtcl ) 3024 3084 ENDIF … … 3026 3086 !-- Deposition 3027 3087 IF ( lsdepo ) THEN 3028 CALL deposition( lo_aero, in_t(k), in_adn(k), in_u(k), in_lad, kvis(k), schmidt_num(k,:), 3088 CALL deposition( lo_aero, in_t(k), in_adn(k), in_u(k), in_lad, kvis(k), schmidt_num(k,:),& 3029 3089 vd(k,:) ) 3030 3090 ENDIF … … 3258 3318 ! Description: 3259 3319 ! ------------ 3260 !> Set logical switches according to the host model state and user-specified 3261 !> NAMELIST options. 3320 !> Set logical switches according to the salsa_parameters options. 3262 3321 !> Juha Tonttila, FMI, 2014 3263 3322 !> Only aerosol processes included, Mona Kurppa, UHel, 2017 … … 3832 3891 ! 3833 3892 !-- Brownian diffusion 3834 v_bd = mag_u * c_brownian_diff * schmidt_num**( -0.66666666_wp ) * 3893 v_bd = mag_u * c_brownian_diff * schmidt_num**( -0.66666666_wp ) * & 3835 3894 ( mag_u * par_l / kvis_a )**( -0.5_wp ) 3836 3895 ! 3837 3896 !-- Interception 3838 v_in = mag_u * c_interception * diameter / par_l * ( 2.0_wp + LOG( 2.0_wp * par_l / diameter ) ) 3897 v_in = mag_u * c_interception * diameter / par_l * & 3898 ( 2.0_wp + LOG( 2.0_wp * par_l / diameter ) ) 3839 3899 ! 3840 3900 !-- Impaction: Petroff (2009) Eq. 18 … … 4179 4239 REAL(wp), DIMENSION(maxspec+1) :: zplusterm !< coagulation gain in a bin (for each 4180 4240 !< chemical compound) 4181 REAL(wp), DIMENSION(nbins_aerosol,nbins_aerosol) :: zcc !< updated coagulation coeff icients(m3/s)4241 REAL(wp), DIMENSION(nbins_aerosol,nbins_aerosol) :: zcc !< updated coagulation coeff. (m3/s) 4182 4242 4183 4243 TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) :: paero !< Aerosol properties … … 4186 4246 zdpart_nn = 0.0_wp 4187 4247 ! 4188 !-- 1) Coagulation to coarse mode calculated in a simplified way: 4189 !-- CoagSink ~ Dp in continuum subrange, thus we calculate 'effective' 4190 !-- number concentration of coarse particles 4248 !-- 1) Coagulation to coarse mode calculated in a simplified way: 4249 !-- CoagSink ~ Dp in continuum subrange --> 'effective' number conc. of coarse particles 4191 4250 4192 4251 !-- 2) Updating coagulation coefficients … … 4418 4477 fmdist = SQRT( tva(1)**2 + tva(2)**2 ) 4419 4478 ! 4420 !-- 5) Coagulation coefficient = coalescence efficiency * collision kernel (m3/s) (eq. 15.33). 4479 !-- 5) Coagulation coefficient = coalescence efficiency * collision kernel (m3/s) (eq. 15.33). 4421 4480 !-- Here assumed coalescence efficiency 1!! 4422 4481 coagc = flux(1) / ( mdiam / ( mdiam + fmdist) + flux(1) / flux(2) ) … … 4632 4691 ! 4633 4692 !-- 5.1.2) Semivolatile organic compound: all bins except subrange 1 4634 zcs_ocsv = SUM( zcolrate(start_subrange_2a:end_subrange_2b) ) !< sink for semi-volatile org anics4693 zcs_ocsv = SUM( zcolrate(start_subrange_2a:end_subrange_2b) ) !< sink for semi-volatile org. 4635 4694 IF ( pcocsv > 1.0E+10_wp .AND. zcs_ocsv > 1.0E-30 .AND. is_used( prtcl,'OC') ) THEN 4636 4695 ! … … 4849 4908 zkocnv = 0.0_wp 4850 4909 4910 zc_h2so4 = pc_sa * 1.0E-6_wp ! sulphuric acid conc. to #/cm3 4911 zc_org = pc_ocnv * 1.0E-6_wp ! conc. of non-volatile OC to #/cm3 4912 zmixnh3 = pc_nh3 * ptemp * argas / ( ppres * avo ) 4913 4851 4914 SELECT CASE ( nsnucl ) 4852 4915 ! … … 4854 4917 CASE(1) 4855 4918 4856 zc_h2so4 = pc_sa * 1.0E-6_wp ! sulphuric acid conc. to #/cm34857 4919 CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv ) 4858 4920 ! 4859 !-- Activation type nucleation 4921 !-- Activation type nucleation (See Riipinen et al. (2007), Atmos. Chem. Phys., 7(8), 1899-1914) 4860 4922 CASE(2) 4861 4862 zc_h2so4 = pc_sa * 1.0E-6_wp ! sulphuric acid conc. to #/cm3 4863 CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv ) 4864 CALL actnucl( pc_sa, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv, act_coeff ) 4923 ! 4924 !-- Nucleation rate (#/(m3 s)) 4925 zc_h2so4 = MAX( zc_h2so4, 1.0E4_wp ) 4926 zc_h2so4 = MIN( zc_h2so4, 1.0E11_wp ) 4927 zjnuc = act_coeff * pc_sa ! (#/(m3 s)) 4928 ! 4929 !-- Organic compounds not involved when kinetic nucleation is assumed. 4930 zdcrit = 7.9375E-10_wp ! (m) 4931 zkocnv = 0.0_wp 4932 zksa = 1.0_wp 4933 znoc = 0.0_wp 4934 znsa = 2.0_wp 4865 4935 ! 4866 4936 !-- Kinetically limited nucleation of (NH4)HSO4 clusters 4937 !-- (See Sihto et al. (2006), Atmos. Chem. Phys., 6(12), 4079-4091.) 4867 4938 CASE(3) 4868 4869 zc_h2so4 = pc_sa * 1.0E-6_wp ! sulphuric acid conc. to #/cm3 4870 CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv ) 4871 CALL kinnucl( zc_h2so4, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv ) 4939 ! 4940 !-- Nucleation rate = coagcoeff*zpcsa**2 (#/(m3 s)) 4941 zc_h2so4 = MAX( zc_h2so4, 1.0E4_wp ) 4942 zc_h2so4 = MIN( zc_h2so4, 1.0E11_wp ) 4943 zjnuc = 5.0E-13_wp * zc_h2so4**2.0_wp * 1.0E+6_wp 4944 ! 4945 !-- Organic compounds not involved when kinetic nucleation is assumed. 4946 zdcrit = 7.9375E-10_wp ! (m) 4947 zkocnv = 0.0_wp 4948 zksa = 1.0_wp 4949 znoc = 0.0_wp 4950 znsa = 2.0_wp 4872 4951 ! 4873 4952 !-- Ternary H2SO4-H2O-NH3 nucleation 4874 4953 CASE(4) 4875 4954 4876 zmixnh3 = pc_nh3 * ptemp * argas / ( ppres * avo )4877 zc_h2so4 = pc_sa * 1.0E-6_wp ! sulphuric acid conc. to #/cm34878 4955 CALL ternucl( zc_h2so4, zmixnh3, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv ) 4879 4956 ! 4880 4957 !-- Organic nucleation, J~[ORG] or J~[ORG]**2 4958 !-- (See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.) 4881 4959 CASE(5) 4882 4883 zc_org = pc_ocnv * 1.0E-6_wp ! conc. of non-volatile OC to #/cm3 4884 zc_h2so4 = pc_sa * 1.0E-6_wp ! sulphuric acid conc. to #/cm3 4885 CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv ) 4886 CALL orgnucl( pc_ocnv, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv ) 4960 ! 4961 !-- Homomolecular nuleation rate 4962 zjnuc = 1.3E-7_wp * pc_ocnv ! (1/s) (Paasonen et al. Table 4: median a_org) 4963 ! 4964 !-- H2SO4 not involved when pure organic nucleation is assumed. 4965 zdcrit = 1.5E-9 ! (m) 4966 zkocnv = 1.0_wp 4967 zksa = 0.0_wp 4968 znoc = 1.0_wp 4969 znsa = 0.0_wp 4887 4970 ! 4888 4971 !-- Sum of H2SO4 and organic activation type nucleation, J~[H2SO4]+[ORG] 4972 !-- (See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242) 4889 4973 CASE(6) 4890 4891 zc_h2so4 = pc_sa * 1.0E-6_wp ! sulphuric acid conc. to #/cm3 4892 CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv ) 4893 CALL sumnucl( pc_sa, pc_ocnv, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv ) 4894 ! 4895 !-- Heteromolecular nucleation, J~[H2SO4]*[ORG] 4974 ! 4975 !-- Nucleation rate (#/m3/s) 4976 zjnuc = 6.1E-7_wp * pc_sa + 0.39E-7_wp * pc_ocnv ! (Paasonen et al. Table 3.) 4977 ! 4978 !-- Both organic compounds and H2SO4 are involved when sumnucleation is assumed. 4979 zdcrit = 1.5E-9_wp ! (m) 4980 zkocnv = 1.0_wp 4981 zksa = 1.0_wp 4982 znoc = 1.0_wp 4983 znsa = 1.0_wp 4984 ! 4985 !-- Heteromolecular nucleation, J~[H2SO4]*[ORG] 4986 !-- (See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.) 4896 4987 CASE(7) 4897 4898 zc_h2so4 = pc_sa * 1.0E-6_wp ! sulphuric acid conc. to #/cm3 4899 zc_org = pc_ocnv * 1.0E-6_wp ! conc. of non-volatile OC to #/cm3 4900 CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv ) 4901 CALL hetnucl( zc_h2so4, zc_org, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv ) 4988 ! 4989 !-- Nucleation rate (#/m3/s) 4990 zjnuc = 4.1E-14_wp * pc_sa * pc_ocnv * 1.0E6_wp ! (Paasonen et al. Table 4: median) 4991 ! 4992 !-- Both organic compounds and H2SO4 are involved when heteromolecular nucleation is assumed 4993 zdcrit = 1.5E-9_wp ! (m) 4994 zkocnv = 1.0_wp 4995 zksa = 1.0_wp 4996 znoc = 1.0_wp 4997 znsa = 1.0_wp 4902 4998 ! 4903 4999 !-- Homomolecular nucleation of H2SO4 and heteromolecular nucleation of H2SO4 and organic vapour, 4904 5000 !-- J~[H2SO4]**2 + [H2SO4]*[ORG] (EUCAARI project) 5001 !-- (See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242) 4905 5002 CASE(8) 4906 zc_h2so4 = pc_sa * 1.0E-6_wp ! sulphuric acid conc. to #/cm3 4907 zc_org = pc_ocnv * 1.0E-6_wp ! conc. of non-volatile OC to #/cm3 4908 CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv ) 4909 CALL SAnucl( zc_h2so4, zc_org, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv ) 5003 ! 5004 !-- Nucleation rate (#/m3/s) 5005 zjnuc = ( 1.1E-14_wp * zc_h2so4**2 + 3.2E-14_wp * zc_h2so4 * zc_org ) * 1.0E+6_wp 5006 ! 5007 !-- Both organic compounds and H2SO4 are involved when SAnucleation is assumed 5008 zdcrit = 1.5E-9_wp ! (m) 5009 zkocnv = 1.0_wp 5010 zksa = 1.0_wp 5011 znoc = 1.0_wp 5012 znsa = 3.0_wp 4910 5013 ! 4911 5014 !-- Homomolecular nucleation of H2SO4 and organic vapour and heteromolecular nucleation of H2SO4 4912 5015 !-- and organic vapour, J~[H2SO4]**2 + [H2SO4]*[ORG]+[ORG]**2 (EUCAARI project) 4913 5016 CASE(9) 4914 4915 zc_h2so4 = pc_sa * 1.0E-6_wp ! sulphuric acid conc. to #/cm3 4916 zc_org = pc_ocnv * 1.0E-6_wp ! conc. of non-volatile OC to #/cm3 4917 CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv ) 4918 CALL SAORGnucl( zc_h2so4, zc_org, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv ) 5017 ! 5018 !-- Nucleation rate (#/m3/s) 5019 zjnuc = ( 1.4E-14_wp * zc_h2so4**2 + 2.6E-14_wp * zc_h2so4 * zc_org + 0.037E-14_wp * & 5020 zc_org**2 ) * 1.0E+6_wp 5021 ! 5022 !-- Both organic compounds and H2SO4 are involved when SAORGnucleation is assumed 5023 zdcrit = 1.5E-9_wp ! (m) 5024 zkocnv = 1.0_wp 5025 zksa = 1.0_wp 5026 znoc = 3.0_wp 5027 znsa = 3.0_wp 4919 5028 4920 5029 END SELECT … … 5010 5119 ! 5011 5120 !-- Defining the parameter m (zm_para) for calculating the coagulation sink onto background 5012 !-- particles (Eq. 5&6 in Lehtinen et al. 2007). The growth is investigated between 5121 !-- particles (Eq. 5&6 in Lehtinen et al. 2007). The growth is investigated between 5013 5122 !-- [d1,reglim(1)] = [zdcrit,3nm] and m = LOG( CoagS_dx / CoagX_zdcrit ) / LOG( reglim / zdcrit ) 5014 5123 !-- (Lehtinen et al. 2007, Eq. 6). 5015 !-- The steps for the coagulation sink for reglim = 3nm and zdcrit ~= 1nm are explained in 5016 !-- Kulmala et al. (2001). The particles of diameter zdcrit ~1.14 nm and reglim = 3nm are both 5124 !-- The steps for the coagulation sink for reglim = 3nm and zdcrit ~= 1nm are explained in 5125 !-- Kulmala et al. (2001). The particles of diameter zdcrit ~1.14 nm and reglim = 3nm are both 5017 5126 !-- in turn the "number 1" variables (Kulmala et al. 2001). 5018 5127 !-- c = critical (1nm), x = 3nm, 2 = wet or mean droplet … … 5129 5238 z_gr_tot = z_gr_clust * 2.77777777E-7_wp + 1.5708E-6_wp * zlambda * zdcrit**3 * & 5130 5239 ( z_n_nuc * 1.0E-6_wp ) * zcv_c * avo * 2.77777777E-7_wp ! (Eq. 3) 5131 zeta = - zcoagstot / ( ( zm_para + 1.0_wp ) * z_gr_tot * ( zdcrit**zm_para ) ) ! (Eq. 5240 zeta = - zcoagstot / ( ( zm_para + 1.0_wp ) * z_gr_tot * ( zdcrit**zm_para ) ) ! (Eq.7b) 5132 5241 ! 5133 5242 !-- Calculate Eq. 7a (Taylor series for the number of particles between [d1,dx]) … … 5139 5248 z_gr_tot = z_gr_clust * 1.0E-9_wp / 3600.0_wp + 1.5708E-6_wp * zlambda * zdcrit**3 * & 5140 5249 ( z_n_nuc * 1.0E-6_wp ) * zcv_c * avo * 1.0E-9_wp / 3600.0_wp !< (m/s) 5141 zj3 = zjnuc * EXP( MIN( 0.0_wp, -zgamma * zdcrit * zcoagstot / z_gr_tot ) ) ! (#/m3s, Eq. 5250 zj3 = zjnuc * EXP( MIN( 0.0_wp, -zgamma * zdcrit * zcoagstot / z_gr_tot ) ) ! (#/m3s, Eq.5a) 5142 5251 5143 5252 ENDIF … … 5145 5254 ! 5146 5255 !-- If J3 very small (< 1 #/cm3), neglect particle formation. In real atmosphere this would mean 5147 !-- that clusters form but coagulate to pre-existing particles who gain sulphate. Since 5148 !-- CoagS ~ CS (4piD*CS'), we do *not* update H2SO4 concentration here but let condensation take 5256 !-- that clusters form but coagulate to pre-existing particles who gain sulphate. Since 5257 !-- CoagS ~ CS (4piD*CS'), we do *not* update H2SO4 concentration here but let condensation take 5149 5258 !-- care of it. Formation mass rate of molecules (molec/m3s) for 1: H2SO4 and 2: organic vapour 5150 5259 pj3n3(1) = zj3 * n3 * pxsa … … 5156 5265 ! Description: 5157 5266 ! ------------ 5158 !> Calculate the nucleation rate and the size of critical clusters assuming 5267 !> Calculate the nucleation rate and the size of critical clusters assuming 5159 5268 !> binary nucleation. 5160 !> Parametrisation according to Vehkamaki et al. (2002), J. Geophys. Res., 5161 !> 107(D22), 4622. Called from subroutine nucleation. 5269 !> Parametrisation according to Vehkamaki et al. (2002), J. Geophys. Res., 5270 !> 107(D22), 4622. Called from subroutine nucleation. 5162 5271 !------------------------------------------------------------------------------! 5163 5272 SUBROUTINE binnucl( pc_sa, ptemp, prh, pnuc_rate, pn_crit_sa, pn_crit_ocnv, pd_crit, pk_sa, & … … 5519 5628 5520 5629 END SUBROUTINE ternucl 5521 5522 !------------------------------------------------------------------------------!5523 ! Description:5524 ! ------------5525 !> Calculate the nucleation rate and the size of critical clusters assuming5526 !> kinetic nucleation. Each sulphuric acid molecule forms an (NH4)HSO4 molecule5527 !> in the atmosphere and two colliding (NH4)HSO4 molecules form a stable5528 !> cluster. See Sihto et al. (2006), Atmos. Chem. Phys., 6(12), 4079-4091.5529 !>5530 !> Below the following assumption have been made:5531 !> nucrate = coagcoeff*zpcsa**25532 !> coagcoeff = 8*sqrt(3*boltz*ptemp*r_abs/dens_abs)5533 !> r_abs = 0.315d-9 radius of bisulphate molecule [m]5534 !> dens_abs = 1465 density of - " - [kg/m3]5535 !------------------------------------------------------------------------------!5536 SUBROUTINE kinnucl( pc_sa, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv )5537 5538 IMPLICIT NONE5539 5540 REAL(wp), INTENT(in) :: pc_sa !< H2SO4 conc. (#/m3)5541 5542 REAL(wp), INTENT(out) :: pd_crit !< critical diameter of clusters (m)5543 REAL(wp), INTENT(out) :: pk_ocnv !< if pk_ocnv = 1, organic compounds participate in nucleation5544 REAL(wp), INTENT(out) :: pk_sa !< if pk_sa = 1, H2SO4 is participate in nucleation5545 REAL(wp), INTENT(out) :: pn_crit_ocnv !< number of organic molecules in cluster (#)5546 REAL(wp), INTENT(out) :: pn_crit_sa !< number of H2SO4 molecules in cluster (#)5547 REAL(wp), INTENT(out) :: pnuc_rate !< nucl. rate (#/(m3 s))5548 !5549 !-- Nucleation rate (#/(m3 s))5550 pnuc_rate = 5.0E-13_wp * pc_sa**2.0_wp * 1.0E+6_wp5551 !5552 !-- Organic compounds not involved when kinetic nucleation is assumed.5553 pn_crit_sa = 2.0_wp5554 pn_crit_ocnv = 0.0_wp5555 pk_sa = 1.0_wp5556 pk_ocnv = 0.0_wp5557 pd_crit = 7.9375E-10_wp ! (m)5558 5559 END SUBROUTINE kinnucl5560 5561 !------------------------------------------------------------------------------!5562 ! Description:5563 ! ------------5564 !> Calculate the nucleation rate and the size of critical clusters assuming5565 !> activation type nucleation.5566 !> See Riipinen et al. (2007), Atmos. Chem. Phys., 7(8), 1899-1914.5567 !------------------------------------------------------------------------------!5568 SUBROUTINE actnucl( psa_conc, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv, activ )5569 5570 IMPLICIT NONE5571 5572 REAL(wp), INTENT(in) :: activ !< activation coefficient (1e-7 by default)5573 REAL(wp), INTENT(in) :: psa_conc !< H2SO4 conc. (#/m3)5574 5575 REAL(wp), INTENT(out) :: pd_crit !< critical diameter of clusters (m)5576 REAL(wp), INTENT(out) :: pk_ocnv !< if pk_ocnv = 1, organic compounds participate in nucleation5577 REAL(wp), INTENT(out) :: pk_sa !< if pk_sa = 1, H2SO4 participate in nucleation5578 REAL(wp), INTENT(out) :: pn_crit_ocnv !< number of organic molecules in cluster (#)5579 REAL(wp), INTENT(out) :: pn_crit_sa !< number of H2SO4 molecules in cluster (#)5580 REAL(wp), INTENT(out) :: pnuc_rate !< nucl. rate (#/(m3 s))5581 !5582 !-- Nucleation rate (#/(m3 s))5583 pnuc_rate = activ * psa_conc ! (#/(m3 s))5584 !5585 !-- Organic compounds not involved when kinetic nucleation is assumed.5586 pn_crit_sa = 2.0_wp5587 pn_crit_ocnv = 0.0_wp5588 pk_sa = 1.0_wp5589 pk_ocnv = 0.0_wp5590 pd_crit = 7.9375E-10_wp ! (m)5591 5592 END SUBROUTINE actnucl5593 5594 !------------------------------------------------------------------------------!5595 ! Description:5596 ! ------------5597 !> Conciders only the organic matter in nucleation. Paasonen et al. (2010)5598 !> determined particle formation rates for 2 nm particles, J2, from different5599 !> kind of combinations of sulphuric acid and organic matter concentration.5600 !> See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.5601 !------------------------------------------------------------------------------!5602 SUBROUTINE orgnucl( pc_org, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, &5603 pk_sa, pk_ocnv )5604 5605 IMPLICIT NONE5606 5607 REAL(wp) :: a_org = 1.3E-7_wp !< (1/s) (Paasonen et al. Table 4: median)5608 5609 REAL(wp), INTENT(in) :: pc_org !< organic vapour concentration (#/m3)5610 5611 REAL(wp), INTENT(out) :: pd_crit !< critical diameter of clusters (m)5612 REAL(wp), INTENT(out) :: pk_ocnv !< if pk_ocnv = 1, organic compounds participate in nucleation5613 REAL(wp), INTENT(out) :: pk_sa !< if pk_sa = 1, H2SO4 participate in nucleation5614 REAL(wp), INTENT(out) :: pn_crit_ocnv !< number of organic molecules in cluster (#)5615 REAL(wp), INTENT(out) :: pn_crit_sa !< number of H2SO4 molecules in cluster (#)5616 REAL(wp), INTENT(out) :: pnuc_rate !< nucl. rate (#/(m3 s))5617 !5618 !-- Homomolecular nuleation rate5619 pnuc_rate = a_org * pc_org5620 !5621 !-- H2SO4 not involved when pure organic nucleation is assumed.5622 pn_crit_sa = 0.0_wp5623 pn_crit_ocnv = 1.0_wp5624 pk_sa = 0.0_wp5625 pk_ocnv = 1.0_wp5626 pd_crit = 1.5E-9_wp ! (m)5627 5628 END SUBROUTINE orgnucl5629 5630 !------------------------------------------------------------------------------!5631 ! Description:5632 ! ------------5633 !> Conciders both the organic vapor and H2SO4 in nucleation - activation type5634 !> of nucleation.5635 !> See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.5636 !------------------------------------------------------------------------------!5637 SUBROUTINE sumnucl( pc_sa, pc_org, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv )5638 5639 IMPLICIT NONE5640 5641 REAL(wp) :: a_s1 = 6.1E-7_wp !< (1/s)5642 REAL(wp) :: a_s2 = 0.39E-7_wp !< (1/s) (Paasonen et al. Table 3.)5643 5644 REAL(wp), INTENT(in) :: pc_org !< organic vapour concentration (#/m3)5645 REAL(wp), INTENT(in) :: pc_sa !< H2SO4 conc. (#/m3)5646 5647 REAL(wp), INTENT(out) :: pd_crit !< critical diameter of clusters (m)5648 REAL(wp), INTENT(out) :: pk_ocnv !< if pk_ocnv = 1, organic compounds participate in nucleation5649 REAL(wp), INTENT(out) :: pk_sa !< if pk_sa = 1, H2SO4 participate in nucleation5650 REAL(wp), INTENT(out) :: pn_crit_ocnv !< number of organic molecules in cluster (#)5651 REAL(wp), INTENT(out) :: pn_crit_sa !< number of H2SO4 molecules in cluster (#)5652 REAL(wp), INTENT(out) :: pnuc_rate !< nucl. rate (#/(m3 s))5653 !5654 !-- Nucleation rate (#/m3/s)5655 pnuc_rate = a_s1 * pc_sa + a_s2 * pc_org5656 !5657 !-- Both organic compounds and H2SO4 are involved when sumnucleation is assumed.5658 pn_crit_sa = 1.0_wp5659 pn_crit_ocnv = 1.0_wp5660 pk_sa = 1.0_wp5661 pk_ocnv = 1.0_wp5662 pd_crit = 1.5E-9_wp ! (m)5663 5664 END SUBROUTINE sumnucl5665 5666 !------------------------------------------------------------------------------!5667 ! Description:5668 ! ------------5669 !> Conciders both the organic vapor and H2SO4 in nucleation - heteromolecular5670 !> nucleation.5671 !> See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.5672 !------------------------------------------------------------------------------!5673 SUBROUTINE hetnucl( pc_sa, pc_org, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv )5674 5675 IMPLICIT NONE5676 5677 REAL(wp) :: z_k_het = 4.1E-14_wp !< (cm3/s) (Paasonen et al. Table 4: median)5678 5679 REAL(wp), INTENT(in) :: pc_org !< organic vapour concentration (#/m3)5680 REAL(wp), INTENT(in) :: pc_sa !< H2SO4 conc. (#/m3)5681 5682 REAL(wp), INTENT(out) :: pd_crit !< critical diameter of clusters (m)5683 REAL(wp), INTENT(out) :: pk_ocnv !< if pk_ocnv = 1, organic compounds participate in nucleation5684 REAL(wp), INTENT(out) :: pk_sa !< if pk_sa = 1, H2SO4 participate in nucleation5685 REAL(wp), INTENT(out) :: pn_crit_ocnv !< number of organic molecules in cluster (#)5686 REAL(wp), INTENT(out) :: pn_crit_sa !< number of H2SO4 molecules in cluster (#)5687 REAL(wp), INTENT(out) :: pnuc_rate !< nucl. rate (#/(m3 s))5688 !5689 !-- Nucleation rate (#/m3/s)5690 pnuc_rate = z_k_het * pc_sa * pc_org * 1.0E6_wp5691 !5692 !-- Both organic compounds and H2SO4 are involved when heteromolecular5693 !-- nucleation is assumed.5694 pn_crit_sa = 1.0_wp5695 pn_crit_ocnv = 1.0_wp5696 pk_sa = 1.0_wp5697 pk_ocnv = 1.0_wp5698 pd_crit = 1.5E-9_wp ! (m)5699 5700 END SUBROUTINE hetnucl5701 5702 !------------------------------------------------------------------------------!5703 ! Description:5704 ! ------------5705 !> Takes into account the homomolecular nucleation of sulphuric acid H2SO4 with5706 !> both of the available vapours.5707 !> See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.5708 !------------------------------------------------------------------------------!5709 SUBROUTINE SAnucl( pc_sa, pc_org, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv )5710 5711 IMPLICIT NONE5712 5713 REAL(wp) :: z_k_sa1 = 1.1E-14_wp !< (cm3/s)5714 REAL(wp) :: z_k_sa2 = 3.2E-14_wp !< (cm3/s) (Paasonen et al. Table 3.)5715 5716 REAL(wp), INTENT(in) :: pc_org !< organic vapour concentration (#/m3)5717 REAL(wp), INTENT(in) :: pc_sa !< H2SO4 conc. (#/m3)5718 5719 REAL(wp), INTENT(out) :: pd_crit !< critical diameter of clusters (m)5720 REAL(wp), INTENT(out) :: pk_ocnv !< if pk_ocnv = 1, organic compounds participate nucleation5721 REAL(wp), INTENT(out) :: pk_sa !< if pk_sa = 1, H2SO4 participate in nucleation5722 REAL(wp), INTENT(out) :: pn_crit_ocnv !< number of organic molecules in cluster (#)5723 REAL(wp), INTENT(out) :: pn_crit_sa !< number of H2SO4 molecules in cluster (#)5724 REAL(wp), INTENT(out) :: pnuc_rate !< nucleation rate (#/(m3 s))5725 !5726 !-- Nucleation rate (#/m3/s)5727 pnuc_rate = ( z_k_sa1 * pc_sa**2 + z_k_sa2 * pc_sa * pc_org ) * 1.0E+6_wp5728 !5729 !-- Both organic compounds and H2SO4 are involved when SAnucleation is assumed.5730 pn_crit_sa = 3.0_wp5731 pn_crit_ocnv = 1.0_wp5732 pk_sa = 1.0_wp5733 pk_ocnv = 1.0_wp5734 pd_crit = 1.5E-9_wp ! (m)5735 5736 END SUBROUTINE SAnucl5737 5738 !------------------------------------------------------------------------------!5739 ! Description:5740 ! ------------5741 !> Takes into account the homomolecular nucleation of both sulphuric acid and5742 !> Lorganic with heteromolecular nucleation.5743 !> See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.5744 !------------------------------------------------------------------------------!5745 SUBROUTINE SAORGnucl( pc_sa, pc_org, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv )5746 5747 IMPLICIT NONE5748 5749 REAL(wp) :: z_k_s1 = 1.4E-14_wp !< (cm3/s])5750 REAL(wp) :: z_k_s2 = 2.6E-14_wp !< (cm3/s])5751 REAL(wp) :: z_k_s3 = 0.037E-14_wp !< (cm3/s]) (Paasonen et al. Table 3.)5752 5753 REAL(wp), INTENT(in) :: pc_org !< organic vapour concentration (#/m3)5754 REAL(wp), INTENT(in) :: pc_sa !< H2SO4 conc. (#/m3)5755 5756 REAL(wp), INTENT(out) :: pd_crit !< critical diameter of clusters (m)5757 REAL(wp), INTENT(out) :: pk_ocnv !< if pk_ocnv = 1, organic compounds participate in nucleation5758 REAL(wp), INTENT(out) :: pk_sa !< if pk_sa = 1, H2SO4 participate in nucleation5759 REAL(wp), INTENT(out) :: pn_crit_ocnv !< number of organic molecules in cluster (#)5760 REAL(wp), INTENT(out) :: pn_crit_sa !< number of H2SO4 molecules in cluster (#)5761 REAL(wp), INTENT(out) :: pnuc_rate !< nucl. rate (#/(m3 s))5762 !5763 !-- Nucleation rate (#/m3/s)5764 pnuc_rate = ( z_k_s1 * pc_sa**2 + z_k_s2 * pc_sa * pc_org + z_k_s3 * pc_org**2 ) * 1.0E+6_wp5765 !5766 !-- Organic compounds not involved when kinetic nucleation is assumed.5767 pn_crit_sa = 3.0_wp5768 pn_crit_ocnv = 3.0_wp5769 pk_sa = 1.0_wp5770 pk_ocnv = 1.0_wp5771 pd_crit = 1.5E-9_wp ! (m)5772 5773 END SUBROUTINE SAORGnucl5774 5630 5775 5631 !------------------------------------------------------------------------------! … … 8237 8093 8238 8094 DO ib = 1, nbins_aerosol 8239 aerosol_number(ib)%conc_p(nzt+1,:,:) = aerosol_number(ib)%conc_p(nzt,:,:) + 8095 aerosol_number(ib)%conc_p(nzt+1,:,:) = aerosol_number(ib)%conc_p(nzt,:,:) + & 8240 8096 bc_an_t_val(ib) * dzu(nzt+1) 8241 8097 DO ic = 1, ncomponents_mass 8242 8098 icc = ( ic - 1 ) * nbins_aerosol + ib 8243 aerosol_mass(icc)%conc_p(nzt+1,:,:) = aerosol_mass(icc)%conc_p(nzt,:,:) + 8099 aerosol_mass(icc)%conc_p(nzt+1,:,:) = aerosol_mass(icc)%conc_p(nzt,:,:) + & 8244 8100 bc_am_t_val(icc) * dzu(nzt+1) 8245 8101 ENDDO … … 8247 8103 IF ( .NOT. salsa_gases_from_chem ) THEN 8248 8104 DO ig = 1, ngases_salsa 8249 salsa_gas(ig)%conc_p(nzt+1,:,:) = salsa_gas(ig)%conc_p(nzt,:,:) + 8105 salsa_gas(ig)%conc_p(nzt+1,:,:) = salsa_gas(ig)%conc_p(nzt,:,:) + & 8250 8106 bc_gt_t_val(ig) * dzu(nzt+1) 8251 8107 ENDDO … … 8550 8406 ONLY: check_existence, close_input_file, get_attribute, get_variable, & 8551 8407 inquire_num_variables, inquire_variable_names, & 8552 netcdf_data_input_get_dimension_length, open_read_file 8408 netcdf_data_input_get_dimension_length, open_read_file, street_type_f 8553 8409 8554 8410 USE surface_mod, & … … 8561 8417 CHARACTER(LEN=25) :: mod_name !< name in the input file 8562 8418 8419 INTEGER(iwp) :: i !< loop index 8563 8420 INTEGER(iwp) :: ib !< loop index: aerosol number bins 8564 8421 INTEGER(iwp) :: ic !< loop index: aerosol chemical components … … 8566 8423 INTEGER(iwp) :: in !< loop index: emission category 8567 8424 INTEGER(iwp) :: inn !< loop index 8425 INTEGER(iwp) :: j !< loop index 8568 8426 INTEGER(iwp) :: ss !< loop index 8569 8427 … … 8579 8437 8580 8438 ! 8581 !-- Set source arrays to zero:8582 DO ib = 1, nbins_aerosol8583 aerosol_number(ib)%source = 0.0_wp8584 ENDDO8585 8586 DO ic = 1, ncomponents_mass * nbins_aerosol8587 aerosol_mass(ic)%source = 0.0_wp8588 ENDDO8589 8590 !8591 8439 !-- Define emissions: 8592 8593 8440 SELECT CASE ( salsa_emission_mode ) 8594 8441 8595 CASE ( 'uniform' )8442 CASE ( 'uniform', 'parameterized' ) 8596 8443 8597 8444 IF ( init ) THEN ! Do only once … … 8605 8452 CALL size_distribution( surface_aerosol_flux, aerosol_flux_dpg, aerosol_flux_sigmag, & 8606 8453 nsect_emission ) 8607 DO ib = 1, nbins_aerosol 8608 source_array(:,:,ib) = nsect_emission(ib) 8609 ENDDO 8454 IF ( salsa_emission_mode == 'uniform' ) THEN 8455 DO ib = 1, nbins_aerosol 8456 source_array(:,:,ib) = nsect_emission(ib) 8457 ENDDO 8458 ELSE 8459 IF ( street_type_f%from_file ) THEN 8460 DO i = nxl, nxr 8461 DO j = nys, nyn 8462 IF ( street_type_f%var(j,i) >= main_street_id .AND. & 8463 street_type_f%var(j,i) < max_street_id ) THEN 8464 source_array(j,i,:) = nsect_emission(:) * emiss_factor_main 8465 ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. & 8466 street_type_f%var(j,i) < main_street_id ) THEN 8467 source_array(j,i,:) = nsect_emission(:) * emiss_factor_side 8468 ENDIF 8469 ENDDO 8470 ENDDO 8471 ELSE 8472 WRITE( message_string, * ) 'salsa_emission_mode = "parameterized" but the '// & 8473 'street_type data is missing.' 8474 CALL message( 'salsa_emission_setup', 'PA0661', 1, 2, 0, 6, 0 ) 8475 ENDIF 8476 ENDIF 8610 8477 ! 8611 8478 !-- Check which chemical components are used … … 8622 8489 aerosol_flux_mass_fracs_a = aerosol_flux_mass_fracs_a / & 8623 8490 SUM( aerosol_flux_mass_fracs_a(1:ncc ) ) 8624 ! 8625 !-- Set uniform fluxes of default horizontal surfaces 8626 CALL set_flux( surf_def_h(0), cc_i2m, aerosol_flux_mass_fracs_a, source_array ) 8491 IF ( salsa_emission_mode == 'uniform' ) THEN 8492 ! 8493 !-- Set uniform fluxes of default horizontal surfaces 8494 CALL set_flux( surf_def_h(0), cc_i2m, aerosol_flux_mass_fracs_a, source_array ) 8495 ELSE 8496 ! 8497 !-- Set fluxes normalised based on the street type on land surfaces 8498 CALL set_flux( surf_lsm_h, cc_i2m, aerosol_flux_mass_fracs_a, source_array ) 8499 ENDIF 8627 8500 8628 8501 DEALLOCATE( nsect_emission, source_array ) 8629 8502 ENDIF 8630 8631 CASE ( 'parameterized' )8632 !8633 !-- TO DO8634 8503 8635 8504 CASE ( 'read_from_file' ) … … 8643 8512 surf_usm_h%amsws = 0.0_wp 8644 8513 8514 ! 8515 !-- Reset source arrays: 8516 DO ib = 1, nbins_aerosol 8517 aerosol_number(ib)%source = 0.0_wp 8518 ENDDO 8519 8520 DO ic = 1, ncomponents_mass * nbins_aerosol 8521 aerosol_mass(ic)%source = 0.0_wp 8522 ENDDO 8523 8645 8524 #if defined( __netcdf ) 8646 8525 ! … … 8658 8537 IF ( init ) THEN 8659 8538 ! 8539 !-- Variable names 8540 CALL inquire_num_variables( id_salsa, aero_emission_att%num_vars ) 8541 ALLOCATE( aero_emission_att%var_names(1:aero_emission_att%num_vars) ) 8542 CALL inquire_variable_names( id_salsa, aero_emission_att%var_names ) 8543 ! 8660 8544 !-- Read the index and name of chemical components 8661 8545 CALL netcdf_data_input_get_dimension_length( id_salsa, aero_emission_att%ncc, & … … 8663 8547 ALLOCATE( aero_emission_att%cc_index(1:aero_emission_att%ncc) ) 8664 8548 CALL get_variable( id_salsa, 'composition_index', aero_emission_att%cc_index ) 8665 CALL get_variable( id_salsa, 'composition_name', aero_emission_att%cc_name, & 8666 aero_emission_att%ncc ) 8549 8550 IF ( check_existence( aero_emission_att%var_names, 'composition_name' ) ) THEN 8551 CALL get_variable( id_salsa, 'composition_name', aero_emission_att%cc_name, & 8552 aero_emission_att%ncc ) 8553 ELSE 8554 message_string = 'Missing composition_name in ' // TRIM( input_file_salsa ) 8555 CALL message( 'salsa_emission_setup', 'PA0657', 1, 2, 0, 6, 0 ) 8556 ENDIF 8667 8557 ! 8668 8558 !-- Find the corresponding chemical components in the model 8669 aero_emission_att%cc_in put_to_model= 08559 aero_emission_att%cc_in2mod = 0 8670 8560 DO ic = 1, aero_emission_att%ncc 8671 8561 in_name = aero_emission_att%cc_name(ic) 8672 8562 SELECT CASE ( TRIM( in_name ) ) 8673 8563 CASE ( 'H2SO4', 'h2so4', 'SO4', 'so4' ) 8674 aero_emission_att%cc_in put_to_model(1) = ic8564 aero_emission_att%cc_in2mod(1) = ic 8675 8565 CASE ( 'OC', 'oc', 'organics' ) 8676 aero_emission_att%cc_in put_to_model(2) = ic8566 aero_emission_att%cc_in2mod(2) = ic 8677 8567 CASE ( 'BC', 'bc' ) 8678 aero_emission_att%cc_in put_to_model(3) = ic8568 aero_emission_att%cc_in2mod(3) = ic 8679 8569 CASE ( 'DU', 'du' ) 8680 aero_emission_att%cc_in put_to_model(4) = ic8570 aero_emission_att%cc_in2mod(4) = ic 8681 8571 CASE ( 'SS', 'ss' ) 8682 aero_emission_att%cc_in put_to_model(5) = ic8572 aero_emission_att%cc_in2mod(5) = ic 8683 8573 CASE ( 'HNO3', 'hno3', 'NO', 'no' ) 8684 aero_emission_att%cc_in put_to_model(6) = ic8574 aero_emission_att%cc_in2mod(6) = ic 8685 8575 CASE ( 'NH3', 'nh3', 'NH', 'nh' ) 8686 aero_emission_att%cc_in put_to_model(7) = ic8576 aero_emission_att%cc_in2mod(7) = ic 8687 8577 END SELECT 8688 8578 8689 8579 ENDDO 8690 8580 8691 IF ( SUM( aero_emission_att%cc_in put_to_model) == 0 ) THEN8581 IF ( SUM( aero_emission_att%cc_in2mod ) == 0 ) THEN 8692 8582 message_string = 'None of the aerosol chemical components in ' // TRIM( & 8693 8583 input_file_salsa ) // ' correspond to the ones applied in SALSA.' … … 8706 8596 CALL get_attribute( id_salsa, 'lod', aero_emission_att%lod, .FALSE., & 8707 8597 'aerosol_emission_values' ) 8708 !8709 !-- Variable names8710 CALL inquire_num_variables( id_salsa, aero_emission_att%num_vars )8711 ALLOCATE( aero_emission_att%var_names(1:aero_emission_att%num_vars) )8712 CALL inquire_variable_names( id_salsa, aero_emission_att%var_names )8713 8598 8714 8599 ! … … 8738 8623 ! 8739 8624 !-- Get emission category names and indices 8740 CALL get_variable( id_salsa, 'emission_category_name', aero_emission_att%cat_name, & 8741 aero_emission_att%ncat) 8625 IF ( check_existence( aero_emission_att%var_names, 'emission_category_name' ) ) THEN 8626 CALL get_variable( id_salsa, 'emission_category_name', & 8627 aero_emission_att%cat_name, aero_emission_att%ncat ) 8628 ELSE 8629 message_string = 'Missing emission_category_name in ' // TRIM( input_file_salsa ) 8630 CALL message( 'salsa_emission_setup', 'PA0658', 1, 2, 0, 6, 0 ) 8631 ENDIF 8742 8632 CALL get_variable( id_salsa, 'emission_category_index', aero_emission_att%cat_index ) 8743 8633 ! … … 8764 8654 !-- For each hour of year: 8765 8655 IF ( check_existence( aero_emission_att%var_names, 'nhoursyear' ) ) THEN 8766 CALL netcdf_data_input_get_dimension_length( id_salsa, aero_emission_att%nhoursyear,&8767 8656 CALL netcdf_data_input_get_dimension_length( id_salsa, & 8657 aero_emission_att%nhoursyear, 'nhoursyear' ) 8768 8658 ALLOCATE( aero_emission_att%etf(1:aero_emission_att%ncat, & 8769 8659 1:aero_emission_att%nhoursyear) ) … … 8790 8680 ! 8791 8681 !-- Get chemical composition (i.e. mass fraction of different species) in aerosols 8792 ALLOCATE( aero_emission%def_mass_fracs(1:aero_emission_att%ncat, & 8793 1:aero_emission_att%ncc) ) 8794 aero_emission%def_mass_fracs = 0.0_wp 8795 CALL get_variable( id_salsa, 'emission_mass_fracs', aero_emission%def_mass_fracs, & 8796 0, aero_emission_att%ncc-1, 0, aero_emission_att%ncat-1 ) 8682 IF ( check_existence( aero_emission_att%var_names, 'emission_mass_fracs' ) ) THEN 8683 ALLOCATE( aero_emission%def_mass_fracs(1:aero_emission_att%ncat, & 8684 1:aero_emission_att%ncc) ) 8685 aero_emission%def_mass_fracs = 0.0_wp 8686 CALL get_variable( id_salsa, 'emission_mass_fracs', aero_emission%def_mass_fracs,& 8687 0, aero_emission_att%ncc-1, 0, aero_emission_att%ncat-1 ) 8688 ELSE 8689 message_string = 'Missing emission_mass_fracs in ' // TRIM( input_file_salsa ) 8690 CALL message( 'salsa_emission_setup', 'PA0659', 1, 2, 0, 6, 0 ) 8691 ENDIF 8797 8692 ! 8798 8693 !-- If the chemical component is not activated, set its mass fraction to 0 to avoid 8799 8694 !-- inbalance between number and mass flux 8800 cc_i2m = aero_emission_att%cc_in put_to_model8801 IF ( index_so4 < 0 .AND. cc_i2m(1) /= 0 )&8695 cc_i2m = aero_emission_att%cc_in2mod 8696 IF ( index_so4 < 0 .AND. cc_i2m(1) > 0 ) & 8802 8697 aero_emission%def_mass_fracs(:,cc_i2m(1)) = 0.0_wp 8803 IF ( index_oc < 0 .AND. cc_i2m(2) /= 0 )&8698 IF ( index_oc < 0 .AND. cc_i2m(2) > 0 ) & 8804 8699 aero_emission%def_mass_fracs(:,cc_i2m(2)) = 0.0_wp 8805 IF ( index_bc < 0 .AND. cc_i2m(3) /= 0 )&8700 IF ( index_bc < 0 .AND. cc_i2m(3) > 0 ) & 8806 8701 aero_emission%def_mass_fracs(:,cc_i2m(3)) = 0.0_wp 8807 IF ( index_du < 0 .AND. cc_i2m(4) /= 0 )&8702 IF ( index_du < 0 .AND. cc_i2m(4) > 0 ) & 8808 8703 aero_emission%def_mass_fracs(:,cc_i2m(4)) = 0.0_wp 8809 IF ( index_ss < 0 .AND. cc_i2m(5) /= 0 )&8704 IF ( index_ss < 0 .AND. cc_i2m(5) > 0 ) & 8810 8705 aero_emission%def_mass_fracs(:,cc_i2m(5)) = 0.0_wp 8811 IF ( index_no < 0 .AND. cc_i2m(6) /= 0 )&8706 IF ( index_no < 0 .AND. cc_i2m(6) > 0 ) & 8812 8707 aero_emission%def_mass_fracs(:,cc_i2m(6)) = 0.0_wp 8813 IF ( index_nh < 0 .AND. cc_i2m(7) /= 0 )&8708 IF ( index_nh < 0 .AND. cc_i2m(7) > 0 ) & 8814 8709 aero_emission%def_mass_fracs(:,cc_i2m(7)) = 0.0_wp 8815 8710 ! … … 8886 8781 ! 8887 8782 !-- Read time stamps: 8888 CALL get_variable( id_salsa, 'time', aero_emission_att%time ) 8783 IF ( check_existence( aero_emission_att%var_names, 'time' ) ) THEN 8784 CALL get_variable( id_salsa, 'time', aero_emission_att%time ) 8785 ELSE 8786 message_string = 'Missing time in ' // TRIM( input_file_salsa ) 8787 CALL message( 'salsa_emission_setup', 'PA0660', 1, 2, 0, 6, 0 ) 8788 ENDIF 8889 8789 ! 8890 8790 !-- Read emission mass fractions 8891 CALL get_variable( id_salsa, 'emission_mass_fracs', aero_emission%preproc_mass_fracs ) 8791 IF ( check_existence( aero_emission_att%var_names, 'emission_mass_fracs' ) ) THEN 8792 CALL get_variable( id_salsa, 'emission_mass_fracs', & 8793 aero_emission%preproc_mass_fracs ) 8794 ELSE 8795 message_string = 'Missing emission_mass_fracs in ' // TRIM( input_file_salsa ) 8796 CALL message( 'salsa_emission_setup', 'PA0659', 1, 2, 0, 6, 0 ) 8797 ENDIF 8892 8798 ! 8893 8799 !-- If the chemical component is not activated, set its mass fraction to 0 8894 cc_i2m = aero_emission_att%cc_in put_to_model8800 cc_i2m = aero_emission_att%cc_in2mod 8895 8801 IF ( index_so4 < 0 .AND. cc_i2m(1) /= 0 ) & 8896 8802 aero_emission%preproc_mass_fracs(cc_i2m(1)) = 0.0_wp … … 8970 8876 !-- only for either default, land or urban surface. 8971 8877 IF ( .NOT. land_surface .AND. .NOT. urban_surface ) THEN 8972 CALL set_flux( surf_def_h(0), aero_emission_att%cc_in put_to_model,&8878 CALL set_flux( surf_def_h(0), aero_emission_att%cc_in2mod, & 8973 8879 aero_emission%def_mass_fracs(in,:), source_array ) 8974 8880 ELSE 8975 CALL set_flux( surf_usm_h, aero_emission_att%cc_in put_to_model,&8881 CALL set_flux( surf_usm_h, aero_emission_att%cc_in2mod, & 8976 8882 aero_emission%def_mass_fracs(in,:), source_array ) 8977 CALL set_flux( surf_lsm_h, aero_emission_att%cc_in put_to_model,&8883 CALL set_flux( surf_lsm_h, aero_emission_att%cc_in2mod, & 8978 8884 aero_emission%def_mass_fracs(in,:), source_array ) 8979 8885 ENDIF … … 9004 8910 !-- for either default, land and urban surface. 9005 8911 IF ( .NOT. land_surface .AND. .NOT. urban_surface ) THEN 9006 CALL set_flux( surf_def_h(0), aero_emission_att%cc_in put_to_model,&8912 CALL set_flux( surf_def_h(0), aero_emission_att%cc_in2mod, & 9007 8913 aero_emission%preproc_mass_fracs, aero_emission%preproc_data ) 9008 8914 ELSE 9009 CALL set_flux( surf_usm_h, aero_emission_att%cc_in put_to_model,&8915 CALL set_flux( surf_usm_h, aero_emission_att%cc_in2mod, & 9010 8916 aero_emission%preproc_mass_fracs, aero_emission%preproc_data ) 9011 CALL set_flux( surf_lsm_h, aero_emission_att%cc_in put_to_model,&8917 CALL set_flux( surf_lsm_h, aero_emission_att%cc_in2mod, & 9012 8918 aero_emission%preproc_mass_fracs, aero_emission%preproc_data ) 9013 8919 ENDIF … … 9249 9155 9250 9156 INTEGER(iwp) :: id_chem !< NetCDF id of chemistry emission file 9157 INTEGER(iwp) :: i !< loop index 9251 9158 INTEGER(iwp) :: ig !< loop index 9252 9159 INTEGER(iwp) :: in !< running index for emission categories 9160 INTEGER(iwp) :: j !< loop index 9253 9161 INTEGER(iwp) :: num_vars !< number of variables 9254 9162 … … 9258 9166 9259 9167 REAL(wp), DIMENSION(:), ALLOCATABLE :: time_factor !< emission time factor 9168 9169 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dum_var_3d !< 9170 9171 REAL(wp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: dum_var_5d !< 9260 9172 9261 9173 ! … … 9280 9192 ! 9281 9193 !-- Read the index and name of chemical components 9282 CALL netcdf_data_input_get_dimension_length( id_chem, chem_emission_att%n spec,&9194 CALL netcdf_data_input_get_dimension_length( id_chem, chem_emission_att%n_emiss_species, & 9283 9195 'nspecies' ) 9284 ALLOCATE( chem_emission_att%species_index(1:chem_emission_att%n spec) )9196 ALLOCATE( chem_emission_att%species_index(1:chem_emission_att%n_emiss_species) ) 9285 9197 CALL get_variable( id_chem, 'emission_index', chem_emission_att%species_index ) 9286 9198 CALL get_variable( id_chem, 'emission_name', chem_emission_att%species_name, & 9287 chem_emission_att%nspec ) 9199 chem_emission_att%n_emiss_species ) 9200 ! 9201 !-- Allocate emission data 9202 ALLOCATE( chem_emission(1:chem_emission_att%n_emiss_species) ) 9288 9203 ! 9289 9204 !-- Find the corresponding indices in the model 9290 9205 emission_index_chem = 0 9291 DO ig = 1, chem_emission_att%n spec9206 DO ig = 1, chem_emission_att%n_emiss_species 9292 9207 in_name = chem_emission_att%species_name(ig) 9293 9208 SELECT CASE ( TRIM( in_name ) ) … … 9366 9281 !-- Allocate and read surface emission data (in total PM) (NOTE that "preprocessed" input data 9367 9282 !-- array is applied now here) 9368 ALLOCATE( chem_emission%preproc_emission_data(nys:nyn,nxl:nxr, 1:chem_emission_att%nspec,& 9369 1:chem_emission_att%ncat) ) 9370 CALL get_variable( id_chem, 'emission_values', chem_emission%preproc_emission_data, & 9371 0, chem_emission_att%ncat-1, 0, chem_emission_att%nspec-1, & 9372 nxl, nxr, nys, nyn ) 9283 ALLOCATE( dum_var_5d(1,nys:nyn,nxl:nxr,1:chem_emission_att%n_emiss_species, & 9284 1:chem_emission_att%ncat) ) 9285 CALL get_variable( id_chem, 'emission_values', dum_var_5d, 0, chem_emission_att%ncat-1, & 9286 0, chem_emission_att%n_emiss_species-1, nxl, nxr, nys, nyn, 0, 0 ) 9287 DO ig = 1, chem_emission_att%n_emiss_species 9288 ALLOCATE( chem_emission(ig)%default_emission_data(nys:nyn,nxl:nxr, & 9289 1:chem_emission_att%ncat) ) 9290 DO in = 1, chem_emission_att%ncat 9291 DO i = nxl, nxr 9292 DO j = nys, nyn 9293 chem_emission(ig)%default_emission_data(j,i,in) = dum_var_5d(1,j,i,ig,in) 9294 ENDDO 9295 ENDDO 9296 ENDDO 9297 ENDDO 9298 DEALLOCATE( dum_var_5d ) 9373 9299 ! 9374 9300 !-- Pre-processed mode: … … 9412 9338 ! 9413 9339 !-- Set gas emissions for each emission category 9340 ALLOCATE( dum_var_3d(nys:nyn,nxl:nxr,1:chem_emission_att%n_emiss_species) ) 9341 9414 9342 DO in = 1, chem_emission_att%ncat 9343 DO ig = 1, chem_emission_att%n_emiss_species 9344 dum_var_3d(:,:,ig) = chem_emission(ig)%default_emission_data(:,:,in) 9345 ENDDO 9415 9346 ! 9416 9347 !-- Set surface fluxes only for either default, land or urban surface 9417 9348 IF ( .NOT. land_surface .AND. .NOT. urban_surface ) THEN 9418 9349 CALL set_gas_flux( surf_def_h(0), emission_index_chem, chem_emission_att%units, & 9419 chem_emission%preproc_emission_data(:,:,:,in), time_factor(in) )9350 dum_var_3d(:,:,in), time_factor(in) ) 9420 9351 ELSE 9421 9352 CALL set_gas_flux( surf_usm_h, emission_index_chem, chem_emission_att%units, & 9422 chem_emission%preproc_emission_data(:,:,:,in), time_factor(in) )9353 dum_var_3d(:,:,in), time_factor(in) ) 9423 9354 CALL set_gas_flux( surf_lsm_h, emission_index_chem, chem_emission_att%units, & 9424 chem_emission%preproc_emission_data(:,:,:,in), time_factor(in) )9355 dum_var_3d(:,:,in), time_factor(in) ) 9425 9356 ENDIF 9426 9357 ENDDO 9358 DEALLOCATE( dum_var_3d ) 9427 9359 ! 9428 9360 !-- The next emission update is again after one hour … … 9438 9370 !-- Allocate the data input array always before reading in the data and deallocate after (NOTE 9439 9371 !-- that "preprocessed" input data array is applied now here) 9440 ALLOCATE( chem_emission%default_emission_data(nys:nyn,nxl:nxr,1:nbins_aerosol) )9372 ALLOCATE( dum_var_5d(1,1,nys:nyn,nxl:nxr,1:chem_emission_att%n_emiss_species) ) 9441 9373 ! 9442 9374 !-- Read in the next time step 9443 CALL get_variable( id_chem, 'emission_values', chem_emission%default_emission_data, & 9444 chem_emission_att%i_hour, 0, chem_emission_att%nspec-1, nxl, nxr, nys, nyn ) 9375 CALL get_variable( id_chem, 'emission_values', dum_var_5d, & 9376 0, chem_emission_att%n_emiss_species-1, nxl, nxr, nys, nyn, 0, 0, & 9377 chem_emission_att%i_hour, chem_emission_att%i_hour ) 9445 9378 ! 9446 9379 !-- Set surface fluxes only for either default, land or urban surface 9447 9380 IF ( .NOT. land_surface .AND. .NOT. urban_surface ) THEN 9448 9381 CALL set_gas_flux( surf_def_h(0), emission_index_chem, chem_emission_att%units, & 9449 chem_emission%default_emission_data)9382 dum_var_5d(1,1,:,:,:) ) 9450 9383 ELSE 9451 9384 CALL set_gas_flux( surf_usm_h, emission_index_chem, chem_emission_att%units, & 9452 chem_emission%default_emission_data)9385 dum_var_5d(1,1,:,:,:) ) 9453 9386 CALL set_gas_flux( surf_lsm_h, emission_index_chem, chem_emission_att%units, & 9454 chem_emission%default_emission_data ) 9455 ENDIF 9387 dum_var_5d(1,1,:,:,:) ) 9388 ENDIF 9389 DEALLOCATE ( dum_var_5d ) 9456 9390 ! 9457 9391 !-- Determine the next emission update 9458 9392 next_gas_emission_update = gas_emission_time(chem_emission_att%i_hour+2) 9459 9393 9460 DEALLOCATE( chem_emission%default_emission_data )9461 9394 ENDIF 9462 9395 ! … … 9506 9439 REAL(wp), DIMENSION(ngases_salsa) :: conv !< unit conversion factor 9507 9440 9508 REAL(wp), DIMENSION(nys:nyn,nxl:nxr,chem_emission_att%n spec), INTENT(in) :: source_array !<9441 REAL(wp), DIMENSION(nys:nyn,nxl:nxr,chem_emission_att%n_emiss_species), INTENT(in) :: source_array !< 9509 9442 9510 9443 TYPE(surf_type), INTENT(inout) :: surface !< respective surface type … … 9581 9514 CHARACTER(LEN=*) :: var !< 9582 9515 9583 SELECT CASE ( TRIM( var ) ) 9584 9585 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 9586 IF ( .NOT. salsa ) THEN 9587 message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.' 9588 CALL message( 'check_parameters', 'PA0652', 1, 2, 0, 6, 0 ) 9589 ENDIF 9590 IF ( salsa_gases_from_chem ) THEN 9591 message_string = 'gases are imported from the chemistry module and thus output of "' & 9592 // TRIM( var ) // '" is not allowed' 9593 CALL message( 'check_parameters', 'PA0653', 1, 2, 0, 6, 0 ) 9594 ENDIF 9516 INTEGER(iwp) :: char_to_int !< for converting character to integer 9517 9518 IF ( var(1:6) /= 'salsa_' ) THEN 9519 unit = 'illegal' 9520 RETURN 9521 ENDIF 9522 ! 9523 !-- Treat bin-specific outputs separately 9524 IF ( var(7:11) == 'N_bin' ) THEN 9525 READ( var(12:),* ) char_to_int 9526 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 9595 9527 unit = '#/m3' 9596 9597 CASE ( 'LDSA' ) 9598 IF ( .NOT. salsa ) THEN 9599 message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.' 9600 CALL message( 'check_parameters', 'PA0646', 1, 2, 0, 6, 0 ) 9601 ENDIF 9602 unit = 'mum2/cm3' 9603 9604 CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4', 'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8', & 9605 'm_bin9', 'm_bin10', 'm_bin11', 'm_bin12', 'PM2.5', 'PM10', 's_BC', 's_DU', & 9606 's_H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 9607 IF ( .NOT. salsa ) THEN 9608 message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.' 9609 CALL message( 'check_parameters', 'PA0647', 1, 2, 0, 6, 0 ) 9610 ENDIF 9528 ELSE 9529 unit = 'illegal' 9530 RETURN 9531 ENDIF 9532 9533 ELSEIF ( var(7:11) == 'm_bin' ) THEN 9534 READ( var(12:),* ) char_to_int 9535 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 9611 9536 unit = 'kg/m3' 9612 9613 CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4', 'N_bin5', 'N_bin6', 'N_bin7', 'N_bin8', & 9614 'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12', 'Ntot' ) 9615 IF ( .NOT. salsa ) THEN 9616 message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.' 9617 CALL message( 'check_parameters', 'PA0645', 1, 2, 0, 6, 0 ) 9618 ENDIF 9619 unit = '#/m3' 9537 ELSE 9538 unit = 'illegal' 9539 RETURN 9540 ENDIF 9541 9542 ELSE 9543 SELECT CASE ( TRIM( var(7:) ) ) 9544 9545 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 9546 IF ( salsa_gases_from_chem ) THEN 9547 message_string = 'gases are imported from the chemistry module and thus output '// & 9548 'of "' // TRIM( var ) // '" is not allowed' 9549 CALL message( 'check_parameters', 'PA0653', 1, 2, 0, 6, 0 ) 9550 ENDIF 9551 unit = '#/m3' 9552 9553 CASE ( 'LDSA' ) 9554 unit = 'mum2/cm3' 9555 9556 CASE ( 'PM0.1', 'PM2.5', 'PM10', 's_BC', 's_DU', 's_H2O', 's_NH', 's_NO', 's_OC', & 9557 's_SO4', 's_SS' ) 9558 unit = 'kg/m3' 9559 9560 CASE ( 'N_UFP', 'Ntot' ) 9561 unit = '#/m3' 9562 9563 CASE DEFAULT 9564 unit = 'illegal' 9565 9566 END SELECT 9567 ENDIF 9568 9569 END SUBROUTINE salsa_check_data_output 9570 9571 !------------------------------------------------------------------------------! 9572 ! Description: 9573 ! ------------ 9574 !> Check profile data output for salsa. Currently only for diagnostic variables 9575 !> Ntot, N_UFP, PM0.1, PM2.5, PM10 and LDSA 9576 !------------------------------------------------------------------------------! 9577 SUBROUTINE salsa_check_data_output_pr( var, var_count, unit, dopr_unit ) 9578 9579 USE arrays_3d, & 9580 ONLY: zu 9581 9582 USE profil_parameter, & 9583 ONLY: dopr_index 9584 9585 USE statistics, & 9586 ONLY: hom, pr_palm, statistic_regions 9587 9588 IMPLICIT NONE 9589 9590 CHARACTER(LEN=*) :: dopr_unit !< 9591 CHARACTER(LEN=*) :: unit !< 9592 CHARACTER(LEN=*) :: var !< 9593 9594 INTEGER(iwp) :: var_count !< 9595 9596 IF ( var(1:6) /= 'salsa_' ) THEN 9597 unit = 'illegal' 9598 RETURN 9599 ENDIF 9600 9601 SELECT CASE ( TRIM( var(7:) ) ) 9602 9603 CASE( 'LDSA' ) 9604 salsa_pr_count = salsa_pr_count + 1 9605 salsa_pr_index(salsa_pr_count) = 1 9606 dopr_index(var_count) = pr_palm + salsa_pr_count 9607 dopr_unit = 'mum2/cm3' 9608 unit = dopr_unit 9609 hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 ) 9610 9611 CASE( 'N_UFP' ) 9612 salsa_pr_count = salsa_pr_count + 1 9613 salsa_pr_index(salsa_pr_count) = 2 9614 dopr_index(var_count) = pr_palm + salsa_pr_count 9615 dopr_unit = '#/m3' 9616 unit = dopr_unit 9617 hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 ) 9618 9619 CASE( 'Ntot' ) 9620 salsa_pr_count = salsa_pr_count + 1 9621 salsa_pr_index(salsa_pr_count) = 3 9622 dopr_index(var_count) = pr_palm + salsa_pr_count 9623 dopr_unit = '#/m3' 9624 unit = dopr_unit 9625 hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 ) 9626 9627 CASE( 'PM0.1' ) 9628 salsa_pr_count = salsa_pr_count + 1 9629 salsa_pr_index(salsa_pr_count) = 4 9630 dopr_index(var_count) = pr_palm + salsa_pr_count 9631 dopr_unit = 'kg/m3' 9632 unit = dopr_unit 9633 hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 ) 9634 9635 CASE( 'PM2.5' ) 9636 salsa_pr_count = salsa_pr_count + 1 9637 salsa_pr_index(salsa_pr_count) = 5 9638 dopr_index(var_count) = pr_palm + salsa_pr_count 9639 dopr_unit = 'kg/m3' 9640 unit = dopr_unit 9641 hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 ) 9642 9643 CASE( 'PM10' ) 9644 salsa_pr_count = salsa_pr_count + 1 9645 salsa_pr_index(salsa_pr_count) = 6 9646 dopr_index(var_count) = pr_palm + salsa_pr_count 9647 dopr_unit = 'kg/m3' 9648 unit = dopr_unit 9649 hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 ) 9620 9650 9621 9651 CASE DEFAULT … … 9624 9654 END SELECT 9625 9655 9626 END SUBROUTINE salsa_check_data_output 9656 9657 END SUBROUTINE salsa_check_data_output_pr 9658 9659 !-------------------------------------------------------------------------------! 9660 !> Description: 9661 !> Calculation of horizontally averaged profiles for salsa. 9662 !-------------------------------------------------------------------------------! 9663 SUBROUTINE salsa_statistics( mode, sr, tn ) 9664 9665 USE control_parameters, & 9666 ONLY: max_pr_user 9667 9668 USE chem_modules, & 9669 ONLY: max_pr_cs 9670 9671 USE statistics, & 9672 ONLY: pr_palm, rmask, sums_l 9673 9674 IMPLICIT NONE 9675 9676 CHARACTER(LEN=*) :: mode !< 9677 9678 INTEGER(iwp) :: i !< loop index 9679 INTEGER(iwp) :: ib !< loop index 9680 INTEGER(iwp) :: ic !< loop index 9681 INTEGER(iwp) :: ii !< loop index 9682 INTEGER(iwp) :: ind !< index in the statistical output 9683 INTEGER(iwp) :: j !< loop index 9684 INTEGER(iwp) :: k !< loop index 9685 INTEGER(iwp) :: sr !< statistical region 9686 INTEGER(iwp) :: tn !< thread number 9687 9688 REAL(wp) :: df !< For calculating LDSA: fraction of particles depositing in the alveolar 9689 !< (or tracheobronchial) region of the lung. Depends on the particle size 9690 REAL(wp) :: mean_d !< Particle diameter in micrometres 9691 REAL(wp) :: temp_bin !< temporary variable 9692 9693 IF ( mode == 'profiles' ) THEN 9694 !$OMP DO 9695 DO ii = 1, salsa_pr_count 9696 9697 ind = pr_palm + max_pr_user + max_pr_cs + ii 9698 9699 SELECT CASE( salsa_pr_index(ii) ) 9700 9701 CASE( 1 ) ! LDSA 9702 DO i = nxl, nxr 9703 DO j = nys, nyn 9704 DO k = nzb, nzt+1 9705 temp_bin = 0.0_wp 9706 DO ib = 1, nbins_aerosol 9707 ! 9708 !-- Diameter in micrometres 9709 mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp 9710 ! 9711 !-- Deposition factor: alveolar 9712 df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) + & 9713 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) - & 9714 1.362_wp )**2 ) ) 9715 ! 9716 !-- Lung-deposited surface area LDSA (units mum2/cm3) 9717 temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp * & 9718 aerosol_number(ib)%conc(k,j,i) 9719 ENDDO 9720 sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr) * & 9721 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 9722 ENDDO 9723 ENDDO 9724 ENDDO 9725 9726 CASE( 2 ) ! N_UFP 9727 DO i = nxl, nxr 9728 DO j = nys, nyn 9729 DO k = nzb, nzt+1 9730 temp_bin = 0.0_wp 9731 DO ib = 1, nbins_aerosol 9732 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp ) & 9733 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 9734 ENDDO 9735 sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr) * & 9736 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 9737 ENDDO 9738 ENDDO 9739 ENDDO 9740 9741 CASE( 3 ) ! Ntot 9742 DO i = nxl, nxr 9743 DO j = nys, nyn 9744 DO k = nzb, nzt+1 9745 temp_bin = 0.0_wp 9746 DO ib = 1, nbins_aerosol 9747 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 9748 ENDDO 9749 sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr) * & 9750 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 9751 ENDDO 9752 ENDDO 9753 ENDDO 9754 9755 CASE( 4 ) ! PM0.1 9756 DO i = nxl, nxr 9757 DO j = nys, nyn 9758 DO k = nzb, nzt+1 9759 temp_bin = 0.0_wp 9760 DO ib = 1, nbins_aerosol 9761 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp ) THEN 9762 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 9763 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 9764 ENDDO 9765 ENDIF 9766 ENDDO 9767 sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr) * & 9768 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 9769 ENDDO 9770 ENDDO 9771 ENDDO 9772 9773 CASE( 5 ) ! PM2.5 9774 DO i = nxl, nxr 9775 DO j = nys, nyn 9776 DO k = nzb, nzt+1 9777 temp_bin = 0.0_wp 9778 DO ib = 1, nbins_aerosol 9779 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp ) THEN 9780 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 9781 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 9782 ENDDO 9783 ENDIF 9784 ENDDO 9785 sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr) * & 9786 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 9787 ENDDO 9788 ENDDO 9789 ENDDO 9790 9791 CASE( 6 ) ! PM10 9792 DO i = nxl, nxr 9793 DO j = nys, nyn 9794 DO k = nzb, nzt+1 9795 temp_bin = 0.0_wp 9796 DO ib = 1, nbins_aerosol 9797 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp ) THEN 9798 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 9799 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 9800 ENDDO 9801 ENDIF 9802 ENDDO 9803 sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr) * & 9804 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) ) 9805 ENDDO 9806 ENDDO 9807 ENDDO 9808 9809 END SELECT 9810 ENDDO 9811 9812 ELSEIF ( mode == 'time_series' ) THEN 9813 ! 9814 !-- TODO 9815 ENDIF 9816 9817 END SUBROUTINE salsa_statistics 9818 9627 9819 9628 9820 !------------------------------------------------------------------------------! … … 9643 9835 CHARACTER(LEN=*) :: variable !< 9644 9836 9837 INTEGER(iwp) :: char_to_int !< for converting character to integer 9645 9838 INTEGER(iwp) :: found_index !< 9646 9839 INTEGER(iwp) :: i !< … … 9653 9846 !< (or tracheobronchial) region of the lung. Depends on the particle size 9654 9847 REAL(wp) :: mean_d !< Particle diameter in micrometres 9655 REAL(wp) :: nc !< Particle number concentration in units 1/cm**39656 9848 REAL(wp) :: temp_bin !< temporary variable 9657 9849 … … 9662 9854 IF ( mode == 'allocate' ) THEN 9663 9855 9664 SELECT CASE ( TRIM( variable ) ) 9665 9666 CASE ( 'g_H2SO4' ) 9667 IF ( .NOT. ALLOCATED( g_h2so4_av ) ) THEN 9668 ALLOCATE( g_h2so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9669 ENDIF 9670 g_h2so4_av = 0.0_wp 9671 9672 CASE ( 'g_HNO3' ) 9673 IF ( .NOT. ALLOCATED( g_hno3_av ) ) THEN 9674 ALLOCATE( g_hno3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9675 ENDIF 9676 g_hno3_av = 0.0_wp 9677 9678 CASE ( 'g_NH3' ) 9679 IF ( .NOT. ALLOCATED( g_nh3_av ) ) THEN 9680 ALLOCATE( g_nh3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9681 ENDIF 9682 g_nh3_av = 0.0_wp 9683 9684 CASE ( 'g_OCNV' ) 9685 IF ( .NOT. ALLOCATED( g_ocnv_av ) ) THEN 9686 ALLOCATE( g_ocnv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9687 ENDIF 9688 g_ocnv_av = 0.0_wp 9689 9690 CASE ( 'g_OCSV' ) 9691 IF ( .NOT. ALLOCATED( g_ocsv_av ) ) THEN 9692 ALLOCATE( g_ocsv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9693 ENDIF 9694 g_ocsv_av = 0.0_wp 9695 9696 CASE ( 'LDSA' ) 9697 IF ( .NOT. ALLOCATED( ldsa_av ) ) THEN 9698 ALLOCATE( ldsa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9699 ENDIF 9700 ldsa_av = 0.0_wp 9701 9702 CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4', 'N_bin5', 'N_bin6', 'N_bin7', 'N_bin8', & 9703 'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12' ) 9704 IF ( .NOT. ALLOCATED( nbins_av ) ) THEN 9705 ALLOCATE( nbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 9706 ENDIF 9707 nbins_av = 0.0_wp 9708 9709 CASE ( 'Ntot' ) 9710 IF ( .NOT. ALLOCATED( ntot_av ) ) THEN 9711 ALLOCATE( ntot_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9712 ENDIF 9713 ntot_av = 0.0_wp 9714 9715 CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4', 'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8', & 9716 'm_bin9', 'm_bin10', 'm_bin11', 'm_bin12' ) 9717 IF ( .NOT. ALLOCATED( mbins_av ) ) THEN 9718 ALLOCATE( mbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 9719 ENDIF 9720 mbins_av = 0.0_wp 9721 9722 CASE ( 'PM2.5' ) 9723 IF ( .NOT. ALLOCATED( pm25_av ) ) THEN 9724 ALLOCATE( pm25_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9725 ENDIF 9726 pm25_av = 0.0_wp 9727 9728 CASE ( 'PM10' ) 9729 IF ( .NOT. ALLOCATED( pm10_av ) ) THEN 9730 ALLOCATE( pm10_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9731 ENDIF 9732 pm10_av = 0.0_wp 9733 9734 CASE ( 's_BC' ) 9735 IF ( .NOT. ALLOCATED( s_bc_av ) ) THEN 9736 ALLOCATE( s_bc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9737 ENDIF 9738 s_bc_av = 0.0_wp 9739 9740 CASE ( 's_DU' ) 9741 IF ( .NOT. ALLOCATED( s_du_av ) ) THEN 9742 ALLOCATE( s_du_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9743 ENDIF 9744 s_du_av = 0.0_wp 9745 9746 CASE ( 's_H2O' ) 9747 IF ( .NOT. ALLOCATED( s_h2o_av ) ) THEN 9748 ALLOCATE( s_h2o_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9749 ENDIF 9750 s_h2o_av = 0.0_wp 9751 9752 CASE ( 's_NH' ) 9753 IF ( .NOT. ALLOCATED( s_nh_av ) ) THEN 9754 ALLOCATE( s_nh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9755 ENDIF 9756 s_nh_av = 0.0_wp 9757 9758 CASE ( 's_NO' ) 9759 IF ( .NOT. ALLOCATED( s_no_av ) ) THEN 9760 ALLOCATE( s_no_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9761 ENDIF 9762 s_no_av = 0.0_wp 9763 9764 CASE ( 's_OC' ) 9765 IF ( .NOT. ALLOCATED( s_oc_av ) ) THEN 9766 ALLOCATE( s_oc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9767 ENDIF 9768 s_oc_av = 0.0_wp 9769 9770 CASE ( 's_SO4' ) 9771 IF ( .NOT. ALLOCATED( s_so4_av ) ) THEN 9772 ALLOCATE( s_so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9773 ENDIF 9774 s_so4_av = 0.0_wp 9775 9776 CASE ( 's_SS' ) 9777 IF ( .NOT. ALLOCATED( s_ss_av ) ) THEN 9778 ALLOCATE( s_ss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9779 ENDIF 9780 s_ss_av = 0.0_wp 9781 9782 CASE DEFAULT 9783 CONTINUE 9784 9785 END SELECT 9856 IF ( variable(7:11) == 'N_bin' ) THEN 9857 IF ( .NOT. ALLOCATED( nbins_av ) ) THEN 9858 ALLOCATE( nbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 9859 ENDIF 9860 nbins_av = 0.0_wp 9861 9862 ELSEIF ( variable(7:11) == 'm_bin' ) THEN 9863 IF ( .NOT. ALLOCATED( mbins_av ) ) THEN 9864 ALLOCATE( mbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 9865 ENDIF 9866 mbins_av = 0.0_wp 9867 9868 ELSE 9869 9870 SELECT CASE ( TRIM( variable(7:) ) ) 9871 9872 CASE ( 'g_H2SO4' ) 9873 IF ( .NOT. ALLOCATED( g_h2so4_av ) ) THEN 9874 ALLOCATE( g_h2so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9875 ENDIF 9876 g_h2so4_av = 0.0_wp 9877 9878 CASE ( 'g_HNO3' ) 9879 IF ( .NOT. ALLOCATED( g_hno3_av ) ) THEN 9880 ALLOCATE( g_hno3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9881 ENDIF 9882 g_hno3_av = 0.0_wp 9883 9884 CASE ( 'g_NH3' ) 9885 IF ( .NOT. ALLOCATED( g_nh3_av ) ) THEN 9886 ALLOCATE( g_nh3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9887 ENDIF 9888 g_nh3_av = 0.0_wp 9889 9890 CASE ( 'g_OCNV' ) 9891 IF ( .NOT. ALLOCATED( g_ocnv_av ) ) THEN 9892 ALLOCATE( g_ocnv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9893 ENDIF 9894 g_ocnv_av = 0.0_wp 9895 9896 CASE ( 'g_OCSV' ) 9897 IF ( .NOT. ALLOCATED( g_ocsv_av ) ) THEN 9898 ALLOCATE( g_ocsv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9899 ENDIF 9900 g_ocsv_av = 0.0_wp 9901 9902 CASE ( 'LDSA' ) 9903 IF ( .NOT. ALLOCATED( ldsa_av ) ) THEN 9904 ALLOCATE( ldsa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9905 ENDIF 9906 ldsa_av = 0.0_wp 9907 9908 CASE ( 'N_UFP' ) 9909 IF ( .NOT. ALLOCATED( nufp_av ) ) THEN 9910 ALLOCATE( nufp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9911 ENDIF 9912 nufp_av = 0.0_wp 9913 9914 CASE ( 'Ntot' ) 9915 IF ( .NOT. ALLOCATED( ntot_av ) ) THEN 9916 ALLOCATE( ntot_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9917 ENDIF 9918 ntot_av = 0.0_wp 9919 9920 CASE ( 'PM0.1' ) 9921 IF ( .NOT. ALLOCATED( pm01_av ) ) THEN 9922 ALLOCATE( pm01_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9923 ENDIF 9924 pm01_av = 0.0_wp 9925 9926 CASE ( 'PM2.5' ) 9927 IF ( .NOT. ALLOCATED( pm25_av ) ) THEN 9928 ALLOCATE( pm25_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9929 ENDIF 9930 pm25_av = 0.0_wp 9931 9932 CASE ( 'PM10' ) 9933 IF ( .NOT. ALLOCATED( pm10_av ) ) THEN 9934 ALLOCATE( pm10_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9935 ENDIF 9936 pm10_av = 0.0_wp 9937 9938 CASE ( 's_BC' ) 9939 IF ( .NOT. ALLOCATED( s_bc_av ) ) THEN 9940 ALLOCATE( s_bc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9941 ENDIF 9942 s_bc_av = 0.0_wp 9943 9944 CASE ( 's_DU' ) 9945 IF ( .NOT. ALLOCATED( s_du_av ) ) THEN 9946 ALLOCATE( s_du_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9947 ENDIF 9948 s_du_av = 0.0_wp 9949 9950 CASE ( 's_H2O' ) 9951 IF ( .NOT. ALLOCATED( s_h2o_av ) ) THEN 9952 ALLOCATE( s_h2o_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9953 ENDIF 9954 s_h2o_av = 0.0_wp 9955 9956 CASE ( 's_NH' ) 9957 IF ( .NOT. ALLOCATED( s_nh_av ) ) THEN 9958 ALLOCATE( s_nh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9959 ENDIF 9960 s_nh_av = 0.0_wp 9961 9962 CASE ( 's_NO' ) 9963 IF ( .NOT. ALLOCATED( s_no_av ) ) THEN 9964 ALLOCATE( s_no_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9965 ENDIF 9966 s_no_av = 0.0_wp 9967 9968 CASE ( 's_OC' ) 9969 IF ( .NOT. ALLOCATED( s_oc_av ) ) THEN 9970 ALLOCATE( s_oc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9971 ENDIF 9972 s_oc_av = 0.0_wp 9973 9974 CASE ( 's_SO4' ) 9975 IF ( .NOT. ALLOCATED( s_so4_av ) ) THEN 9976 ALLOCATE( s_so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9977 ENDIF 9978 s_so4_av = 0.0_wp 9979 9980 CASE ( 's_SS' ) 9981 IF ( .NOT. ALLOCATED( s_ss_av ) ) THEN 9982 ALLOCATE( s_ss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 9983 ENDIF 9984 s_ss_av = 0.0_wp 9985 9986 CASE DEFAULT 9987 CONTINUE 9988 9989 END SELECT 9990 9991 ENDIF 9786 9992 9787 9993 ELSEIF ( mode == 'sum' ) THEN 9788 9994 9789 SELECT CASE ( TRIM( variable ) ) 9790 9791 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 9792 9793 vari = TRIM( variable(3:) ) 9794 9795 SELECT CASE( vari ) 9796 9797 CASE( 'H2SO4' ) 9798 found_index = 1 9799 to_be_resorted => g_h2so4_av 9800 9801 CASE( 'HNO3' ) 9802 found_index = 2 9803 to_be_resorted => g_hno3_av 9804 9805 CASE( 'NH3' ) 9806 found_index = 3 9807 to_be_resorted => g_nh3_av 9808 9809 CASE( 'OCNV' ) 9810 found_index = 4 9811 to_be_resorted => g_ocnv_av 9812 9813 CASE( 'OCSN' ) 9814 found_index = 5 9815 to_be_resorted => g_ocsv_av 9816 9817 END SELECT 9818 9995 IF ( variable(7:11) == 'N_bin' ) THEN 9996 READ( variable(12:),* ) char_to_int 9997 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 9998 ib = char_to_int 9819 9999 DO i = nxlg, nxrg 9820 10000 DO j = nysg, nyng 9821 10001 DO k = nzb, nzt+1 9822 to_be_resorted(k,j,i) = to_be_resorted(k,j,i) + & 9823 salsa_gas(found_index)%conc(k,j,i) 10002 nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) + aerosol_number(ib)%conc(k,j,i) 9824 10003 ENDDO 9825 10004 ENDDO 9826 10005 ENDDO 9827 9828 CASE ( 'LDSA' ) 10006 ENDIF 10007 10008 ELSEIF ( variable(7:11) == 'm_bin' ) THEN 10009 READ( variable(12:),* ) char_to_int 10010 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 10011 ib = char_to_int 9829 10012 DO i = nxlg, nxrg 9830 10013 DO j = nysg, nyng 9831 10014 DO k = nzb, nzt+1 9832 10015 temp_bin = 0.0_wp 9833 DO ib = 1, nbins_aerosol 9834 ! 9835 !-- Diameter in micrometres 9836 mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp 9837 ! 9838 !-- Deposition factor: alveolar (use ra_dry) 9839 df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) + & 9840 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) - & 9841 1.362_wp )**2 ) ) 9842 ! 9843 !-- Number concentration in 1/cm3 9844 nc = 1.0E-6_wp * aerosol_number(ib)%conc(k,j,i) 9845 ! 9846 !-- Lung-deposited surface area LDSA (units mum2/cm3) 9847 temp_bin = temp_bin + pi * mean_d**2 * df * nc 10016 DO ic = ib, nbins_aerosol * ncomponents_mass, nbins_aerosol 10017 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 9848 10018 ENDDO 9849 ldsa_av(k,j,i) = ldsa_av(k,j,i) + temp_bin10019 mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) + temp_bin 9850 10020 ENDDO 9851 10021 ENDDO 9852 10022 ENDDO 9853 9854 CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4', 'N_bin5', 'N_bin6', 'N_bin7', 'N_bin8', & 9855 'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12' ) 9856 DO i = nxlg, nxrg 9857 DO j = nysg, nyng 9858 DO k = nzb, nzt+1 9859 DO ib = 1, nbins_aerosol 9860 nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) + aerosol_number(ib)%conc(k,j,i) 10023 ENDIF 10024 ELSE 10025 10026 SELECT CASE ( TRIM( variable(7:) ) ) 10027 10028 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 10029 10030 vari = TRIM( variable(9:) ) ! remove salsa_g_ from beginning 10031 10032 SELECT CASE( vari ) 10033 10034 CASE( 'H2SO4' ) 10035 found_index = 1 10036 to_be_resorted => g_h2so4_av 10037 10038 CASE( 'HNO3' ) 10039 found_index = 2 10040 to_be_resorted => g_hno3_av 10041 10042 CASE( 'NH3' ) 10043 found_index = 3 10044 to_be_resorted => g_nh3_av 10045 10046 CASE( 'OCNV' ) 10047 found_index = 4 10048 to_be_resorted => g_ocnv_av 10049 10050 CASE( 'OCSV' ) 10051 found_index = 5 10052 to_be_resorted => g_ocsv_av 10053 10054 END SELECT 10055 10056 DO i = nxlg, nxrg 10057 DO j = nysg, nyng 10058 DO k = nzb, nzt+1 10059 to_be_resorted(k,j,i) = to_be_resorted(k,j,i) + & 10060 salsa_gas(found_index)%conc(k,j,i) 9861 10061 ENDDO 9862 10062 ENDDO 9863 10063 ENDDO 9864 ENDDO 9865 9866 CASE ( 'Ntot' ) 9867 DO i = nxlg, nxrg 9868 DO j = nysg, nyng 9869 DO k = nzb, nzt+1 9870 DO ib = 1, nbins_aerosol 9871 ntot_av(k,j,i) = ntot_av(k,j,i) + aerosol_number(ib)%conc(k,j,i) 10064 10065 CASE ( 'LDSA' ) 10066 DO i = nxlg, nxrg 10067 DO j = nysg, nyng 10068 DO k = nzb, nzt+1 10069 temp_bin = 0.0_wp 10070 DO ib = 1, nbins_aerosol 10071 ! 10072 !-- Diameter in micrometres 10073 mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp 10074 ! 10075 !-- Deposition factor: alveolar (use ra_dry) 10076 df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) + & 10077 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) - & 10078 1.362_wp )**2 ) ) 10079 ! 10080 !-- Lung-deposited surface area LDSA (units mum2/cm3) 10081 temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp * & 10082 aerosol_number(ib)%conc(k,j,i) 10083 ENDDO 10084 ldsa_av(k,j,i) = ldsa_av(k,j,i) + temp_bin 9872 10085 ENDDO 9873 10086 ENDDO 9874 10087 ENDDO 9875 ENDDO 9876 9877 CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4', 'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8', & 9878 'm_bin9', 'm_bin10', 'm_bin11', 'm_bin12' ) 9879 DO i = nxlg, nxrg 9880 DO j = nysg, nyng 9881 DO k = nzb, nzt+1 9882 DO ib = 1, nbins_aerosol 9883 DO ic = ib, nbins_aerosol * ncomponents_mass, nbins_aerosol 9884 mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) + aerosol_mass(ic)%conc(k,j,i) 10088 10089 CASE ( 'N_UFP' ) 10090 DO i = nxlg, nxrg 10091 DO j = nysg, nyng 10092 DO k = nzb, nzt+1 10093 temp_bin = 0.0_wp 10094 DO ib = 1, nbins_aerosol 10095 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp ) THEN 10096 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 10097 ENDIF 10098 ENDDO 10099 nufp_av(k,j,i) = nufp_av(k,j,i) + temp_bin 10100 ENDDO 10101 ENDDO 10102 ENDDO 10103 10104 CASE ( 'Ntot' ) 10105 DO i = nxlg, nxrg 10106 DO j = nysg, nyng 10107 DO k = nzb, nzt+1 10108 DO ib = 1, nbins_aerosol 10109 ntot_av(k,j,i) = ntot_av(k,j,i) + aerosol_number(ib)%conc(k,j,i) 9885 10110 ENDDO 9886 10111 ENDDO 9887 10112 ENDDO 9888 10113 ENDDO 9889 ENDDO 9890 9891 CASE ( 'PM2.5' ) 9892 DO i = nxlg, nxrg 9893 DO j = nysg, nyng 9894 DO k = nzb, nzt+1 9895 temp_bin = 0.0_wp 9896 DO ib = 1, nbins_aerosol 9897 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp ) THEN 9898 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 9899 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 9900 ENDDO 9901 ENDIF 10114 10115 CASE ( 'PM0.1' ) 10116 DO i = nxlg, nxrg 10117 DO j = nysg, nyng 10118 DO k = nzb, nzt+1 10119 temp_bin = 0.0_wp 10120 DO ib = 1, nbins_aerosol 10121 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp ) THEN 10122 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 10123 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10124 ENDDO 10125 ENDIF 10126 ENDDO 10127 pm01_av(k,j,i) = pm01_av(k,j,i) + temp_bin 9902 10128 ENDDO 9903 pm25_av(k,j,i) = pm25_av(k,j,i) + temp_bin9904 10129 ENDDO 9905 10130 ENDDO 9906 ENDDO 9907 9908 CASE ( 'PM10' ) 9909 DO i = nxlg, nxrg 9910 DO j = nysg, nyng 9911 DO k = nzb, nzt+1 9912 temp_bin = 0.0_wp 9913 DO ib = 1, nbins_aerosol 9914 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp ) THEN 9915 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 9916 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 9917 ENDDO 9918 ENDIF 10131 10132 CASE ( 'PM2.5' ) 10133 DO i = nxlg, nxrg 10134 DO j = nysg, nyng 10135 DO k = nzb, nzt+1 10136 temp_bin = 0.0_wp 10137 DO ib = 1, nbins_aerosol 10138 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp ) THEN 10139 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 10140 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10141 ENDDO 10142 ENDIF 10143 ENDDO 10144 pm25_av(k,j,i) = pm25_av(k,j,i) + temp_bin 9919 10145 ENDDO 9920 pm10_av(k,j,i) = pm10_av(k,j,i) + temp_bin9921 10146 ENDDO 9922 10147 ENDDO 9923 ENDDO 9924 9925 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 9926 IF ( is_used( prtcl, TRIM( variable(3:) ) ) ) THEN 9927 found_index = get_index( prtcl, TRIM( variable(3:) ) ) 9928 IF ( TRIM( variable(3:) ) == 'BC' ) to_be_resorted => s_bc_av 9929 IF ( TRIM( variable(3:) ) == 'DU' ) to_be_resorted => s_du_av 9930 IF ( TRIM( variable(3:) ) == 'NH' ) to_be_resorted => s_nh_av 9931 IF ( TRIM( variable(3:) ) == 'NO' ) to_be_resorted => s_no_av 9932 IF ( TRIM( variable(3:) ) == 'OC' ) to_be_resorted => s_oc_av 9933 IF ( TRIM( variable(3:) ) == 'SO4' ) to_be_resorted => s_so4_av 9934 IF ( TRIM( variable(3:) ) == 'SS' ) to_be_resorted => s_ss_av 10148 10149 CASE ( 'PM10' ) 10150 DO i = nxlg, nxrg 10151 DO j = nysg, nyng 10152 DO k = nzb, nzt+1 10153 temp_bin = 0.0_wp 10154 DO ib = 1, nbins_aerosol 10155 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp ) THEN 10156 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 10157 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10158 ENDDO 10159 ENDIF 10160 ENDDO 10161 pm10_av(k,j,i) = pm10_av(k,j,i) + temp_bin 10162 ENDDO 10163 ENDDO 10164 ENDDO 10165 10166 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 10167 IF ( is_used( prtcl, TRIM( variable(9:) ) ) ) THEN ! 9: remove salsa_s_ 10168 found_index = get_index( prtcl, TRIM( variable(9:) ) ) 10169 IF ( TRIM( variable(9:) ) == 'BC' ) to_be_resorted => s_bc_av 10170 IF ( TRIM( variable(9:) ) == 'DU' ) to_be_resorted => s_du_av 10171 IF ( TRIM( variable(9:) ) == 'NH' ) to_be_resorted => s_nh_av 10172 IF ( TRIM( variable(9:) ) == 'NO' ) to_be_resorted => s_no_av 10173 IF ( TRIM( variable(9:) ) == 'OC' ) to_be_resorted => s_oc_av 10174 IF ( TRIM( variable(9:) ) == 'SO4' ) to_be_resorted => s_so4_av 10175 IF ( TRIM( variable(9:) ) == 'SS' ) to_be_resorted => s_ss_av 10176 DO i = nxlg, nxrg 10177 DO j = nysg, nyng 10178 DO k = nzb, nzt+1 10179 DO ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol 10180 to_be_resorted(k,j,i) = to_be_resorted(k,j,i) + & 10181 aerosol_mass(ic)%conc(k,j,i) 10182 ENDDO 10183 ENDDO 10184 ENDDO 10185 ENDDO 10186 ENDIF 10187 10188 CASE ( 's_H2O' ) 10189 found_index = get_index( prtcl,'H2O' ) 10190 to_be_resorted => s_h2o_av 9935 10191 DO i = nxlg, nxrg 9936 10192 DO j = nysg, nyng 9937 10193 DO k = nzb, nzt+1 9938 10194 DO ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol 9939 to_be_resorted(k,j,i) = to_be_resorted(k,j,i) + & 9940 aerosol_mass(ic)%conc(k,j,i) 10195 s_h2o_av(k,j,i) = s_h2o_av(k,j,i) + aerosol_mass(ic)%conc(k,j,i) 9941 10196 ENDDO 9942 10197 ENDDO 9943 10198 ENDDO 9944 10199 ENDDO 9945 ENDIF 9946 9947 CASE ( 's_H2O' ) 9948 found_index = get_index( prtcl,'H2O' ) 9949 to_be_resorted => s_h2o_av 10200 10201 CASE DEFAULT 10202 CONTINUE 10203 10204 END SELECT 10205 10206 ENDIF 10207 10208 ELSEIF ( mode == 'average' ) THEN 10209 10210 IF ( variable(7:11) == 'N_bin' ) THEN 10211 READ( variable(12:),* ) char_to_int 10212 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 10213 ib = char_to_int 9950 10214 DO i = nxlg, nxrg 9951 10215 DO j = nysg, nyng 9952 10216 DO k = nzb, nzt+1 9953 DO ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol 9954 to_be_resorted(k,j,i) = to_be_resorted(k,j,i) + & 9955 aerosol_mass(ic)%conc(k,j,i) 9956 ENDDO 10217 nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp ) 9957 10218 ENDDO 9958 10219 ENDDO 9959 10220 ENDDO 9960 9961 CASE DEFAULT 9962 CONTINUE 9963 9964 END SELECT 9965 9966 ELSEIF ( mode == 'average' ) THEN 9967 9968 SELECT CASE ( TRIM( variable ) ) 9969 9970 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 9971 IF ( TRIM( variable(3:) ) == 'H2SO4' ) THEN 9972 found_index = 1 9973 to_be_resorted => g_h2so4_av 9974 ELSEIF ( TRIM( variable(3:) ) == 'HNO3' ) THEN 9975 found_index = 2 9976 to_be_resorted => g_hno3_av 9977 ELSEIF ( TRIM( variable(3:) ) == 'NH3' ) THEN 9978 found_index = 3 9979 to_be_resorted => g_nh3_av 9980 ELSEIF ( TRIM( variable(3:) ) == 'OCNV' ) THEN 9981 found_index = 4 9982 to_be_resorted => g_ocnv_av 9983 ELSEIF ( TRIM( variable(3:) ) == 'OCSV' ) THEN 9984 found_index = 5 9985 to_be_resorted => g_ocsv_av 9986 ENDIF 10221 ENDIF 10222 10223 ELSEIF ( variable(7:11) == 'm_bin' ) THEN 10224 READ( variable(12:),* ) char_to_int 10225 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 10226 ib = char_to_int 9987 10227 DO i = nxlg, nxrg 9988 10228 DO j = nysg, nyng 9989 10229 DO k = nzb, nzt+1 9990 to_be_resorted(k,j,i) = to_be_resorted(k,j,i) / & 9991 REAL( average_count_3d, KIND=wp ) 10230 mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp) 9992 10231 ENDDO 9993 10232 ENDDO 9994 10233 ENDDO 9995 9996 CASE ( 'LDSA' ) 9997 DO i = nxlg, nxrg 9998 DO j = nysg, nyng 9999 DO k = nzb, nzt+1 10000 ldsa_av(k,j,i) = ldsa_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 10001 ENDDO 10002 ENDDO 10003 ENDDO 10004 10005 CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4', 'N_bin5', 'N_bin6', 'N_bin7', 'N_bin8', & 10006 'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12' ) 10007 DO i = nxlg, nxrg 10008 DO j = nysg, nyng 10009 DO k = nzb, nzt+1 10010 DO ib = 1, nbins_aerosol 10011 nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp ) 10012 ENDDO 10013 ENDDO 10014 ENDDO 10015 ENDDO 10016 10017 CASE ( 'Ntot' ) 10018 DO i = nxlg, nxrg 10019 DO j = nysg, nyng 10020 DO k = nzb, nzt+1 10021 ntot_av(k,j,i) = ntot_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 10022 ENDDO 10023 ENDDO 10024 ENDDO 10025 10026 CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4', 'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8', & 10027 'm_bin9', 'm_bin10', 'm_bin11', 'm_bin12' ) 10028 DO i = nxlg, nxrg 10029 DO j = nysg, nyng 10030 DO k = nzb, nzt+1 10031 DO ib = 1, nbins_aerosol 10032 DO ic = ib, nbins_aerosol * ncomponents_mass, nbins_aerosol 10033 mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) / & 10034 REAL( average_count_3d, KIND=wp) 10035 ENDDO 10036 ENDDO 10037 ENDDO 10038 ENDDO 10039 ENDDO 10040 10041 CASE ( 'PM2.5' ) 10042 DO i = nxlg, nxrg 10043 DO j = nysg, nyng 10044 DO k = nzb, nzt+1 10045 pm25_av(k,j,i) = pm25_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 10046 ENDDO 10047 ENDDO 10048 ENDDO 10049 10050 CASE ( 'PM10' ) 10051 DO i = nxlg, nxrg 10052 DO j = nysg, nyng 10053 DO k = nzb, nzt+1 10054 pm10_av(k,j,i) = pm10_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 10055 ENDDO 10056 ENDDO 10057 ENDDO 10058 10059 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 10060 IF ( is_used( prtcl, TRIM( variable(3:) ) ) ) THEN 10061 IF ( TRIM( variable(3:) ) == 'BC' ) to_be_resorted => s_bc_av 10062 IF ( TRIM( variable(3:) ) == 'DU' ) to_be_resorted => s_du_av 10063 IF ( TRIM( variable(3:) ) == 'NH' ) to_be_resorted => s_nh_av 10064 IF ( TRIM( variable(3:) ) == 'NO' ) to_be_resorted => s_no_av 10065 IF ( TRIM( variable(3:) ) == 'OC' ) to_be_resorted => s_oc_av 10066 IF ( TRIM( variable(3:) ) == 'SO4' ) to_be_resorted => s_so4_av 10067 IF ( TRIM( variable(3:) ) == 'SS' ) to_be_resorted => s_ss_av 10234 ENDIF 10235 ELSE 10236 10237 SELECT CASE ( TRIM( variable(7:) ) ) 10238 10239 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 10240 IF ( TRIM( variable(9:) ) == 'H2SO4' ) THEN ! 9: remove salsa_g_ from beginning 10241 found_index = 1 10242 to_be_resorted => g_h2so4_av 10243 ELSEIF ( TRIM( variable(9:) ) == 'HNO3' ) THEN 10244 found_index = 2 10245 to_be_resorted => g_hno3_av 10246 ELSEIF ( TRIM( variable(9:) ) == 'NH3' ) THEN 10247 found_index = 3 10248 to_be_resorted => g_nh3_av 10249 ELSEIF ( TRIM( variable(9:) ) == 'OCNV' ) THEN 10250 found_index = 4 10251 to_be_resorted => g_ocnv_av 10252 ELSEIF ( TRIM( variable(9:) ) == 'OCSV' ) THEN 10253 found_index = 5 10254 to_be_resorted => g_ocsv_av 10255 ENDIF 10068 10256 DO i = nxlg, nxrg 10069 10257 DO j = nysg, nyng … … 10074 10262 ENDDO 10075 10263 ENDDO 10076 ENDIF 10077 10078 CASE ( 's_H2O' ) 10079 to_be_resorted => s_h2o_av 10080 DO i = nxlg, nxrg 10081 DO j = nysg, nyng 10082 DO k = nzb, nzt+1 10083 to_be_resorted(k,j,i) = to_be_resorted(k,j,i) / & 10084 REAL( average_count_3d, KIND=wp ) 10264 10265 CASE ( 'LDSA' ) 10266 DO i = nxlg, nxrg 10267 DO j = nysg, nyng 10268 DO k = nzb, nzt+1 10269 ldsa_av(k,j,i) = ldsa_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 10270 ENDDO 10085 10271 ENDDO 10086 10272 ENDDO 10087 ENDDO 10088 10089 END SELECT 10090 10273 10274 CASE ( 'N_UFP' ) 10275 DO i = nxlg, nxrg 10276 DO j = nysg, nyng 10277 DO k = nzb, nzt+1 10278 nufp_av(k,j,i) = nufp_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 10279 ENDDO 10280 ENDDO 10281 ENDDO 10282 10283 CASE ( 'Ntot' ) 10284 DO i = nxlg, nxrg 10285 DO j = nysg, nyng 10286 DO k = nzb, nzt+1 10287 ntot_av(k,j,i) = ntot_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 10288 ENDDO 10289 ENDDO 10290 ENDDO 10291 10292 10293 CASE ( 'PM0.1' ) 10294 DO i = nxlg, nxrg 10295 DO j = nysg, nyng 10296 DO k = nzb, nzt+1 10297 pm01_av(k,j,i) = pm01_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 10298 ENDDO 10299 ENDDO 10300 ENDDO 10301 10302 CASE ( 'PM2.5' ) 10303 DO i = nxlg, nxrg 10304 DO j = nysg, nyng 10305 DO k = nzb, nzt+1 10306 pm25_av(k,j,i) = pm25_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 10307 ENDDO 10308 ENDDO 10309 ENDDO 10310 10311 CASE ( 'PM10' ) 10312 DO i = nxlg, nxrg 10313 DO j = nysg, nyng 10314 DO k = nzb, nzt+1 10315 pm10_av(k,j,i) = pm10_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 10316 ENDDO 10317 ENDDO 10318 ENDDO 10319 10320 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 10321 IF ( is_used( prtcl, TRIM( variable(9:) ) ) ) THEN ! 9: remove salsa_s_ 10322 IF ( TRIM( variable(9:) ) == 'BC' ) to_be_resorted => s_bc_av 10323 IF ( TRIM( variable(9:) ) == 'DU' ) to_be_resorted => s_du_av 10324 IF ( TRIM( variable(9:) ) == 'NH' ) to_be_resorted => s_nh_av 10325 IF ( TRIM( variable(9:) ) == 'NO' ) to_be_resorted => s_no_av 10326 IF ( TRIM( variable(9:) ) == 'OC' ) to_be_resorted => s_oc_av 10327 IF ( TRIM( variable(9:) ) == 'SO4' ) to_be_resorted => s_so4_av 10328 IF ( TRIM( variable(9:) ) == 'SS' ) to_be_resorted => s_ss_av 10329 DO i = nxlg, nxrg 10330 DO j = nysg, nyng 10331 DO k = nzb, nzt+1 10332 to_be_resorted(k,j,i) = to_be_resorted(k,j,i) / & 10333 REAL( average_count_3d, KIND=wp ) 10334 ENDDO 10335 ENDDO 10336 ENDDO 10337 ENDIF 10338 10339 CASE ( 's_H2O' ) 10340 to_be_resorted => s_h2o_av 10341 DO i = nxlg, nxrg 10342 DO j = nysg, nyng 10343 DO k = nzb, nzt+1 10344 to_be_resorted(k,j,i) = to_be_resorted(k,j,i) / & 10345 REAL( average_count_3d, KIND=wp ) 10346 ENDDO 10347 ENDDO 10348 ENDDO 10349 10350 END SELECT 10351 10352 ENDIF 10091 10353 ENDIF 10092 10354 … … 10115 10377 10116 10378 INTEGER(iwp) :: av !< 10379 INTEGER(iwp) :: char_to_int !< for converting character to integer 10117 10380 INTEGER(iwp) :: found_index !< index of a chemical compound 10118 10381 INTEGER(iwp) :: i !< … … 10132 10395 REAL(wp) :: fill_value = -9999.0_wp !< value for the _FillValue attribute 10133 10396 REAL(wp) :: mean_d !< Particle diameter in micrometres 10134 REAL(wp) :: nc !< Particle number concentration in units 1/cm**310135 10397 REAL(wp) :: temp_bin !< temporary array for calculating output variables 10136 10398 … … 10145 10407 temp_bin = 0.0_wp 10146 10408 10147 SELECT CASE ( TRIM( variable( 1:LEN( TRIM( variable ) ) - 3 ) ) ) ! cut out _xy, _xz or _yz 10148 10149 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 10150 vari = TRIM( variable( 3:LEN( TRIM( variable ) ) - 3 ) ) 10151 IF ( av == 0 ) THEN 10152 IF ( vari == 'H2SO4') found_index = 1 10153 IF ( vari == 'HNO3') found_index = 2 10154 IF ( vari == 'NH3') found_index = 3 10155 IF ( vari == 'OCNV') found_index = 4 10156 IF ( vari == 'OCSV') found_index = 5 10157 DO i = nxl, nxr 10158 DO j = nys, nyn 10159 DO k = nzb_do, nzt_do 10160 local_pf(i,j,k) = MERGE( salsa_gas(found_index)%conc(k,j,i), REAL( fill_value, & 10161 KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 10162 ENDDO 10163 ENDDO 10164 ENDDO 10165 ELSE 10166 IF ( vari == 'H2SO4' ) to_be_resorted => g_h2so4_av 10167 IF ( vari == 'HNO3' ) to_be_resorted => g_hno3_av 10168 IF ( vari == 'NH3' ) to_be_resorted => g_nh3_av 10169 IF ( vari == 'OCNV' ) to_be_resorted => g_ocnv_av 10170 IF ( vari == 'OCSV' ) to_be_resorted => g_ocsv_av 10171 DO i = nxl, nxr 10172 DO j = nys, nyn 10173 DO k = nzb_do, nzt_do 10174 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, & 10175 KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 10176 ENDDO 10177 ENDDO 10178 ENDDO 10179 ENDIF 10180 10181 IF ( mode == 'xy' ) grid = 'zu' 10182 10183 CASE ( 'LDSA' ) 10184 IF ( av == 0 ) THEN 10185 DO i = nxl, nxr 10186 DO j = nys, nyn 10187 DO k = nzb_do, nzt_do 10188 temp_bin = 0.0_wp 10189 DO ib = 1, nbins_aerosol 10190 ! 10191 !-- Diameter in micrometres 10192 mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp 10193 ! 10194 !-- Deposition factor: alveolar 10195 df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) + & 10196 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) - & 10197 1.362_wp )**2 ) ) 10198 ! 10199 !-- Number concentration in 1/cm3 10200 nc = 1.0E-6_wp * aerosol_number(ib)%conc(k,j,i) 10201 ! 10202 !-- Lung-deposited surface area LDSA (units mum2/cm3) 10203 temp_bin = temp_bin + pi * mean_d**2 * df * nc 10204 ENDDO 10205 10206 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10207 BTEST( wall_flags_0(k,j,i), 0 ) ) 10208 ENDDO 10209 ENDDO 10210 ENDDO 10211 ELSE 10212 DO i = nxl, nxr 10213 DO j = nys, nyn 10214 DO k = nzb_do, nzt_do 10215 local_pf(i,j,k) = MERGE( ldsa_av(k,j,i), REAL( fill_value, KIND = wp ), & 10216 BTEST( wall_flags_0(k,j,i), 0 ) ) 10217 ENDDO 10218 ENDDO 10219 ENDDO 10220 ENDIF 10221 10222 IF ( mode == 'xy' ) grid = 'zu' 10223 10224 CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4', 'N_bin5', 'N_bin6', 'N_bin7', 'N_bin8', & 10225 'N_bin9', 'N_bin10' , 'N_bin11', 'N_bin12' ) 10226 vari = TRIM( variable( 6:LEN( TRIM( variable ) ) - 3 ) ) 10227 10228 IF ( vari == '1' ) ib = 1 10229 IF ( vari == '2' ) ib = 2 10230 IF ( vari == '3' ) ib = 3 10231 IF ( vari == '4' ) ib = 4 10232 IF ( vari == '5' ) ib = 5 10233 IF ( vari == '6' ) ib = 6 10234 IF ( vari == '7' ) ib = 7 10235 IF ( vari == '8' ) ib = 8 10236 IF ( vari == '9' ) ib = 9 10237 IF ( vari == '10' ) ib = 10 10238 IF ( vari == '11' ) ib = 11 10239 IF ( vari == '12' ) ib = 12 10240 10409 IF ( variable(7:11) == 'N_bin' ) THEN 10410 10411 READ( variable( 12:LEN( TRIM( variable ) ) - 3 ), * ) char_to_int 10412 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 10413 10414 ib = char_to_int 10241 10415 IF ( av == 0 ) THEN 10242 10416 DO i = nxl, nxr … … 10258 10432 ENDDO 10259 10433 ENDIF 10260 10261 10434 IF ( mode == 'xy' ) grid = 'zu' 10262 10263 CASE ( 'Ntot' ) 10264 10265 IF ( av == 0 ) THEN 10266 DO i = nxl, nxr 10267 DO j = nys, nyn 10268 DO k = nzb_do, nzt_do 10269 temp_bin = 0.0_wp 10270 DO ib = 1, nbins_aerosol 10271 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 10272 ENDDO 10273 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10274 BTEST( wall_flags_0(k,j,i), 0 ) ) 10275 ENDDO 10276 ENDDO 10277 ENDDO 10278 ELSE 10279 DO i = nxl, nxr 10280 DO j = nys, nyn 10281 DO k = nzb_do, nzt_do 10282 local_pf(i,j,k) = MERGE( ntot_av(k,j,i), REAL( fill_value, KIND = wp ), & 10283 BTEST( wall_flags_0(k,j,i), 0 ) ) 10284 ENDDO 10285 ENDDO 10286 ENDDO 10287 ENDIF 10288 10289 IF ( mode == 'xy' ) grid = 'zu' 10290 10291 CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4', 'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8', & 10292 'm_bin9', 'm_bin10' , 'm_bin11', 'm_bin12' ) 10293 vari = TRIM( variable( 6:LEN( TRIM( variable ) ) - 3 ) ) 10294 10295 IF ( vari == '1' ) ib = 1 10296 IF ( vari == '2' ) ib = 2 10297 IF ( vari == '3' ) ib = 3 10298 IF ( vari == '4' ) ib = 4 10299 IF ( vari == '5' ) ib = 5 10300 IF ( vari == '6' ) ib = 6 10301 IF ( vari == '7' ) ib = 7 10302 IF ( vari == '8' ) ib = 8 10303 IF ( vari == '9' ) ib = 9 10304 IF ( vari == '10' ) ib = 10 10305 IF ( vari == '11' ) ib = 11 10306 IF ( vari == '12' ) ib = 12 10307 10435 ENDIF 10436 10437 ELSEIF ( variable(7:11) == 'm_bin' ) THEN 10438 10439 READ( variable( 12:LEN( TRIM( variable ) ) - 3 ), * ) char_to_int 10440 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 10441 10442 ib = char_to_int 10308 10443 IF ( av == 0 ) THEN 10309 10444 DO i = nxl, nxr … … 10329 10464 ENDDO 10330 10465 ENDIF 10331 10332 10466 IF ( mode == 'xy' ) grid = 'zu' 10333 10334 CASE ( 'PM2.5' ) 10335 IF ( av == 0 ) THEN 10336 DO i = nxl, nxr 10337 DO j = nys, nyn 10338 DO k = nzb_do, nzt_do 10339 temp_bin = 0.0_wp 10340 DO ib = 1, nbins_aerosol 10341 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp ) THEN 10342 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 10467 ENDIF 10468 10469 ELSE 10470 10471 SELECT CASE ( TRIM( variable( 7:LEN( TRIM( variable ) ) - 3 ) ) ) ! cut out _xy, _xz or _yz 10472 10473 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 10474 vari = TRIM( variable( 9:LEN( TRIM( variable ) ) - 3 ) ) ! 9: remove salsa_g_ 10475 IF ( av == 0 ) THEN 10476 IF ( vari == 'H2SO4') found_index = 1 10477 IF ( vari == 'HNO3') found_index = 2 10478 IF ( vari == 'NH3') found_index = 3 10479 IF ( vari == 'OCNV') found_index = 4 10480 IF ( vari == 'OCSV') found_index = 5 10481 DO i = nxl, nxr 10482 DO j = nys, nyn 10483 DO k = nzb_do, nzt_do 10484 local_pf(i,j,k) = MERGE( salsa_gas(found_index)%conc(k,j,i), & 10485 REAL( fill_value, KIND = wp ), & 10486 BTEST( wall_flags_0(k,j,i), 0 ) ) 10487 ENDDO 10488 ENDDO 10489 ENDDO 10490 ELSE 10491 IF ( vari == 'H2SO4' ) to_be_resorted => g_h2so4_av 10492 IF ( vari == 'HNO3' ) to_be_resorted => g_hno3_av 10493 IF ( vari == 'NH3' ) to_be_resorted => g_nh3_av 10494 IF ( vari == 'OCNV' ) to_be_resorted => g_ocnv_av 10495 IF ( vari == 'OCSV' ) to_be_resorted => g_ocsv_av 10496 DO i = nxl, nxr 10497 DO j = nys, nyn 10498 DO k = nzb_do, nzt_do 10499 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, & 10500 KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 10501 ENDDO 10502 ENDDO 10503 ENDDO 10504 ENDIF 10505 10506 IF ( mode == 'xy' ) grid = 'zu' 10507 10508 CASE ( 'LDSA' ) 10509 IF ( av == 0 ) THEN 10510 DO i = nxl, nxr 10511 DO j = nys, nyn 10512 DO k = nzb_do, nzt_do 10513 temp_bin = 0.0_wp 10514 DO ib = 1, nbins_aerosol 10515 ! 10516 !-- Diameter in micrometres 10517 mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp 10518 ! 10519 !-- Deposition factor: alveolar 10520 df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) + & 10521 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) - & 10522 1.362_wp )**2 ) ) 10523 ! 10524 !-- Lung-deposited surface area LDSA (units mum2/cm3) 10525 temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp * & 10526 aerosol_number(ib)%conc(k,j,i) 10527 ENDDO 10528 10529 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10530 BTEST( wall_flags_0(k,j,i), 0 ) ) 10531 ENDDO 10532 ENDDO 10533 ENDDO 10534 ELSE 10535 DO i = nxl, nxr 10536 DO j = nys, nyn 10537 DO k = nzb_do, nzt_do 10538 local_pf(i,j,k) = MERGE( ldsa_av(k,j,i), REAL( fill_value, KIND = wp ), & 10539 BTEST( wall_flags_0(k,j,i), 0 ) ) 10540 ENDDO 10541 ENDDO 10542 ENDDO 10543 ENDIF 10544 10545 IF ( mode == 'xy' ) grid = 'zu' 10546 10547 CASE ( 'N_UFP' ) 10548 10549 IF ( av == 0 ) THEN 10550 DO i = nxl, nxr 10551 DO j = nys, nyn 10552 DO k = nzb_do, nzt_do 10553 temp_bin = 0.0_wp 10554 DO ib = 1, nbins_aerosol 10555 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp ) THEN 10556 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 10557 ENDIF 10558 ENDDO 10559 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10560 BTEST( wall_flags_0(k,j,i), 0 ) ) 10561 ENDDO 10562 ENDDO 10563 ENDDO 10564 ELSE 10565 DO i = nxl, nxr 10566 DO j = nys, nyn 10567 DO k = nzb_do, nzt_do 10568 local_pf(i,j,k) = MERGE( nufp_av(k,j,i), REAL( fill_value, KIND = wp ), & 10569 BTEST( wall_flags_0(k,j,i), 0 ) ) 10570 ENDDO 10571 ENDDO 10572 ENDDO 10573 ENDIF 10574 10575 IF ( mode == 'xy' ) grid = 'zu' 10576 10577 CASE ( 'Ntot' ) 10578 10579 IF ( av == 0 ) THEN 10580 DO i = nxl, nxr 10581 DO j = nys, nyn 10582 DO k = nzb_do, nzt_do 10583 temp_bin = 0.0_wp 10584 DO ib = 1, nbins_aerosol 10585 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 10586 ENDDO 10587 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10588 BTEST( wall_flags_0(k,j,i), 0 ) ) 10589 ENDDO 10590 ENDDO 10591 ENDDO 10592 ELSE 10593 DO i = nxl, nxr 10594 DO j = nys, nyn 10595 DO k = nzb_do, nzt_do 10596 local_pf(i,j,k) = MERGE( ntot_av(k,j,i), REAL( fill_value, KIND = wp ), & 10597 BTEST( wall_flags_0(k,j,i), 0 ) ) 10598 ENDDO 10599 ENDDO 10600 ENDDO 10601 ENDIF 10602 10603 IF ( mode == 'xy' ) grid = 'zu' 10604 10605 CASE ( 'PM0.1' ) 10606 IF ( av == 0 ) THEN 10607 DO i = nxl, nxr 10608 DO j = nys, nyn 10609 DO k = nzb_do, nzt_do 10610 temp_bin = 0.0_wp 10611 DO ib = 1, nbins_aerosol 10612 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp ) THEN 10613 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 10614 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10615 ENDDO 10616 ENDIF 10617 ENDDO 10618 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10619 BTEST( wall_flags_0(k,j,i), 0 ) ) 10620 ENDDO 10621 ENDDO 10622 ENDDO 10623 ELSE 10624 DO i = nxl, nxr 10625 DO j = nys, nyn 10626 DO k = nzb_do, nzt_do 10627 local_pf(i,j,k) = MERGE( pm01_av(k,j,i), REAL( fill_value, KIND = wp ), & 10628 BTEST( wall_flags_0(k,j,i), 0 ) ) 10629 ENDDO 10630 ENDDO 10631 ENDDO 10632 ENDIF 10633 10634 IF ( mode == 'xy' ) grid = 'zu' 10635 10636 CASE ( 'PM2.5' ) 10637 IF ( av == 0 ) THEN 10638 DO i = nxl, nxr 10639 DO j = nys, nyn 10640 DO k = nzb_do, nzt_do 10641 temp_bin = 0.0_wp 10642 DO ib = 1, nbins_aerosol 10643 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp ) THEN 10644 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 10645 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10646 ENDDO 10647 ENDIF 10648 ENDDO 10649 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10650 BTEST( wall_flags_0(k,j,i), 0 ) ) 10651 ENDDO 10652 ENDDO 10653 ENDDO 10654 ELSE 10655 DO i = nxl, nxr 10656 DO j = nys, nyn 10657 DO k = nzb_do, nzt_do 10658 local_pf(i,j,k) = MERGE( pm25_av(k,j,i), REAL( fill_value, KIND = wp ), & 10659 BTEST( wall_flags_0(k,j,i), 0 ) ) 10660 ENDDO 10661 ENDDO 10662 ENDDO 10663 ENDIF 10664 10665 IF ( mode == 'xy' ) grid = 'zu' 10666 10667 CASE ( 'PM10' ) 10668 IF ( av == 0 ) THEN 10669 DO i = nxl, nxr 10670 DO j = nys, nyn 10671 DO k = nzb_do, nzt_do 10672 temp_bin = 0.0_wp 10673 DO ib = 1, nbins_aerosol 10674 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp ) THEN 10675 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 10676 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10677 ENDDO 10678 ENDIF 10679 ENDDO 10680 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10681 BTEST( wall_flags_0(k,j,i), 0 ) ) 10682 ENDDO 10683 ENDDO 10684 ENDDO 10685 ELSE 10686 DO i = nxl, nxr 10687 DO j = nys, nyn 10688 DO k = nzb_do, nzt_do 10689 local_pf(i,j,k) = MERGE( pm10_av(k,j,i), REAL( fill_value, KIND = wp ), & 10690 BTEST( wall_flags_0(k,j,i), 0 ) ) 10691 ENDDO 10692 ENDDO 10693 ENDDO 10694 ENDIF 10695 10696 IF ( mode == 'xy' ) grid = 'zu' 10697 10698 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 10699 vari = TRIM( variable( 9:LEN( TRIM( variable ) ) - 3 ) ) ! 9: remove salsa_s_ 10700 IF ( is_used( prtcl, vari ) ) THEN 10701 found_index = get_index( prtcl, vari ) 10702 IF ( av == 0 ) THEN 10703 DO i = nxl, nxr 10704 DO j = nys, nyn 10705 DO k = nzb_do, nzt_do 10706 temp_bin = 0.0_wp 10707 DO ic = ( found_index-1 ) * nbins_aerosol+1, found_index * nbins_aerosol 10343 10708 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10344 10709 ENDDO 10345 ENDIF 10710 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10711 BTEST( wall_flags_0(k,j,i), 0 ) ) 10712 ENDDO 10346 10713 ENDDO 10347 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), &10348 BTEST( wall_flags_0(k,j,i), 0 ) )10349 10714 ENDDO 10350 ENDDO 10351 ENDDO 10352 ELSE 10353 DO i = nxl, nxr 10354 DO j = nys, nyn 10355 DO k = nzb_do, nzt_do 10356 local_pf(i,j,k) = MERGE( pm25_av(k,j,i), REAL( fill_value, KIND = wp ), & 10357 BTEST( wall_flags_0(k,j,i), 0 ) ) 10715 ELSE 10716 IF ( vari == 'BC' ) to_be_resorted => s_bc_av 10717 IF ( vari == 'DU' ) to_be_resorted => s_du_av 10718 IF ( vari == 'NH' ) to_be_resorted => s_nh_av 10719 IF ( vari == 'NO' ) to_be_resorted => s_no_av 10720 IF ( vari == 'OC' ) to_be_resorted => s_oc_av 10721 IF ( vari == 'SO4' ) to_be_resorted => s_so4_av 10722 IF ( vari == 'SS' ) to_be_resorted => s_ss_av 10723 DO i = nxl, nxr 10724 DO j = nys, nyn 10725 DO k = nzb_do, nzt_do 10726 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, & 10727 KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 10728 ENDDO 10729 ENDDO 10358 10730 ENDDO 10359 ENDDO 10360 ENDDO 10361 ENDIF 10362 10363 IF ( mode == 'xy' ) grid = 'zu' 10364 10365 CASE ( 'PM10' ) 10366 IF ( av == 0 ) THEN 10367 DO i = nxl, nxr 10368 DO j = nys, nyn 10369 DO k = nzb_do, nzt_do 10370 temp_bin = 0.0_wp 10371 DO ib = 1, nbins_aerosol 10372 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp ) THEN 10373 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 10374 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10375 ENDDO 10376 ENDIF 10377 ENDDO 10378 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10379 BTEST( wall_flags_0(k,j,i), 0 ) ) 10380 ENDDO 10381 ENDDO 10382 ENDDO 10383 ELSE 10384 DO i = nxl, nxr 10385 DO j = nys, nyn 10386 DO k = nzb_do, nzt_do 10387 local_pf(i,j,k) = MERGE( pm10_av(k,j,i), REAL( fill_value, KIND = wp ), & 10388 BTEST( wall_flags_0(k,j,i), 0 ) ) 10389 ENDDO 10390 ENDDO 10391 ENDDO 10392 ENDIF 10393 10394 IF ( mode == 'xy' ) grid = 'zu' 10395 10396 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 10397 vari = TRIM( variable( 3:LEN( TRIM( variable ) ) - 3 ) ) 10398 IF ( is_used( prtcl, vari ) ) THEN 10399 found_index = get_index( prtcl, vari ) 10731 ENDIF 10732 ELSE 10733 local_pf = fill_value 10734 ENDIF 10735 10736 IF ( mode == 'xy' ) grid = 'zu' 10737 10738 CASE ( 's_H2O' ) 10739 found_index = get_index( prtcl, 'H2O' ) 10400 10740 IF ( av == 0 ) THEN 10401 10741 DO i = nxl, nxr … … 10412 10752 ENDDO 10413 10753 ELSE 10414 IF ( vari == 'BC' ) to_be_resorted => s_bc_av 10415 IF ( vari == 'DU' ) to_be_resorted => s_du_av 10416 IF ( vari == 'NH' ) to_be_resorted => s_nh_av 10417 IF ( vari == 'NO' ) to_be_resorted => s_no_av 10418 IF ( vari == 'OC' ) to_be_resorted => s_oc_av 10419 IF ( vari == 'SO4' ) to_be_resorted => s_so4_av 10420 IF ( vari == 'SS' ) to_be_resorted => s_ss_av 10754 to_be_resorted => s_h2o_av 10421 10755 DO i = nxl, nxr 10422 10756 DO j = nys, nyn … … 10428 10762 ENDDO 10429 10763 ENDIF 10430 ELSE 10431 local_pf = fill_value 10432 ENDIF 10433 10434 IF ( mode == 'xy' ) grid = 'zu' 10435 10436 CASE ( 's_H2O' ) 10437 found_index = get_index( prtcl, 'H2O' ) 10438 IF ( av == 0 ) THEN 10439 DO i = nxl, nxr 10440 DO j = nys, nyn 10441 DO k = nzb_do, nzt_do 10442 temp_bin = 0.0_wp 10443 DO ic = ( found_index-1 ) * nbins_aerosol+1, found_index * nbins_aerosol 10444 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10445 ENDDO 10446 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10447 BTEST( wall_flags_0(k,j,i), 0 ) ) 10448 ENDDO 10449 ENDDO 10450 ENDDO 10451 ELSE 10452 to_be_resorted => s_h2o_av 10453 DO i = nxl, nxr 10454 DO j = nys, nyn 10455 DO k = nzb_do, nzt_do 10456 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, & 10457 KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 10458 ENDDO 10459 ENDDO 10460 ENDDO 10461 ENDIF 10462 10463 IF ( mode == 'xy' ) grid = 'zu' 10464 10465 CASE DEFAULT 10466 found = .FALSE. 10467 grid = 'none' 10468 10469 END SELECT 10764 10765 IF ( mode == 'xy' ) grid = 'zu' 10766 10767 CASE DEFAULT 10768 found = .FALSE. 10769 grid = 'none' 10770 10771 END SELECT 10772 10773 ENDIF 10470 10774 10471 10775 END SUBROUTINE salsa_data_output_2d … … 10489 10793 10490 10794 INTEGER(iwp) :: av !< 10795 INTEGER(iwp) :: char_to_int !< for converting character to integer 10491 10796 INTEGER(iwp) :: found_index !< index of a chemical compound 10492 10797 INTEGER(iwp) :: ib !< running index: size bins … … 10505 10810 REAL(wp) :: fill_value = -9999.0_wp !< value for the _FillValue attribute 10506 10811 REAL(wp) :: mean_d !< Particle diameter in micrometres 10507 REAL(wp) :: nc !< Particle number concentration in units 1/cm**310508 10812 REAL(wp) :: temp_bin !< temporary array for calculating output variables 10509 10813 … … 10515 10819 temp_bin = 0.0_wp 10516 10820 10517 SELECT CASE ( TRIM( variable ) ) 10518 10519 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 10520 IF ( av == 0 ) THEN 10521 IF ( TRIM( variable ) == 'g_H2SO4') found_index = 1 10522 IF ( TRIM( variable ) == 'g_HNO3') found_index = 2 10523 IF ( TRIM( variable ) == 'g_NH3') found_index = 3 10524 IF ( TRIM( variable ) == 'g_OCNV') found_index = 4 10525 IF ( TRIM( variable ) == 'g_OCSV') found_index = 5 10526 10527 DO i = nxl, nxr 10528 DO j = nys, nyn 10529 DO k = nzb_do, nzt_do 10530 local_pf(i,j,k) = MERGE( salsa_gas(found_index)%conc(k,j,i), & 10531 REAL( fill_value, KIND = wp ), & 10532 BTEST( wall_flags_0(k,j,i), 0 ) ) 10533 ENDDO 10534 ENDDO 10535 ENDDO 10536 ELSE 10537 IF ( TRIM( variable(3:) ) == 'H2SO4' ) to_be_resorted => g_h2so4_av 10538 IF ( TRIM( variable(3:) ) == 'HNO3' ) to_be_resorted => g_hno3_av 10539 IF ( TRIM( variable(3:) ) == 'NH3' ) to_be_resorted => g_nh3_av 10540 IF ( TRIM( variable(3:) ) == 'OCNV' ) to_be_resorted => g_ocnv_av 10541 IF ( TRIM( variable(3:) ) == 'OCSV' ) to_be_resorted => g_ocsv_av 10542 DO i = nxl, nxr 10543 DO j = nys, nyn 10544 DO k = nzb_do, nzt_do 10545 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, & 10546 KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 10547 ENDDO 10548 ENDDO 10549 ENDDO 10550 ENDIF 10551 10552 CASE ( 'LDSA' ) 10553 IF ( av == 0 ) THEN 10554 DO i = nxl, nxr 10555 DO j = nys, nyn 10556 DO k = nzb_do, nzt_do 10557 temp_bin = 0.0_wp 10558 DO ib = 1, nbins_aerosol 10559 ! 10560 !-- Diameter in micrometres 10561 mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp 10562 ! 10563 !-- Deposition factor: alveolar 10564 df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) + & 10565 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) - & 10566 1.362_wp )**2 ) ) 10567 ! 10568 !-- Number concentration in 1/cm3 10569 nc = 1.0E-6_wp * aerosol_number(ib)%conc(k,j,i) 10570 ! 10571 !-- Lung-deposited surface area LDSA (units mum2/cm3) 10572 temp_bin = temp_bin + pi * mean_d**2 * df * nc 10573 ENDDO 10574 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10575 BTEST( wall_flags_0(k,j,i), 0 ) ) 10576 ENDDO 10577 ENDDO 10578 ENDDO 10579 ELSE 10580 DO i = nxl, nxr 10581 DO j = nys, nyn 10582 DO k = nzb_do, nzt_do 10583 local_pf(i,j,k) = MERGE( ldsa_av(k,j,i), REAL( fill_value, KIND = wp ), & 10584 BTEST( wall_flags_0(k,j,i), 0 ) ) 10585 ENDDO 10586 ENDDO 10587 ENDDO 10588 ENDIF 10589 10590 CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4', 'N_bin5', 'N_bin6', 'N_bin7', 'N_bin8', & 10591 'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12' ) 10592 IF ( TRIM( variable(6:) ) == '1' ) ib = 1 10593 IF ( TRIM( variable(6:) ) == '2' ) ib = 2 10594 IF ( TRIM( variable(6:) ) == '3' ) ib = 3 10595 IF ( TRIM( variable(6:) ) == '4' ) ib = 4 10596 IF ( TRIM( variable(6:) ) == '5' ) ib = 5 10597 IF ( TRIM( variable(6:) ) == '6' ) ib = 6 10598 IF ( TRIM( variable(6:) ) == '7' ) ib = 7 10599 IF ( TRIM( variable(6:) ) == '8' ) ib = 8 10600 IF ( TRIM( variable(6:) ) == '9' ) ib = 9 10601 IF ( TRIM( variable(6:) ) == '10' ) ib = 10 10602 IF ( TRIM( variable(6:) ) == '11' ) ib = 11 10603 IF ( TRIM( variable(6:) ) == '12' ) ib = 12 10604 10821 IF ( variable(7:11) == 'N_bin' ) THEN 10822 READ( variable(12:),* ) char_to_int 10823 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 10824 10825 ib = char_to_int 10605 10826 IF ( av == 0 ) THEN 10606 10827 DO i = nxl, nxr … … 10622 10843 ENDDO 10623 10844 ENDIF 10624 10625 CASE ( 'Ntot' ) 10626 IF ( av == 0 ) THEN 10627 DO i = nxl, nxr 10628 DO j = nys, nyn 10629 DO k = nzb_do, nzt_do 10630 temp_bin = 0.0_wp 10631 DO ib = 1, nbins_aerosol 10632 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 10633 ENDDO 10634 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10635 BTEST( wall_flags_0(k,j,i), 0 ) ) 10636 ENDDO 10637 ENDDO 10638 ENDDO 10639 ELSE 10640 DO i = nxl, nxr 10641 DO j = nys, nyn 10642 DO k = nzb_do, nzt_do 10643 local_pf(i,j,k) = MERGE( ntot_av(k,j,i), REAL( fill_value, KIND = wp ), & 10644 BTEST( wall_flags_0(k,j,i), 0 ) ) 10645 ENDDO 10646 ENDDO 10647 ENDDO 10648 ENDIF 10649 10650 CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4', 'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8', & 10651 'm_bin9', 'm_bin10' , 'm_bin11', 'm_bin12' ) 10652 IF ( TRIM( variable(6:) ) == '1' ) ib = 1 10653 IF ( TRIM( variable(6:) ) == '2' ) ib = 2 10654 IF ( TRIM( variable(6:) ) == '3' ) ib = 3 10655 IF ( TRIM( variable(6:) ) == '4' ) ib = 4 10656 IF ( TRIM( variable(6:) ) == '5' ) ib = 5 10657 IF ( TRIM( variable(6:) ) == '6' ) ib = 6 10658 IF ( TRIM( variable(6:) ) == '7' ) ib = 7 10659 IF ( TRIM( variable(6:) ) == '8' ) ib = 8 10660 IF ( TRIM( variable(6:) ) == '9' ) ib = 9 10661 IF ( TRIM( variable(6:) ) == '10' ) ib = 10 10662 IF ( TRIM( variable(6:) ) == '11' ) ib = 11 10663 IF ( TRIM( variable(6:) ) == '12' ) ib = 12 10664 10845 ENDIF 10846 10847 ELSEIF ( variable(7:11) == 'm_bin' ) THEN 10848 READ( variable(12:),* ) char_to_int 10849 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 10850 10851 ib = char_to_int 10665 10852 IF ( av == 0 ) THEN 10666 10853 DO i = nxl, nxr … … 10686 10873 ENDDO 10687 10874 ENDIF 10688 10689 CASE ( 'PM2.5' ) 10690 IF ( av == 0 ) THEN 10691 DO i = nxl, nxr 10692 DO j = nys, nyn 10693 DO k = nzb_do, nzt_do 10694 temp_bin = 0.0_wp 10695 DO ib = 1, nbins_aerosol 10696 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp ) THEN 10697 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 10875 ENDIF 10876 10877 ELSE 10878 SELECT CASE ( TRIM( variable(7:) ) ) 10879 10880 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 10881 IF ( av == 0 ) THEN 10882 IF ( TRIM( variable(7:) ) == 'g_H2SO4') found_index = 1 10883 IF ( TRIM( variable(7:) ) == 'g_HNO3') found_index = 2 10884 IF ( TRIM( variable(7:) ) == 'g_NH3') found_index = 3 10885 IF ( TRIM( variable(7:) ) == 'g_OCNV') found_index = 4 10886 IF ( TRIM( variable(7:) ) == 'g_OCSV') found_index = 5 10887 10888 DO i = nxl, nxr 10889 DO j = nys, nyn 10890 DO k = nzb_do, nzt_do 10891 local_pf(i,j,k) = MERGE( salsa_gas(found_index)%conc(k,j,i), & 10892 REAL( fill_value, KIND = wp ), & 10893 BTEST( wall_flags_0(k,j,i), 0 ) ) 10894 ENDDO 10895 ENDDO 10896 ENDDO 10897 ELSE 10898 ! 10899 !-- 9: remove salsa_g_ from the beginning 10900 IF ( TRIM( variable(9:) ) == 'H2SO4' ) to_be_resorted => g_h2so4_av 10901 IF ( TRIM( variable(9:) ) == 'HNO3' ) to_be_resorted => g_hno3_av 10902 IF ( TRIM( variable(9:) ) == 'NH3' ) to_be_resorted => g_nh3_av 10903 IF ( TRIM( variable(9:) ) == 'OCNV' ) to_be_resorted => g_ocnv_av 10904 IF ( TRIM( variable(9:) ) == 'OCSV' ) to_be_resorted => g_ocsv_av 10905 DO i = nxl, nxr 10906 DO j = nys, nyn 10907 DO k = nzb_do, nzt_do 10908 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, & 10909 KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 10910 ENDDO 10911 ENDDO 10912 ENDDO 10913 ENDIF 10914 10915 CASE ( 'LDSA' ) 10916 IF ( av == 0 ) THEN 10917 DO i = nxl, nxr 10918 DO j = nys, nyn 10919 DO k = nzb_do, nzt_do 10920 temp_bin = 0.0_wp 10921 DO ib = 1, nbins_aerosol 10922 ! 10923 !-- Diameter in micrometres 10924 mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp 10925 ! 10926 !-- Deposition factor: alveolar 10927 df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) + & 10928 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) - & 10929 1.362_wp )**2 ) ) 10930 ! 10931 !-- Lung-deposited surface area LDSA (units mum2/cm3) 10932 temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp * & 10933 aerosol_number(ib)%conc(k,j,i) 10934 ENDDO 10935 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10936 BTEST( wall_flags_0(k,j,i), 0 ) ) 10937 ENDDO 10938 ENDDO 10939 ENDDO 10940 ELSE 10941 DO i = nxl, nxr 10942 DO j = nys, nyn 10943 DO k = nzb_do, nzt_do 10944 local_pf(i,j,k) = MERGE( ldsa_av(k,j,i), REAL( fill_value, KIND = wp ), & 10945 BTEST( wall_flags_0(k,j,i), 0 ) ) 10946 ENDDO 10947 ENDDO 10948 ENDDO 10949 ENDIF 10950 10951 CASE ( 'N_UFP' ) 10952 IF ( av == 0 ) THEN 10953 DO i = nxl, nxr 10954 DO j = nys, nyn 10955 DO k = nzb_do, nzt_do 10956 temp_bin = 0.0_wp 10957 DO ib = 1, nbins_aerosol 10958 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp ) THEN 10959 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 10960 ENDIF 10961 ENDDO 10962 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10963 BTEST( wall_flags_0(k,j,i), 0 ) ) 10964 ENDDO 10965 ENDDO 10966 ENDDO 10967 ELSE 10968 DO i = nxl, nxr 10969 DO j = nys, nyn 10970 DO k = nzb_do, nzt_do 10971 local_pf(i,j,k) = MERGE( nufp_av(k,j,i), REAL( fill_value, KIND = wp ), & 10972 BTEST( wall_flags_0(k,j,i), 0 ) ) 10973 ENDDO 10974 ENDDO 10975 ENDDO 10976 ENDIF 10977 10978 CASE ( 'Ntot' ) 10979 IF ( av == 0 ) THEN 10980 DO i = nxl, nxr 10981 DO j = nys, nyn 10982 DO k = nzb_do, nzt_do 10983 temp_bin = 0.0_wp 10984 DO ib = 1, nbins_aerosol 10985 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 10986 ENDDO 10987 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10988 BTEST( wall_flags_0(k,j,i), 0 ) ) 10989 ENDDO 10990 ENDDO 10991 ENDDO 10992 ELSE 10993 DO i = nxl, nxr 10994 DO j = nys, nyn 10995 DO k = nzb_do, nzt_do 10996 local_pf(i,j,k) = MERGE( ntot_av(k,j,i), REAL( fill_value, KIND = wp ), & 10997 BTEST( wall_flags_0(k,j,i), 0 ) ) 10998 ENDDO 10999 ENDDO 11000 ENDDO 11001 ENDIF 11002 11003 CASE ( 'PM0.1' ) 11004 IF ( av == 0 ) THEN 11005 DO i = nxl, nxr 11006 DO j = nys, nyn 11007 DO k = nzb_do, nzt_do 11008 temp_bin = 0.0_wp 11009 DO ib = 1, nbins_aerosol 11010 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp ) THEN 11011 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 11012 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 11013 ENDDO 11014 ENDIF 11015 ENDDO 11016 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 11017 BTEST( wall_flags_0(k,j,i), 0 ) ) 11018 ENDDO 11019 ENDDO 11020 ENDDO 11021 ELSE 11022 DO i = nxl, nxr 11023 DO j = nys, nyn 11024 DO k = nzb_do, nzt_do 11025 local_pf(i,j,k) = MERGE( pm01_av(k,j,i), REAL( fill_value, KIND = wp ), & 11026 BTEST( wall_flags_0(k,j,i), 0 ) ) 11027 ENDDO 11028 ENDDO 11029 ENDDO 11030 ENDIF 11031 11032 CASE ( 'PM2.5' ) 11033 IF ( av == 0 ) THEN 11034 DO i = nxl, nxr 11035 DO j = nys, nyn 11036 DO k = nzb_do, nzt_do 11037 temp_bin = 0.0_wp 11038 DO ib = 1, nbins_aerosol 11039 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp ) THEN 11040 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 11041 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 11042 ENDDO 11043 ENDIF 11044 ENDDO 11045 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 11046 BTEST( wall_flags_0(k,j,i), 0 ) ) 11047 ENDDO 11048 ENDDO 11049 ENDDO 11050 ELSE 11051 DO i = nxl, nxr 11052 DO j = nys, nyn 11053 DO k = nzb_do, nzt_do 11054 local_pf(i,j,k) = MERGE( pm25_av(k,j,i), REAL( fill_value, KIND = wp ), & 11055 BTEST( wall_flags_0(k,j,i), 0 ) ) 11056 ENDDO 11057 ENDDO 11058 ENDDO 11059 ENDIF 11060 11061 CASE ( 'PM10' ) 11062 IF ( av == 0 ) THEN 11063 DO i = nxl, nxr 11064 DO j = nys, nyn 11065 DO k = nzb_do, nzt_do 11066 temp_bin = 0.0_wp 11067 DO ib = 1, nbins_aerosol 11068 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp ) THEN 11069 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 11070 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 11071 ENDDO 11072 ENDIF 11073 ENDDO 11074 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 11075 BTEST( wall_flags_0(k,j,i), 0 ) ) 11076 ENDDO 11077 ENDDO 11078 ENDDO 11079 ELSE 11080 DO i = nxl, nxr 11081 DO j = nys, nyn 11082 DO k = nzb_do, nzt_do 11083 local_pf(i,j,k) = MERGE( pm10_av(k,j,i), REAL( fill_value, KIND = wp ), & 11084 BTEST( wall_flags_0(k,j,i), 0 ) ) 11085 ENDDO 11086 ENDDO 11087 ENDDO 11088 ENDIF 11089 11090 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 11091 IF ( is_used( prtcl, TRIM( variable(9:) ) ) ) THEN ! 9: remove salsa_s_ 11092 found_index = get_index( prtcl, TRIM( variable(9:) ) ) 11093 IF ( av == 0 ) THEN 11094 DO i = nxl, nxr 11095 DO j = nys, nyn 11096 DO k = nzb_do, nzt_do 11097 temp_bin = 0.0_wp 11098 DO ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol 10698 11099 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10699 11100 ENDDO 10700 ENDIF 11101 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 11102 BTEST( wall_flags_0(k,j,i), 0 ) ) 11103 ENDDO 10701 11104 ENDDO 10702 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), &10703 BTEST( wall_flags_0(k,j,i), 0 ) )10704 11105 ENDDO 10705 ENDDO 10706 ENDDO 10707 ELSE 10708 DO i = nxl, nxr 10709 DO j = nys, nyn 10710 DO k = nzb_do, nzt_do 10711 local_pf(i,j,k) = MERGE( pm25_av(k,j,i), REAL( fill_value, KIND = wp ), & 10712 BTEST( wall_flags_0(k,j,i), 0 ) ) 11106 ELSE 11107 ! 11108 !-- 9: remove salsa_s_ from the beginning 11109 IF ( TRIM( variable(9:) ) == 'BC' ) to_be_resorted => s_bc_av 11110 IF ( TRIM( variable(9:) ) == 'DU' ) to_be_resorted => s_du_av 11111 IF ( TRIM( variable(9:) ) == 'NH' ) to_be_resorted => s_nh_av 11112 IF ( TRIM( variable(9:) ) == 'NO' ) to_be_resorted => s_no_av 11113 IF ( TRIM( variable(9:) ) == 'OC' ) to_be_resorted => s_oc_av 11114 IF ( TRIM( variable(9:) ) == 'SO4' ) to_be_resorted => s_so4_av 11115 IF ( TRIM( variable(9:) ) == 'SS' ) to_be_resorted => s_ss_av 11116 DO i = nxl, nxr 11117 DO j = nys, nyn 11118 DO k = nzb_do, nzt_do 11119 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, & 11120 KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 11121 ENDDO 11122 ENDDO 10713 11123 ENDDO 10714 ENDDO 10715 ENDDO 10716 ENDIF 10717 10718 CASE ( 'PM10' ) 10719 IF ( av == 0 ) THEN 10720 DO i = nxl, nxr 10721 DO j = nys, nyn 10722 DO k = nzb_do, nzt_do 10723 temp_bin = 0.0_wp 10724 DO ib = 1, nbins_aerosol 10725 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp ) THEN 10726 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 10727 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10728 ENDDO 10729 ENDIF 10730 ENDDO 10731 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10732 BTEST( wall_flags_0(k,j,i), 0 ) ) 10733 ENDDO 10734 ENDDO 10735 ENDDO 10736 ELSE 10737 DO i = nxl, nxr 10738 DO j = nys, nyn 10739 DO k = nzb_do, nzt_do 10740 local_pf(i,j,k) = MERGE( pm10_av(k,j,i), REAL( fill_value, KIND = wp ), & 10741 BTEST( wall_flags_0(k,j,i), 0 ) ) 10742 ENDDO 10743 ENDDO 10744 ENDDO 10745 ENDIF 10746 10747 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 10748 IF ( is_used( prtcl, TRIM( variable(3:) ) ) ) THEN 10749 found_index = get_index( prtcl, TRIM( variable(3:) ) ) 11124 ENDIF 11125 ENDIF 11126 11127 CASE ( 's_H2O' ) 11128 found_index = get_index( prtcl, 'H2O' ) 10750 11129 IF ( av == 0 ) THEN 10751 11130 DO i = nxl, nxr … … 10762 11141 ENDDO 10763 11142 ELSE 10764 IF ( TRIM( variable(3:) ) == 'BC' ) to_be_resorted => s_bc_av 10765 IF ( TRIM( variable(3:) ) == 'DU' ) to_be_resorted => s_du_av 10766 IF ( TRIM( variable(3:) ) == 'NH' ) to_be_resorted => s_nh_av 10767 IF ( TRIM( variable(3:) ) == 'NO' ) to_be_resorted => s_no_av 10768 IF ( TRIM( variable(3:) ) == 'OC' ) to_be_resorted => s_oc_av 10769 IF ( TRIM( variable(3:) ) == 'SO4' ) to_be_resorted => s_so4_av 10770 IF ( TRIM( variable(3:) ) == 'SS' ) to_be_resorted => s_ss_av 11143 to_be_resorted => s_h2o_av 10771 11144 DO i = nxl, nxr 10772 11145 DO j = nys, nyn … … 10778 11151 ENDDO 10779 11152 ENDIF 10780 ENDIF 10781 10782 CASE ( 's_H2O' ) 10783 found_index = get_index( prtcl, 'H2O' ) 10784 IF ( av == 0 ) THEN 10785 DO i = nxl, nxr 10786 DO j = nys, nyn 10787 DO k = nzb_do, nzt_do 10788 temp_bin = 0.0_wp 10789 DO ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol 10790 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 10791 ENDDO 10792 local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ), & 10793 BTEST( wall_flags_0(k,j,i), 0 ) ) 10794 ENDDO 10795 ENDDO 10796 ENDDO 10797 ELSE 10798 to_be_resorted => s_h2o_av 10799 DO i = nxl, nxr 10800 DO j = nys, nyn 10801 DO k = nzb_do, nzt_do 10802 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, & 10803 KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 10804 ENDDO 10805 ENDDO 10806 ENDDO 10807 ENDIF 10808 10809 CASE DEFAULT 10810 found = .FALSE. 10811 10812 END SELECT 11153 11154 CASE DEFAULT 11155 found = .FALSE. 11156 11157 END SELECT 11158 ENDIF 10813 11159 10814 11160 END SUBROUTINE salsa_data_output_3d … … 10837 11183 CHARACTER(LEN=7) :: vari !< trimmed format of variable 10838 11184 10839 INTEGER(iwp) :: av !< 10840 INTEGER(iwp) :: found_index !< index of a chemical compound 10841 INTEGER(iwp) :: ib !< loop index for aerosol size number bins 10842 INTEGER(iwp) :: ic !< loop index for chemical components 10843 INTEGER(iwp) :: i !< loop index in x-direction 10844 INTEGER(iwp) :: j !< loop index in y-direction 10845 INTEGER(iwp) :: k !< loop index in z-direction 11185 INTEGER(iwp) :: av !< 11186 INTEGER(iwp) :: char_to_int !< for converting character to integer 11187 INTEGER(iwp) :: found_index !< index of a chemical compound 11188 INTEGER(iwp) :: ib !< loop index for aerosol size number bins 11189 INTEGER(iwp) :: ic !< loop index for chemical components 11190 INTEGER(iwp) :: i !< loop index in x-direction 11191 INTEGER(iwp) :: j !< loop index in y-direction 11192 INTEGER(iwp) :: k !< loop index in z-direction 10846 11193 INTEGER(iwp) :: mid !< masked output running index 10847 INTEGER(iwp) :: topo_top_ind 11194 INTEGER(iwp) :: topo_top_ind !< k index of highest horizontal surface 10848 11195 10849 11196 LOGICAL :: found !< 10850 11197 LOGICAL :: resorted !< 10851 11198 10852 REAL(wp) :: df !< For calculating LDSA: fraction of particles 10853 !< depositing in the alveolar (or tracheobronchial) 10854 !< region of the lung. Depends on the particle size 11199 REAL(wp) :: df !< For calculating LDSA: fraction of particles depositing in the alveolar 11200 !< (or tracheobronchial) region of the lung. Depends on the particle size 10855 11201 REAL(wp) :: mean_d !< Particle diameter in micrometres 10856 REAL(wp) :: nc !< Particle number concentration in units 1/cm**310857 11202 REAL(wp) :: temp_bin !< temporary array for calculating output variables 10858 11203 10859 11204 REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: local_pf !< 10860 11205 11206 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), TARGET :: temp_array !< temporary array 11207 10861 11208 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< pointer 10862 11209 10863 found = .TRUE. 10864 resorted = .FALSE. 10865 grid = 's' 10866 temp_bin = 0.0_wp 10867 10868 SELECT CASE ( TRIM( variable ) ) 10869 10870 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 10871 vari = TRIM( variable ) 10872 IF ( av == 0 ) THEN 10873 IF ( vari == 'g_H2SO4') to_be_resorted => salsa_gas(1)%conc 10874 IF ( vari == 'g_HNO3') to_be_resorted => salsa_gas(2)%conc 10875 IF ( vari == 'g_NH3') to_be_resorted => salsa_gas(3)%conc 10876 IF ( vari == 'g_OCNV') to_be_resorted => salsa_gas(4)%conc 10877 IF ( vari == 'g_OCSV') to_be_resorted => salsa_gas(5)%conc 10878 ELSE 10879 IF ( vari == 'g_H2SO4') to_be_resorted => g_h2so4_av 10880 IF ( vari == 'g_HNO3') to_be_resorted => g_hno3_av 10881 IF ( vari == 'g_NH3') to_be_resorted => g_nh3_av 10882 IF ( vari == 'g_OCNV') to_be_resorted => g_ocnv_av 10883 IF ( vari == 'g_OCSV') to_be_resorted => g_ocsv_av 10884 ENDIF 10885 10886 CASE ( 'LDSA' ) 10887 IF ( av == 0 ) THEN 10888 DO i = nxl, nxr 10889 DO j = nys, nyn 10890 DO k = nzb, nz_do3d 10891 temp_bin = 0.0_wp 10892 DO ib = 1, nbins_aerosol 10893 ! 10894 !-- Diameter in micrometres 10895 mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp 10896 ! 10897 !-- Deposition factor: alveolar 10898 df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) + & 10899 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) - & 10900 1.362_wp )**2 ) ) 10901 ! 10902 !-- Number concentration in 1/cm3 10903 nc = 1.0E-6_wp * aerosol_number(ib)%conc(k,j,i) 10904 ! 10905 !-- Lung-deposited surface area LDSA (units mum2/cm3) 10906 temp_bin = temp_bin + pi * mean_d**2 * df * nc 10907 ENDDO 10908 tend(k,j,i) = temp_bin 10909 ENDDO 10910 ENDDO 10911 ENDDO 10912 IF ( .NOT. mask_surface(mid) ) THEN 10913 DO i = 1, mask_size_l(mid,1) 10914 DO j = 1, mask_size_l(mid,2) 10915 DO k = 1, mask_size_l(mid,3) 10916 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 10917 ENDDO 10918 ENDDO 10919 ENDDO 10920 ELSE 10921 DO i = 1, mask_size_l(mid,1) 10922 DO j = 1, mask_size_l(mid,2) 10923 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 10924 grid ) 10925 DO k = 1, mask_size_l(mid,3) 10926 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 10927 mask_j(mid,j), mask_i(mid,i) ) 10928 ENDDO 10929 ENDDO 10930 ENDDO 10931 ENDIF 10932 resorted = .TRUE. 10933 ELSE 10934 to_be_resorted => ldsa_av 10935 ENDIF 10936 10937 CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4', 'N_bin5', 'N_bin6', 'N_bin7', 'N_bin8', & 10938 'N_bin9', 'N_bin10' , 'N_bin11', 'N_bin12' ) 10939 IF ( TRIM( variable(6:) ) == '1' ) ib = 1 10940 IF ( TRIM( variable(6:) ) == '2' ) ib = 2 10941 IF ( TRIM( variable(6:) ) == '3' ) ib = 3 10942 IF ( TRIM( variable(6:) ) == '4' ) ib = 4 10943 IF ( TRIM( variable(6:) ) == '5' ) ib = 5 10944 IF ( TRIM( variable(6:) ) == '6' ) ib = 6 10945 IF ( TRIM( variable(6:) ) == '7' ) ib = 7 10946 IF ( TRIM( variable(6:) ) == '8' ) ib = 8 10947 IF ( TRIM( variable(6:) ) == '9' ) ib = 9 10948 IF ( TRIM( variable(6:) ) == '10' ) ib = 10 10949 IF ( TRIM( variable(6:) ) == '11' ) ib = 11 10950 IF ( TRIM( variable(6:) ) == '12' ) ib = 12 10951 11210 found = .TRUE. 11211 resorted = .FALSE. 11212 grid = 's' 11213 temp_array = 0.0_wp 11214 temp_bin = 0.0_wp 11215 11216 IF ( variable(7:11) == 'N_bin' ) THEN 11217 READ( variable(12:),* ) char_to_int 11218 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 11219 ib = char_to_int 10952 11220 IF ( av == 0 ) THEN 10953 11221 IF ( .NOT. mask_surface(mid) ) THEN … … 10975 11243 resorted = .TRUE. 10976 11244 ELSE 10977 to_be_resorted => nbins_av(:,:,:,ib) 11245 temp_array = nbins_av(:,:,:,ib) 11246 to_be_resorted => temp_array 10978 11247 ENDIF 10979 10980 CASE ( 'Ntot' ) 10981 IF ( av == 0 ) THEN 10982 DO i = nxl, nxr 10983 DO j = nys, nyn 10984 DO k = nzb, nz_do3d 10985 temp_bin = 0.0_wp 10986 DO ib = 1, nbins_aerosol 10987 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 10988 ENDDO 10989 tend(k,j,i) = temp_bin 10990 ENDDO 10991 ENDDO 10992 ENDDO 10993 IF ( .NOT. mask_surface(mid) ) THEN 10994 DO i = 1, mask_size_l(mid,1) 10995 DO j = 1, mask_size_l(mid,2) 10996 DO k = 1, mask_size_l(mid,3) 10997 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 10998 ENDDO 10999 ENDDO 11000 ENDDO 11001 ELSE 11002 DO i = 1, mask_size_l(mid,1) 11003 DO j = 1, mask_size_l(mid,2) 11004 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11005 grid ) 11006 DO k = 1, mask_size_l(mid,3) 11007 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11008 mask_j(mid,j), mask_i(mid,i) ) 11009 ENDDO 11010 ENDDO 11011 ENDDO 11012 ENDIF 11013 resorted = .TRUE. 11014 ELSE 11015 to_be_resorted => ntot_av 11016 ENDIF 11017 11018 CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4', 'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8', & 11019 'm_bin9', 'm_bin10', 'm_bin11', 'm_bin12' ) 11020 IF ( TRIM( variable(6:) ) == '1' ) ib = 1 11021 IF ( TRIM( variable(6:) ) == '2' ) ib = 2 11022 IF ( TRIM( variable(6:) ) == '3' ) ib = 3 11023 IF ( TRIM( variable(6:) ) == '4' ) ib = 4 11024 IF ( TRIM( variable(6:) ) == '5' ) ib = 5 11025 IF ( TRIM( variable(6:) ) == '6' ) ib = 6 11026 IF ( TRIM( variable(6:) ) == '7' ) ib = 7 11027 IF ( TRIM( variable(6:) ) == '8' ) ib = 8 11028 IF ( TRIM( variable(6:) ) == '9' ) ib = 9 11029 IF ( TRIM( variable(6:) ) == '10' ) ib = 10 11030 IF ( TRIM( variable(6:) ) == '11' ) ib = 11 11031 IF ( TRIM( variable(6:) ) == '12' ) ib = 12 11032 11248 ENDIF 11249 11250 ELSEIF ( variable(7:11) == 'm_bin' ) THEN 11251 11252 READ( variable(12:),* ) char_to_int 11253 IF ( char_to_int >= 1 .AND. char_to_int <= SUM( nbin ) ) THEN 11254 11255 ib = char_to_int 11033 11256 IF ( av == 0 ) THEN 11034 11257 DO i = nxl, nxr … … 11065 11288 resorted = .TRUE. 11066 11289 ELSE 11067 to_be_resorted => mbins_av(:,:,:,ib) 11290 temp_array = mbins_av(:,:,:,ib) 11291 to_be_resorted => temp_array 11068 11292 ENDIF 11069 11070 CASE ( 'PM2.5' ) 11071 IF ( av == 0 ) THEN 11072 DO i = nxl, nxr 11073 DO j = nys, nyn 11074 DO k = nzb, nz_do3d 11075 temp_bin = 0.0_wp 11076 DO ib = 1, nbins_aerosol 11077 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp ) THEN 11078 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 11079 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 11080 ENDDO 11081 ENDIF 11082 ENDDO 11083 tend(k,j,i) = temp_bin 11084 ENDDO 11085 ENDDO 11086 ENDDO 11087 IF ( .NOT. mask_surface(mid) ) THEN 11088 DO i = 1, mask_size_l(mid,1) 11089 DO j = 1, mask_size_l(mid,2) 11090 DO k = 1, mask_size_l(mid,3) 11091 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11293 ENDIF 11294 11295 ELSE 11296 SELECT CASE ( TRIM( variable(7:) ) ) 11297 11298 CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' ) 11299 vari = TRIM( variable(7:) ) 11300 IF ( av == 0 ) THEN 11301 IF ( vari == 'g_H2SO4') to_be_resorted => salsa_gas(1)%conc 11302 IF ( vari == 'g_HNO3') to_be_resorted => salsa_gas(2)%conc 11303 IF ( vari == 'g_NH3') to_be_resorted => salsa_gas(3)%conc 11304 IF ( vari == 'g_OCNV') to_be_resorted => salsa_gas(4)%conc 11305 IF ( vari == 'g_OCSV') to_be_resorted => salsa_gas(5)%conc 11306 ELSE 11307 IF ( vari == 'g_H2SO4') to_be_resorted => g_h2so4_av 11308 IF ( vari == 'g_HNO3') to_be_resorted => g_hno3_av 11309 IF ( vari == 'g_NH3') to_be_resorted => g_nh3_av 11310 IF ( vari == 'g_OCNV') to_be_resorted => g_ocnv_av 11311 IF ( vari == 'g_OCSV') to_be_resorted => g_ocsv_av 11312 ENDIF 11313 11314 CASE ( 'LDSA' ) 11315 IF ( av == 0 ) THEN 11316 DO i = nxl, nxr 11317 DO j = nys, nyn 11318 DO k = nzb, nz_do3d 11319 temp_bin = 0.0_wp 11320 DO ib = 1, nbins_aerosol 11321 ! 11322 !-- Diameter in micrometres 11323 mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp 11324 ! 11325 !-- Deposition factor: alveolar 11326 df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) + & 11327 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) - & 11328 1.362_wp )**2 ) ) 11329 ! 11330 !-- Lung-deposited surface area LDSA (units mum2/cm3) 11331 temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp * & 11332 aerosol_number(ib)%conc(k,j,i) 11333 ENDDO 11334 tend(k,j,i) = temp_bin 11092 11335 ENDDO 11093 11336 ENDDO 11094 11337 ENDDO 11095 ELSE 11096 DO i = 1, mask_size_l(mid,1) 11097 DO j = 1, mask_size_l(mid,2) 11098 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11099 grid ) 11100 DO k = 1, mask_size_l(mid,3) 11101 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11102 mask_j(mid,j), mask_i(mid,i) ) 11338 IF ( .NOT. mask_surface(mid) ) THEN 11339 DO i = 1, mask_size_l(mid,1) 11340 DO j = 1, mask_size_l(mid,2) 11341 DO k = 1, mask_size_l(mid,3) 11342 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11343 ENDDO 11103 11344 ENDDO 11104 11345 ENDDO 11105 ENDDO 11346 ELSE 11347 DO i = 1, mask_size_l(mid,1) 11348 DO j = 1, mask_size_l(mid,2) 11349 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11350 grid ) 11351 DO k = 1, mask_size_l(mid,3) 11352 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11353 mask_j(mid,j), mask_i(mid,i) ) 11354 ENDDO 11355 ENDDO 11356 ENDDO 11357 ENDIF 11358 resorted = .TRUE. 11359 ELSE 11360 to_be_resorted => ldsa_av 11106 11361 ENDIF 11107 resorted = .TRUE. 11108 ELSE 11109 to_be_resorted => pm25_av 11110 ENDIF 11111 11112 CASE ( 'PM10' ) 11113 IF ( av == 0 ) THEN 11114 DO i = nxl, nxr 11115 DO j = nys, nyn 11116 DO k = nzb, nz_do3d 11117 temp_bin = 0.0_wp 11118 DO ib = 1, nbins_aerosol 11119 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp ) THEN 11120 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 11362 11363 CASE ( 'N_UFP' ) 11364 IF ( av == 0 ) THEN 11365 DO i = nxl, nxr 11366 DO j = nys, nyn 11367 DO k = nzb, nz_do3d 11368 temp_bin = 0.0_wp 11369 DO ib = 1, nbins_aerosol 11370 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp ) THEN 11371 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 11372 ENDIF 11373 ENDDO 11374 tend(k,j,i) = temp_bin 11375 ENDDO 11376 ENDDO 11377 ENDDO 11378 IF ( .NOT. mask_surface(mid) ) THEN 11379 DO i = 1, mask_size_l(mid,1) 11380 DO j = 1, mask_size_l(mid,2) 11381 DO k = 1, mask_size_l(mid,3) 11382 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11383 ENDDO 11384 ENDDO 11385 ENDDO 11386 ELSE 11387 DO i = 1, mask_size_l(mid,1) 11388 DO j = 1, mask_size_l(mid,2) 11389 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11390 grid ) 11391 DO k = 1, mask_size_l(mid,3) 11392 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11393 mask_j(mid,j), mask_i(mid,i) ) 11394 ENDDO 11395 ENDDO 11396 ENDDO 11397 ENDIF 11398 resorted = .TRUE. 11399 ELSE 11400 to_be_resorted => nufp_av 11401 ENDIF 11402 11403 CASE ( 'Ntot' ) 11404 IF ( av == 0 ) THEN 11405 DO i = nxl, nxr 11406 DO j = nys, nyn 11407 DO k = nzb, nz_do3d 11408 temp_bin = 0.0_wp 11409 DO ib = 1, nbins_aerosol 11410 temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i) 11411 ENDDO 11412 tend(k,j,i) = temp_bin 11413 ENDDO 11414 ENDDO 11415 ENDDO 11416 IF ( .NOT. mask_surface(mid) ) THEN 11417 DO i = 1, mask_size_l(mid,1) 11418 DO j = 1, mask_size_l(mid,2) 11419 DO k = 1, mask_size_l(mid,3) 11420 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11421 ENDDO 11422 ENDDO 11423 ENDDO 11424 ELSE 11425 DO i = 1, mask_size_l(mid,1) 11426 DO j = 1, mask_size_l(mid,2) 11427 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11428 grid ) 11429 DO k = 1, mask_size_l(mid,3) 11430 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11431 mask_j(mid,j), mask_i(mid,i) ) 11432 ENDDO 11433 ENDDO 11434 ENDDO 11435 ENDIF 11436 resorted = .TRUE. 11437 ELSE 11438 to_be_resorted => ntot_av 11439 ENDIF 11440 11441 CASE ( 'PM0.1' ) 11442 IF ( av == 0 ) THEN 11443 DO i = nxl, nxr 11444 DO j = nys, nyn 11445 DO k = nzb, nz_do3d 11446 temp_bin = 0.0_wp 11447 DO ib = 1, nbins_aerosol 11448 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp ) THEN 11449 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 11450 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 11451 ENDDO 11452 ENDIF 11453 ENDDO 11454 tend(k,j,i) = temp_bin 11455 ENDDO 11456 ENDDO 11457 ENDDO 11458 IF ( .NOT. mask_surface(mid) ) THEN 11459 DO i = 1, mask_size_l(mid,1) 11460 DO j = 1, mask_size_l(mid,2) 11461 DO k = 1, mask_size_l(mid,3) 11462 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11463 ENDDO 11464 ENDDO 11465 ENDDO 11466 ELSE 11467 DO i = 1, mask_size_l(mid,1) 11468 DO j = 1, mask_size_l(mid,2) 11469 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11470 grid ) 11471 DO k = 1, mask_size_l(mid,3) 11472 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11473 mask_j(mid,j), mask_i(mid,i) ) 11474 ENDDO 11475 ENDDO 11476 ENDDO 11477 ENDIF 11478 resorted = .TRUE. 11479 ELSE 11480 to_be_resorted => pm01_av 11481 ENDIF 11482 11483 CASE ( 'PM2.5' ) 11484 IF ( av == 0 ) THEN 11485 DO i = nxl, nxr 11486 DO j = nys, nyn 11487 DO k = nzb, nz_do3d 11488 temp_bin = 0.0_wp 11489 DO ib = 1, nbins_aerosol 11490 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp ) THEN 11491 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 11492 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 11493 ENDDO 11494 ENDIF 11495 ENDDO 11496 tend(k,j,i) = temp_bin 11497 ENDDO 11498 ENDDO 11499 ENDDO 11500 IF ( .NOT. mask_surface(mid) ) THEN 11501 DO i = 1, mask_size_l(mid,1) 11502 DO j = 1, mask_size_l(mid,2) 11503 DO k = 1, mask_size_l(mid,3) 11504 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11505 ENDDO 11506 ENDDO 11507 ENDDO 11508 ELSE 11509 DO i = 1, mask_size_l(mid,1) 11510 DO j = 1, mask_size_l(mid,2) 11511 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11512 grid ) 11513 DO k = 1, mask_size_l(mid,3) 11514 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11515 mask_j(mid,j), mask_i(mid,i) ) 11516 ENDDO 11517 ENDDO 11518 ENDDO 11519 ENDIF 11520 resorted = .TRUE. 11521 ELSE 11522 to_be_resorted => pm25_av 11523 ENDIF 11524 11525 CASE ( 'PM10' ) 11526 IF ( av == 0 ) THEN 11527 DO i = nxl, nxr 11528 DO j = nys, nyn 11529 DO k = nzb, nz_do3d 11530 temp_bin = 0.0_wp 11531 DO ib = 1, nbins_aerosol 11532 IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp ) THEN 11533 DO ic = ib, nbins_aerosol * ncc, nbins_aerosol 11534 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 11535 ENDDO 11536 ENDIF 11537 ENDDO 11538 tend(k,j,i) = temp_bin 11539 ENDDO 11540 ENDDO 11541 ENDDO 11542 IF ( .NOT. mask_surface(mid) ) THEN 11543 DO i = 1, mask_size_l(mid,1) 11544 DO j = 1, mask_size_l(mid,2) 11545 DO k = 1, mask_size_l(mid,3) 11546 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11547 ENDDO 11548 ENDDO 11549 ENDDO 11550 ELSE 11551 DO i = 1, mask_size_l(mid,1) 11552 DO j = 1, mask_size_l(mid,2) 11553 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11554 grid ) 11555 DO k = 1, mask_size_l(mid,3) 11556 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11557 mask_j(mid,j), mask_i(mid,i) ) 11558 ENDDO 11559 ENDDO 11560 ENDDO 11561 ENDIF 11562 resorted = .TRUE. 11563 ELSE 11564 to_be_resorted => pm10_av 11565 ENDIF 11566 11567 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 11568 IF ( av == 0 ) THEN 11569 IF ( is_used( prtcl, TRIM( variable(3:) ) ) ) THEN 11570 found_index = get_index( prtcl, TRIM( variable(3:) ) ) 11571 DO i = nxl, nxr 11572 DO j = nys, nyn 11573 DO k = nzb, nz_do3d 11574 temp_bin = 0.0_wp 11575 DO ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol 11121 11576 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 11122 11577 ENDDO 11123 ENDIF 11124 ENDDO 11125 tend(k,j,i) = temp_bin 11126 ENDDO 11127 ENDDO 11128 ENDDO 11129 IF ( .NOT. mask_surface(mid) ) THEN 11130 DO i = 1, mask_size_l(mid,1) 11131 DO j = 1, mask_size_l(mid,2) 11132 DO k = 1, mask_size_l(mid,3) 11133 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11578 tend(k,j,i) = temp_bin 11579 ENDDO 11134 11580 ENDDO 11135 11581 ENDDO 11136 E NDDO11137 ELSE11138 DO i = 1, mask_size_l(mid,1)11139 DO j = 1, mask_size_l(mid,2)11140 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), &11141 grid)11142 DO k = 1, mask_size_l(mid,3)11143 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), &11144 mask_j(mid,j), mask_i(mid,i) )11582 ELSE 11583 tend = 0.0_wp 11584 ENDIF 11585 IF ( .NOT. mask_surface(mid) ) THEN 11586 DO i = 1, mask_size_l(mid,1) 11587 DO j = 1, mask_size_l(mid,2) 11588 DO k = 1, mask_size_l(mid,3) 11589 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11590 ENDDO 11145 11591 ENDDO 11146 11592 ENDDO 11147 ENDDO 11593 ELSE 11594 DO i = 1, mask_size_l(mid,1) 11595 DO j = 1, mask_size_l(mid,2) 11596 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11597 grid ) 11598 DO k = 1, mask_size_l(mid,3) 11599 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11600 mask_j(mid,j), mask_i(mid,i) ) 11601 ENDDO 11602 ENDDO 11603 ENDDO 11604 ENDIF 11605 resorted = .TRUE. 11606 ELSE 11607 ! 11608 !-- 9: remove salsa_s_ from the beginning 11609 IF ( TRIM( variable(9:) ) == 'BC' ) to_be_resorted => s_bc_av 11610 IF ( TRIM( variable(9:) ) == 'DU' ) to_be_resorted => s_du_av 11611 IF ( TRIM( variable(9:) ) == 'NH' ) to_be_resorted => s_nh_av 11612 IF ( TRIM( variable(9:) ) == 'NO' ) to_be_resorted => s_no_av 11613 IF ( TRIM( variable(9:) ) == 'OC' ) to_be_resorted => s_oc_av 11614 IF ( TRIM( variable(9:) ) == 'SO4' ) to_be_resorted => s_so4_av 11615 IF ( TRIM( variable(9:) ) == 'SS' ) to_be_resorted => s_ss_av 11148 11616 ENDIF 11149 resorted = .TRUE. 11150 ELSE 11151 to_be_resorted => pm10_av 11152 ENDIF 11153 11154 CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' ) 11155 IF ( av == 0 ) THEN 11156 IF ( is_used( prtcl, TRIM( variable(3:) ) ) ) THEN 11157 found_index = get_index( prtcl, TRIM( variable(3:) ) ) 11617 11618 CASE ( 's_H2O' ) 11619 IF ( av == 0 ) THEN 11620 found_index = get_index( prtcl, 'H2O' ) 11158 11621 DO i = nxl, nxr 11159 11622 DO j = nys, nyn … … 11167 11630 ENDDO 11168 11631 ENDDO 11169 ELSE 11170 tend = 0.0_wp 11171 ENDIF 11172 IF ( .NOT. mask_surface(mid) ) THEN 11173 DO i = 1, mask_size_l(mid,1) 11174 DO j = 1, mask_size_l(mid,2) 11175 DO k = 1, mask_size_l(mid,3) 11176 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11632 IF ( .NOT. mask_surface(mid) ) THEN 11633 DO i = 1, mask_size_l(mid,1) 11634 DO j = 1, mask_size_l(mid,2) 11635 DO k = 1, mask_size_l(mid,3) 11636 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11637 ENDDO 11177 11638 ENDDO 11178 11639 ENDDO 11179 E NDDO11180 ELSE11181 DO i = 1, mask_size_l(mid,1)11182 DO j = 1, mask_size_l(mid,2)11183 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), &11184 grid)11185 DO k = 1, mask_size_l(mid,3)11186 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), &11187 mask_j(mid,j), mask_i(mid,i) )11640 ELSE 11641 DO i = 1, mask_size_l(mid,1) 11642 DO j = 1, mask_size_l(mid,2) 11643 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11644 grid ) 11645 DO k = 1, mask_size_l(mid,3) 11646 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11647 mask_j(mid,j), mask_i(mid,i) ) 11648 ENDDO 11188 11649 ENDDO 11189 11650 ENDDO 11190 ENDDO 11651 ENDIF 11652 resorted = .TRUE. 11653 ELSE 11654 to_be_resorted => s_h2o_av 11191 11655 ENDIF 11192 resorted = .TRUE. 11193 ELSE 11194 IF ( TRIM( variable(3:) ) == 'BC' ) to_be_resorted => s_bc_av 11195 IF ( TRIM( variable(3:) ) == 'DU' ) to_be_resorted => s_du_av 11196 IF ( TRIM( variable(3:) ) == 'NH' ) to_be_resorted => s_nh_av 11197 IF ( TRIM( variable(3:) ) == 'NO' ) to_be_resorted => s_no_av 11198 IF ( TRIM( variable(3:) ) == 'OC' ) to_be_resorted => s_oc_av 11199 IF ( TRIM( variable(3:) ) == 'SO4' ) to_be_resorted => s_so4_av 11200 IF ( TRIM( variable(3:) ) == 'SS' ) to_be_resorted => s_ss_av 11201 ENDIF 11202 11203 CASE ( 's_H2O' ) 11204 IF ( av == 0 ) THEN 11205 found_index = get_index( prtcl, 'H2O' ) 11206 DO i = nxl, nxr 11207 DO j = nys, nyn 11208 DO k = nzb, nz_do3d 11209 temp_bin = 0.0_wp 11210 DO ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol 11211 temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i) 11212 ENDDO 11213 tend(k,j,i) = temp_bin 11214 ENDDO 11215 ENDDO 11216 ENDDO 11217 IF ( .NOT. mask_surface(mid) ) THEN 11218 DO i = 1, mask_size_l(mid,1) 11219 DO j = 1, mask_size_l(mid,2) 11220 DO k = 1, mask_size_l(mid,3) 11221 local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11222 ENDDO 11223 ENDDO 11224 ENDDO 11225 ELSE 11226 DO i = 1, mask_size_l(mid,1) 11227 DO j = 1, mask_size_l(mid,2) 11228 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), & 11229 grid ) 11230 DO k = 1, mask_size_l(mid,3) 11231 local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ), & 11232 mask_j(mid,j), mask_i(mid,i) ) 11233 ENDDO 11234 ENDDO 11235 ENDDO 11236 ENDIF 11237 resorted = .TRUE. 11238 ELSE 11239 to_be_resorted => s_h2o_av 11240 ENDIF 11241 11242 CASE DEFAULT 11243 found = .FALSE. 11244 11245 END SELECT 11656 11657 CASE DEFAULT 11658 found = .FALSE. 11659 11660 END SELECT 11661 ENDIF 11246 11662 11247 11663 IF ( .NOT. resorted ) THEN … … 11252 11668 DO j = 1, mask_size_l(mid,2) 11253 11669 DO k = 1, mask_size_l(mid,3) 11254 local_pf(i,j,k) = to_be_resorted( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) )11670 local_pf(i,j,k) = to_be_resorted( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) 11255 11671 ENDDO 11256 11672 ENDDO … … 11293 11709 INTEGER(iwp), INTENT(inout) :: ncomp !< Number of components 11294 11710 11295 CHARACTER( len=3), INTENT(in) :: listcomp(nlist) !< List cof component names11711 CHARACTER(LEN=3), INTENT(in) :: listcomp(nlist) !< List cof component names 11296 11712 11297 11713 TYPE(component_index), INTENT(inout) :: self !< Object containing the indices of different … … 11329 11745 IMPLICIT NONE 11330 11746 11331 CHARACTER( len=*), INTENT(in) :: incomp !< Component name11747 CHARACTER(LEN=*), INTENT(in) :: incomp !< Component name 11332 11748 11333 11749 INTEGER(iwp) :: ii !< index … … 11359 11775 IMPLICIT NONE 11360 11776 11361 CHARACTER( len=*), INTENT(in) :: icomp !< Component name11777 CHARACTER(LEN=*), INTENT(in) :: icomp !< Component name 11362 11778 11363 11779 TYPE(component_index), INTENT(in) :: self !< Object containing the indices of different
Note: See TracChangeset
for help on using the changeset viewer.
