Ignore:
Timestamp:
May 31, 2019 3:19:05 PM (2 years ago)
Author:
monakurppa
Message:

remove salsa_util_mod.f90 and correct some bugs in salsa and salsa test case

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/salsa_mod.f90

    r3956 r4012  
    2626! -----------------
    2727! $Id$
     28! Merge salsa branch to trunk. List of changes:
     29! - Error corrected in distr_update that resulted in the aerosol number size
     30!   distribution not converging if the concentration was nclim.
     31! - Added a separate output for aerosol liquid water (s_H2O)
     32! - aerosol processes for a size bin are now calculated only if the aerosol
     33!   number of concentration of that bin is > 2*nclim
     34! - An initialisation error in the subroutine "deposition" corrected and the
     35!   subroutine reformatted.
     36! - stuff from salsa_util_mod.f90 moved into salsa_mod.f90
     37! - calls for closing the netcdf input files added
     38!
     39! 3956 2019-05-07 12:32:52Z monakurppa
    2840! - Conceptual bug in depo_surf correct for urban and land surface model
    2941! - Subroutine salsa_tendency_ij optimized.
     
    147159!> @todo emission mode "parameterized", i.e. based on street type
    148160!> @todo Allow insoluble emissions
    149 !> @todo two-way nesting is not working properly
     161!> @todo Apply flux limiter in prognostic equations
    150162!------------------------------------------------------------------------------!
    151163 MODULE salsa_mod
    152164
    153     USE basic_constants_and_equations_mod,                                     &
     165    USE basic_constants_and_equations_mod,                                                         &
    154166        ONLY:  c_p, g, p_0, pi, r_d
    155167
    156     USE chem_gasphase_mod,                                                     &
     168    USE chem_gasphase_mod,                                                                         &
    157169        ONLY:  nspec, nvar, spc_names
    158170
    159     USE chem_modules,                                                          &
     171    USE chem_modules,                                                                              &
    160172        ONLY:  call_chem_at_all_substeps, chem_gasphase_on, chem_species
    161173
    162174    USE control_parameters
    163175
    164     USE indices,                                                               &
    165         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
    166                nzb_s_inner, nz, nzt, wall_flags_0
     176    USE indices,                                                                                   &
     177        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nz, nzt, wall_flags_0
    167178
    168179    USE kinds
    169180
     181    USE netcdf_data_input_mod,                                                                     &
     182        ONLY:  chem_emis_att_type, chem_emis_val_type
     183
    170184    USE pegrid
    171185
    172     USE salsa_util_mod
    173 
    174     USE statistics,                                                            &
     186    USE statistics,                                                                                &
    175187        ONLY:  sums_salsa_ws_l
    176188
     
    452464                                                               (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
    453465    REAL(wp), DIMENSION(nmod) ::  surface_aerosol_flux = &
    454                              (/1.04e+11_wp, 3.23E+10_wp, 5.4E+6_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
     466                                 (/1.0E+8_wp, 1.0E+9_wp, 1.0E+5_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
    455467
    456468    REAL(wp), DIMENSION(:), ALLOCATABLE ::  bin_low_limits     !< to deliver information about
     
    464476!
    465477!-- SALSA derived datatypes:
     478!
     479!-- Component index
     480    TYPE component_index
     481       CHARACTER(len=3), ALLOCATABLE ::  comp(:)  !< Component name
     482       INTEGER(iwp) ::  ncomp  !< Number of components
     483       INTEGER(iwp), ALLOCATABLE ::  ind(:)  !< Component index
     484    END TYPE component_index
    466485!
    467486!-- For matching LSM and USM surface types and the deposition module surface types
     
    607626    TYPE(salsa_emission_value_type)     ::  aero_emission      !< emission values
    608627    TYPE(salsa_emission_mode_type)      ::  def_modes          !< default emission modes
     628
     629    TYPE(chem_emis_att_type) ::  chem_emission_att  !< chemistry emission attributes
     630    TYPE(chem_emis_val_type) ::  chem_emission      !< chemistry emission values
    609631
    610632    TYPE(t_section), DIMENSION(:), ALLOCATABLE ::  aero  !< local aerosol properties
     
    9791001 SUBROUTINE salsa_header( io )
    9801002
     1003    USE indices,                                                                                   &
     1004        ONLY:  nx, ny, nz
     1005
    9811006    IMPLICIT NONE
    9821007 
     
    9871012    WRITE( io, 2 ) skip_time_do_salsa
    9881013    WRITE( io, 3 ) dt_salsa
    989     WRITE( io, 4 )  SHAPE( aerosol_number(1)%conc ), nbins_aerosol
     1014    WRITE( io, 4 )  nz, ny, nx, nbins_aerosol
    9901015    IF ( advect_particle_water )  THEN
    991        WRITE( io, 5 )  SHAPE( aerosol_mass(1)%conc ), ncomponents_mass*nbins_aerosol,             &
     1016       WRITE( io, 5 )  SHAPE( aerosol_mass(1)%conc ), ncomponents_mass*nbins_aerosol,              &
    9921017                        advect_particle_water
    9931018    ELSE
     
    11541179                 aerosol_number(i)%init(nzb:nzt+1),                                                &
    11551180                 aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     1181       aerosol_number(i)%init = nclim
    11561182       IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
    11571183          ALLOCATE( aerosol_number(i)%source(nys:nyn,nxl:nxr) )
     
    11791205                 aerosol_mass(i)%init(nzb:nzt+1),                                                  &
    11801206                 aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1)  )
     1207       aerosol_mass(i)%init = mclim
    11811208       IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
    11821209          ALLOCATE( aerosol_mass(i)%source(nys:nyn,nxl:nxr) )
     
    12831310                    salsa_gas(i)%init(nzb:nzt+1),                              &
    12841311                    salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     1312          salsa_gas(i)%init = nclim
    12851313          IF ( include_emission )  ALLOCATE( salsa_gas(i)%source(nys:nys,nxl:nxr) )
    12861314       ENDDO
     
    14021430    IF ( nldepo )  sedim_vd = 0.0_wp
    14031431
    1404     DO  ib = 1, nbins_aerosol
    1405        IF ( .NOT. read_restart_data_salsa )  aerosol_number(ib)%conc = nclim
    1406        aerosol_number(ib)%conc_p    = 0.0_wp
    1407        aerosol_number(ib)%tconc_m   = 0.0_wp
    1408        aerosol_number(ib)%flux_s    = 0.0_wp
    1409        aerosol_number(ib)%diss_s    = 0.0_wp
    1410        aerosol_number(ib)%flux_l    = 0.0_wp
    1411        aerosol_number(ib)%diss_l    = 0.0_wp
    1412        aerosol_number(ib)%init      = nclim
    1413        aerosol_number(ib)%sums_ws_l = 0.0_wp
    1414     ENDDO
    1415     DO  ic = 1, ncomponents_mass*nbins_aerosol
    1416        IF ( .NOT. read_restart_data_salsa )  aerosol_mass(ic)%conc = mclim
    1417        aerosol_mass(ic)%conc_p    = 0.0_wp
    1418        aerosol_mass(ic)%tconc_m   = 0.0_wp
    1419        aerosol_mass(ic)%flux_s    = 0.0_wp
    1420        aerosol_mass(ic)%diss_s    = 0.0_wp
    1421        aerosol_mass(ic)%flux_l    = 0.0_wp
    1422        aerosol_mass(ic)%diss_l    = 0.0_wp
    1423        aerosol_mass(ic)%init      = mclim
    1424        aerosol_mass(ic)%sums_ws_l = 0.0_wp
    1425     ENDDO
    1426 
    14271432    IF ( .NOT. salsa_gases_from_chem )  THEN
     1433       IF ( .NOT. read_restart_data_salsa )  THEN
     1434          salsa_gas(1)%conc = h2so4_init
     1435          salsa_gas(2)%conc = hno3_init
     1436          salsa_gas(3)%conc = nh3_init
     1437          salsa_gas(4)%conc = ocnv_init
     1438          salsa_gas(5)%conc = ocsv_init
     1439       ENDIF
    14281440       DO  ig = 1, ngases_salsa
    14291441          salsa_gas(ig)%conc_p    = 0.0_wp
     
    14341446          salsa_gas(ig)%diss_l    = 0.0_wp
    14351447          salsa_gas(ig)%sums_ws_l = 0.0_wp
     1448          salsa_gas(ig)%conc_p    = salsa_gas(ig)%conc
    14361449       ENDDO
    1437        IF ( .NOT. read_restart_data_salsa )  THEN
    1438           salsa_gas(1)%conc = h2so4_init
    1439           salsa_gas(2)%conc = hno3_init
    1440           salsa_gas(3)%conc = nh3_init
    1441           salsa_gas(4)%conc = ocnv_init
    1442           salsa_gas(5)%conc = ocsv_init
    1443        ENDIF
    1444 !
    1445 !--    Set initial value for gas compound tracers and initial values
     1450!
     1451!--    Set initial value for gas compound tracer
    14461452       salsa_gas(1)%init = h2so4_init
    14471453       salsa_gas(2)%init = hno3_init
     
    14661472!-- Initialise location-dependent aerosol size distributions and chemical compositions:
    14671473    CALL aerosol_init
    1468 !
     1474
    14691475!-- Initalisation run of SALSA + calculate the vertical top index of the topography
    14701476    DO  i = nxl, nxr
     
    14771483       ENDDO
    14781484    ENDDO
     1485
     1486    DO  ib = 1, nbins_aerosol
     1487       aerosol_number(ib)%conc_p    = aerosol_number(ib)%conc
     1488       aerosol_number(ib)%tconc_m   = 0.0_wp
     1489       aerosol_number(ib)%flux_s    = 0.0_wp
     1490       aerosol_number(ib)%diss_s    = 0.0_wp
     1491       aerosol_number(ib)%flux_l    = 0.0_wp
     1492       aerosol_number(ib)%diss_l    = 0.0_wp
     1493       aerosol_number(ib)%sums_ws_l = 0.0_wp
     1494    ENDDO
     1495    DO  ic = 1, ncomponents_mass*nbins_aerosol
     1496       aerosol_mass(ic)%conc_p    = aerosol_mass(ic)%conc
     1497       aerosol_mass(ic)%tconc_m   = 0.0_wp
     1498       aerosol_mass(ic)%flux_s    = 0.0_wp
     1499       aerosol_mass(ic)%diss_s    = 0.0_wp
     1500       aerosol_mass(ic)%flux_l    = 0.0_wp
     1501       aerosol_mass(ic)%diss_l    = 0.0_wp
     1502       aerosol_mass(ic)%sums_ws_l = 0.0_wp
     1503    ENDDO
     1504!
    14791505!
    14801506!-- Initialise the deposition scheme and surface types
     
    15871613
    15881614    USE netcdf_data_input_mod,                                                                     &
    1589         ONLY:  get_attribute, get_variable, netcdf_data_input_get_dimension_length, open_read_file
     1615        ONLY:  close_input_file, get_attribute, get_variable,                                      &
     1616               netcdf_data_input_get_dimension_length, open_read_file
    15901617
    15911618    IMPLICIT NONE
     
    16151642
    16161643    REAL(wp), DIMENSION(nbins_aerosol) ::  core   !< size of the bin mid aerosol particle
    1617     REAL(wp), DIMENSION(nbins_aerosol) ::  nsect  !< size distribution (#/m3)
    16181644
    16191645    REAL(wp), DIMENSION(0:nz+1) ::  pnf2a   !< number fraction in 2a
     
    16381664!
    16391665!-- Set concentrations to zero
    1640     nsect(:)     = 0.0_wp
    16411666    pndist(:,:)  = 0.0_wp
    16421667    pnf2a(:)     = nf2a
     
    16671692!
    16681693!--       Allocate memory
    1669           ALLOCATE( pr_z(1:pr_nz), pr_mass_fracs_a(nzb:nzt+1,pr_ncc),                            &
     1694          ALLOCATE( pr_z(1:pr_nz), pr_mass_fracs_a(nzb:nzt+1,pr_ncc),                              &
    16701695                    pr_mass_fracs_b(nzb:nzt+1,pr_ncc) )
    16711696          pr_mass_fracs_a = 0.0_wp
     
    18061831                           ' for SALSA missing!'
    18071832          CALL message( 'salsa_mod: aerosol_init', 'PA0607', 1, 2, 0, 6, 0 )
    1808 
     1833!
     1834!--       Close input file
     1835          CALL close_input_file( id_dyn )
    18091836       ENDIF   ! netcdf_extend
    18101837
     
    18921919             salsa_gas(ig)%init(nzb)   =  salsa_gas(ig)%init(nzb+1)
    18931920             salsa_gas(ig)%init(nzt+1) =  salsa_gas(ig)%init(nzt)
    1894              DO  k = nzb, nzt+1
    1895                 salsa_gas(ig)%conc(k,:,:) = salsa_gas(ig)%init(k)
    1896              ENDDO
     1921             IF ( .NOT. read_restart_data_salsa )  THEN
     1922                DO  k = nzb, nzt+1
     1923                   salsa_gas(ig)%conc(k,:,:) = salsa_gas(ig)%init(k)
     1924                ENDDO
     1925             ENDIF
    18971926          ENDDO
    18981927
     
    19011930                           ' for SALSA missing!'
    19021931          CALL message( 'salsa_mod: aerosol_init', 'PA0610', 1, 2, 0, 6, 0 )
     1932!
     1933!--       Close input file
     1934          CALL close_input_file( id_dyn )
    19031935       ENDIF   ! netcdf_extend
    19041936#else
     
    19401972!--          Region 1:
    19411973             DO  ib = start_subrange_1a, end_subrange_1a
    1942                 aerosol_number(ib)%conc(k,j,i) = pndist(k,ib) * flag
     1974                IF ( .NOT. read_restart_data_salsa )  THEN
     1975                   aerosol_number(ib)%conc(k,j,i) = pndist(k,ib) * flag
     1976                ENDIF
    19431977                IF ( prunmode == 1 )  THEN
    19441978                   aerosol_number(ib)%init = pndist(:,ib)
     
    19491983             IF ( nreg > 1 )  THEN
    19501984                DO  ib = start_subrange_2a, end_subrange_2a
    1951                    aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) * pndist(k,ib) * flag
     1985                   IF ( .NOT. read_restart_data_salsa )  THEN
     1986                      aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) * pndist(k,ib) * flag
     1987                   ENDIF
    19521988                   IF ( prunmode == 1 )  THEN
    19531989                      aerosol_number(ib)%init = MAX( 0.0_wp, nf2a ) * pndist(:,ib)
     
    19571993                   DO  ib = start_subrange_2b, end_subrange_2b
    19581994                      IF ( pnf2a(k) < 1.0_wp )  THEN
    1959                          aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pnf2a(k) ) *       &
    1960                                                           pndist(k,ib) * flag
     1995                         IF ( .NOT. read_restart_data_salsa )  THEN
     1996                            aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pnf2a(k) ) *    &
     1997                                                             pndist(k,ib) * flag
     1998                         ENDIF
    19611999                         IF ( prunmode == 1 )  THEN
    19622000                            aerosol_number(ib)%init = MAX( 0.0_wp, 1.0_wp - nf2a ) * pndist(:,ib)
     
    19762014                ib = start_subrange_1a
    19772015                DO  ic = ss, ee
    1978                    aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) * pndist(k,ib)&
    1979                                                   * core(ib) * arhoh2so4 * flag
     2016                   IF ( .NOT. read_restart_data_salsa )  THEN
     2017                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) *          &
     2018                                                     pndist(k,ib) * core(ib) * arhoh2so4 * flag
     2019                   ENDIF
    19802020                   IF ( prunmode == 1 )  THEN
    19812021                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) * pndist(k,ib) &
     
    19912031                ee = ( index_oc - 1 ) * nbins_aerosol + end_subrange_1a !< end
    19922032                ib = start_subrange_1a
    1993                 DO  ic = ss, ee
    1994                    aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *       &
    1995                                                   core(ib) * arhooc * flag
     2033                DO  ic = ss, ee
     2034                   IF ( .NOT. read_restart_data_salsa )  THEN
     2035                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *    &
     2036                                                     core(ib) * arhooc * flag
     2037                   ENDIF
    19962038                   IF ( prunmode == 1 )  THEN
    19972039                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *        &
     
    21432185             ib = start_subrange_2a
    21442186             DO ic = ss, ee
    2145                 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib) * &
    2146                                               pcore(ib) * prho * flag
     2187                IF ( .NOT. read_restart_data_salsa )  THEN
     2188                   aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib)&
     2189                                                  * pcore(ib) * prho * flag
     2190                ENDIF
    21472191                IF ( prunmode == 1 )  THEN
    21482192                   aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib) *  &
     
    21582202                ib = start_subrange_2a
    21592203                DO ic = ss, ee
    2160                    aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k) ) *&
    2161                                                   pndist(k,ib) * pcore(ib) * prho * flag
     2204                   IF ( .NOT. read_restart_data_salsa )  THEN
     2205                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k))&
     2206                                                     * pndist(k,ib) * pcore(ib) * prho * flag
     2207                   ENDIF
    21622208                   IF ( prunmode == 1 )  THEN
    21632209                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k) ) * &
     
    28102856       ENDDO
    28112857!
    2812 !--    On EACH call of salsa_driver, calculate the ambient sizes of
    2813 !--    particles by equilibrating soluble fraction of particles with water
    2814 !--    using the ZSR method.
     2858!--    Calculate the ambient sizes of particles by equilibrating soluble fraction of particles with
     2859!--    water using the ZSR method.
    28152860       in_rh = in_cw(k) / in_cs(k)
    28162861       IF ( prunmode==1  .OR.  .NOT. advect_particle_water )  THEN
     
    28652910!--    Calculate changes in concentrations
    28662911       DO ib = 1, nbins_aerosol
    2867           aerosol_number(ib)%conc(k,j,i) = aerosol_number(ib)%conc(k,j,i) + ( lo_aero(ib)%numc -      &
     2912          aerosol_number(ib)%conc(k,j,i) = aerosol_number(ib)%conc(k,j,i) + ( lo_aero(ib)%numc -   &
    28682913                                           aero_old(ib)%numc ) * flag
    28692914       ENDDO
     
    28752920          ic = 1
    28762921          DO ss = str, endi
    2877              aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
     2922             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
    28782923                                            aero_old(ic)%volc(vc) ) * arhoh2so4 * flag
    28792924             ic = ic+1
     
    28872932          ic = 1
    28882933          DO ss = str, endi
    2889              aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
     2934             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
    28902935                                            aero_old(ic)%volc(vc) ) * arhooc * flag
    28912936             ic = ic+1
     
    28992944          ic = 1 + end_subrange_1a
    29002945          DO ss = str, endi
    2901              aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
     2946             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
    29022947                                            aero_old(ic)%volc(vc) ) * arhobc * flag
    29032948             ic = ic+1
     
    29112956          ic = 1 + end_subrange_1a
    29122957          DO ss = str, endi
    2913              aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
     2958             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
    29142959                                            aero_old(ic)%volc(vc) ) * arhodu * flag
    29152960             ic = ic+1
     
    29232968          ic = 1 + end_subrange_1a
    29242969          DO ss = str, endi
    2925              aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
     2970             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
    29262971                                            aero_old(ic)%volc(vc) ) * arhoss * flag
    29272972             ic = ic+1
     
    29352980          ic = 1
    29362981          DO ss = str, endi
    2937              aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
     2982             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
    29382983                                            aero_old(ic)%volc(vc) ) * arhohno3 * flag
    29392984             ic = ic+1
     
    29472992          ic = 1
    29482993          DO ss = str, endi
    2949              aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
     2994             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
    29502995                                            aero_old(ic)%volc(vc) ) * arhonh3 * flag
    29512996             ic = ic+1
     
    29603005          ic = 1
    29613006          DO ss = str, endi
    2962              aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
     3007             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
    29633008                                            aero_old(ic)%volc(vc) ) * arhoh2o * flag
    2964              IF ( prunmode == 1 )  THEN
    2965                 aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k),                          &
    2966                                                 aerosol_mass(ss)%conc(k,j,i) )
    2967                 IF ( k == nzb+1 )  THEN
    2968                    aerosol_mass(ss)%init(k-1) = 0.0_wp
    2969                 ELSEIF ( k == nzt  )  THEN
    2970                    aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k)
    2971                 ENDIF
     3009             ic = ic+1
     3010          ENDDO
     3011       ENDIF
     3012       IF ( prunmode == 1 )  THEN
     3013          nc_h2o = get_index( prtcl,'H2O' )
     3014          vc = 8
     3015          str = ( nc_h2o-1 ) * nbins_aerosol + 1
     3016          endi = nc_h2o * nbins_aerosol
     3017          ic = 1
     3018          DO ss = str, endi
     3019             aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), ( lo_aero(ic)%volc(vc) - &
     3020                                             aero_old(ic)%volc(vc) ) * arhoh2o )
     3021             IF ( k == nzb+1 )  THEN
     3022                aerosol_mass(ss)%init(k-1) = aerosol_mass(ss)%init(k)
     3023             ELSEIF ( k == nzt  )  THEN
     3024                aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k)
     3025                aerosol_mass(ss)%conc(k+1,j,i) = aerosol_mass(ss)%init(k)
    29723026             ENDIF
    29733027             ic = ic+1
     
    34463500    REAL(wp) ::  avis              !< molecular viscocity of air (kg/(m*s))
    34473501    REAL(wp) ::  beta_im           !< parameter for turbulent impaction
    3448     REAL(wp) ::  beta              !< Cunningham slip-flow correction factor
    34493502    REAL(wp) ::  c_brownian_diff   !< coefficient for Brownian diffusion
    34503503    REAL(wp) ::  c_impaction       !< coefficient for inertial impaction
     
    34533506    REAL(wp) ::  depo              !< deposition velocity (m/s)
    34543507    REAL(wp) ::  gamma             !< parameter, Table 3 in Z01
    3455     REAL(wp) ::  Kn                !< Knudsen number
    34563508    REAL(wp) ::  lambda            !< molecular mean free path (m)
    34573509    REAL(wp) ::  mdiff             !< particle diffusivity coefficient
     
    34613513    REAL(wp) ::  pdn               !< particle density (kg/m3)
    34623514    REAL(wp) ::  ustar             !< friction velocity (m/s)
    3463     REAL(wp) ::  va                !< thermal speed of an air molecule (m/s)
    3464     REAL(wp) ::  zdwet             !< wet diameter (m)
     3515    REAL(wp) ::  va                !< thermal speed of an air molecule (m/s)
    34653516
    34663517    REAL(wp), INTENT(in) ::  adn    !< air density (kg/m3)
     
    34713522    REAL(wp), INTENT(inout) ::  kvis   !< kinematic viscosity of air (m2/s)
    34723523
     3524    REAL(wp), DIMENSION(nbins_aerosol) ::  beta   !< Cunningham slip-flow correction factor
     3525    REAL(wp), DIMENSION(nbins_aerosol) ::  Kn     !< Knudsen number
     3526    REAL(wp), DIMENSION(nbins_aerosol) ::  zdwet  !< wet diameter (m)
     3527
    34733528    REAL(wp), DIMENSION(:), INTENT(inout) ::  schmidt_num  !< particle Schmidt number
    34743529    REAL(wp), DIMENSION(:), INTENT(inout) ::  vc  !< critical fall speed i.e. settling velocity of
     
    34943549    lambda = 2.0_wp * avis / ( adn * va )
    34953550!
    3496 !-- Parameters for the land use category 'deciduous broadleaf trees'(Table 3)
    3497     alpha   = alpha_z01(depo_pcm_type_num)
    3498     gamma   = gamma_z01(depo_pcm_type_num)
    3499     par_a   = A_z01(depo_pcm_type_num, season) * 1.0E-3_wp
    3500 !
    3501 !-- Deposition efficiencies from Table 1. Constants from Table 2.
    3502     par_l            = l_p10(depo_pcm_type_num) * 0.01_wp
    3503     c_brownian_diff  = c_b_p10(depo_pcm_type_num)
    3504     c_interception   = c_in_p10(depo_pcm_type_num)
    3505     c_impaction      = c_im_p10(depo_pcm_type_num)
    3506     beta_im          = beta_im_p10(depo_pcm_type_num)
    3507     c_turb_impaction = c_it_p10(depo_pcm_type_num)
    3508 
    3509     DO  ib = 1, nbins_aerosol
    3510 
    3511        IF ( paero(ib)%numc < nclim )  CYCLE
    3512        zdwet = paero(ib)%dwet
    3513 !
    3514 !--    Knudsen number (Eq. 15.23)
    3515        Kn = MAX( 1.0E-2_wp, lambda / ( zdwet * 0.5_wp ) ) ! To avoid underflow
    3516 !
    3517 !--    Cunningham slip-flow correction (Eq. 15.30)
    3518        beta = 1.0_wp + Kn * ( 1.249_wp + 0.42_wp * EXP( -0.87_wp / Kn ) )
    3519 
    3520 !--    Particle diffusivity coefficient (Eq. 15.29)
    3521        mdiff = ( abo * tk * beta ) / ( 3.0_wp * pi * avis * zdwet )
    3522 !
    3523 !--    Particle Schmidt number (Eq. 15.36)
    3524        schmidt_num(ib) = kvis / mdiff
    3525 !
    3526 !--    Critical fall speed i.e. settling velocity  (Eq. 20.4)
    3527        vc(ib) = MIN( 1.0_wp, terminal_vel( 0.5_wp * zdwet, pdn, adn, avis, beta) )
    3528 !
    3529 !--    Friction velocity for deposition on vegetation. Calculated following Prandtl (1925):
    3530        IF ( lsdepo_pcm  .AND.  plant_canopy  .AND.  lad > 0.0_wp )  THEN
     3551!-- Particle wet diameter (m)
     3552    zdwet = paero(:)%dwet
     3553!
     3554!-- Knudsen number (Eq. 15.23)
     3555    Kn = MAX( 1.0E-2_wp, lambda / ( zdwet * 0.5_wp ) ) ! To avoid underflow
     3556!
     3557!-- Cunningham slip-flow correction (Eq. 15.30)
     3558    beta = 1.0_wp + Kn * ( 1.249_wp + 0.42_wp * EXP( -0.87_wp / Kn ) )
     3559!
     3560!-- Critical fall speed i.e. settling velocity  (Eq. 20.4)
     3561    vc = MIN( 1.0_wp, zdwet**2 * ( pdn - adn ) * g * beta / ( 18.0_wp * avis ) )
     3562!
     3563!-- Deposition on vegetation
     3564    IF ( lsdepo_pcm  .AND.  plant_canopy  .AND.  lad > 0.0_wp )  THEN
     3565!
     3566!--    Parameters for the land use category 'deciduous broadleaf trees'(Table 3)
     3567       alpha   = alpha_z01(depo_pcm_type_num)
     3568       gamma   = gamma_z01(depo_pcm_type_num)
     3569       par_a   = A_z01(depo_pcm_type_num, season) * 1.0E-3_wp
     3570!
     3571!--    Deposition efficiencies from Table 1. Constants from Table 2.
     3572       par_l            = l_p10(depo_pcm_type_num) * 0.01_wp
     3573       c_brownian_diff  = c_b_p10(depo_pcm_type_num)
     3574       c_interception   = c_in_p10(depo_pcm_type_num)
     3575       c_impaction      = c_im_p10(depo_pcm_type_num)
     3576       beta_im          = beta_im_p10(depo_pcm_type_num)
     3577       c_turb_impaction = c_it_p10(depo_pcm_type_num)
     3578
     3579       DO  ib = 1, nbins_aerosol
     3580
     3581          IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
     3582
     3583!--       Particle diffusivity coefficient (Eq. 15.29)
     3584          mdiff = ( abo * tk * beta(ib) ) / ( 3.0_wp * pi * avis * zdwet(ib) )
     3585!
     3586!--       Particle Schmidt number (Eq. 15.36)
     3587          schmidt_num(ib) = kvis / mdiff
     3588!
     3589!--       Friction velocity for deposition on vegetation. Calculated following Prandtl (1925):
    35313590          ustar = SQRT( cdc ) * mag_u
    35323591          SELECT CASE ( depo_pcm_par_num )
     
    35463605             paero(ib)%volc(ic) = paero(ib)%volc(ic) - depo * lad * paero(ib)%volc(ic) * dt_salsa
    35473606          ENDDO
    3548        ENDIF
    3549     ENDDO
     3607       ENDDO
     3608
     3609    ENDIF
    35503610
    35513611 END SUBROUTINE deposition
     
    37793839
    37803840             DO  ib = 1, nbins_aerosol
    3781                 IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp ) &
    3782                    CYCLE
     3841                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                    &
     3842                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
    37833843
    37843844                SELECT CASE ( depo_surf_par_num )
     
    38103870
    38113871             DO  ib = 1, nbins_aerosol
    3812                 IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp ) &
    3813                    CYCLE
     3872                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                    &
     3873                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
    38143874
    38153875                SELECT CASE ( depo_surf_par_num )
     
    38413901
    38423902             DO  ib = 1, nbins_aerosol
    3843                 IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp ) &
    3844                    CYCLE
     3903                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                    &
     3904                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
    38453905
    38463906                SELECT CASE ( depo_surf_par_num )
     
    38603920
    38613921          DO  ib = 1, nbins_aerosol
     3922             IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim ) )  CYCLE
    38623923!
    38633924!--          Calculate changes in surface fluxes due to dry deposition
     
    38913952
    38923953          DO  ib = 1, nbins_aerosol
    3893              IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp )    &
    3894                 CYCLE
     3954             IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                        &
     3955                  schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
    38953956
    38963957             SELECT CASE ( depo_surf_par_num )
     
    39323993! Description:
    39333994! ------------
    3934 ! Function for calculating terminal velocities for different particles sizes.
    3935 !------------------------------------------------------------------------------!
    3936  REAL(wp) FUNCTION terminal_vel( radius, rhop, rhoa, visc, beta )
    3937 
    3938     IMPLICIT NONE
    3939 
    3940     REAL(wp), INTENT(in) ::  beta    !< Cunningham correction factor
    3941     REAL(wp), INTENT(in) ::  radius  !< particle radius (m)
    3942     REAL(wp), INTENT(in) ::  rhop    !< particle density (kg/m3)
    3943     REAL(wp), INTENT(in) ::  rhoa    !< air density (kg/m3)
    3944     REAL(wp), INTENT(in) ::  visc    !< molecular viscosity of air (kg/(m*s))
    3945 
    3946 !
    3947 !-- Stokes law with Cunningham slip correction factor
    3948     terminal_vel = 4.0_wp * radius**2 * ( rhop - rhoa ) * g * beta / ( 18.0_wp * visc ) ! (m/s)
    3949 
    3950  END FUNCTION terminal_vel
    3951 
    3952 !------------------------------------------------------------------------------!
    3953 ! Description:
    3954 ! ------------
    39553995!> Calculates particle loss and change in size distribution due to (Brownian)
    39563996!> coagulation. Only for particles with dwet < 30 micrometres.
     
    40344074!-- Aero-aero coagulation
    40354075    DO  mm = 1, end_subrange_2b   ! smaller colliding particle
    4036        IF ( paero(mm)%numc < nclim )  CYCLE
     4076       IF ( paero(mm)%numc < ( 2.0_wp * nclim ) )  CYCLE
    40374077       DO  nn = mm, end_subrange_2b   ! larger colliding particle
    4038           IF ( paero(nn)%numc < nclim )  CYCLE
     4078          IF ( paero(nn)%numc < ( 2.0_wp * nclim ) )  CYCLE
    40394079
    40404080          zdpart_mm = MIN( paero(mm)%dwet, 30.0E-6_wp )     ! Limit to 30 um
     
    40534093!-- Aerosols in subrange 1a:
    40544094    DO  ib = start_subrange_1a, end_subrange_1a
    4055        IF ( paero(ib)%numc < nclim )  CYCLE
     4095       IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
    40564096       zminusterm   = 0.0_wp
    40574097       zplusterm(:) = 0.0_wp
     
    40804120!-- Aerosols in subrange 2a:
    40814121    DO  ib = start_subrange_2a, end_subrange_2a
    4082        IF ( paero(ib)%numc < nclim )  CYCLE
     4122       IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
    40834123       zminusterm   = 0.0_wp
    40844124       zplusterm(:) = 0.0_wp
     
    41284168    IF ( .NOT. no_insoluble )  THEN
    41294169       DO  ib = start_subrange_2b, end_subrange_2b
    4130           IF ( paero(ib)%numc < nclim )  CYCLE
     4170          IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
    41314171          zminusterm   = 0.0_wp
    41324172          zplusterm(:) = 0.0_wp
     
    43374377    REAL(wp), INTENT(inout) ::  pcw      !< Water vapor concentration (kg/m3)
    43384378
    4339     REAL(wp), DIMENSION(nbins_aerosol)               ::  zbeta          !< transitional correction factor
    4340     REAL(wp), DIMENSION(nbins_aerosol)               ::  zcolrate       !< collision rate (1/s)
    4341     REAL(wp), DIMENSION(nbins_aerosol)               ::  zcolrate_ocnv  !< collision rate of non-vol. OC (1/s)
    4342     REAL(wp), DIMENSION(start_subrange_1a+1) ::  zdfpart        !< particle diffusion coefficient (m2/s)
    4343     REAL(wp), DIMENSION(nbins_aerosol)               ::  zdvoloc        !< change of organics volume
    4344     REAL(wp), DIMENSION(nbins_aerosol)               ::  zdvolsa        !< change of sulphate volume
     4379    REAL(wp), DIMENSION(nbins_aerosol)       ::  zbeta          !< transitional correction factor
     4380    REAL(wp), DIMENSION(nbins_aerosol)       ::  zcolrate       !< collision rate (1/s)
     4381    REAL(wp), DIMENSION(nbins_aerosol)       ::  zcolrate_ocnv  !< collision rate of OCNV (1/s)
     4382    REAL(wp), DIMENSION(start_subrange_1a+1) ::  zdfpart        !< particle diffusion coef. (m2/s)
     4383    REAL(wp), DIMENSION(nbins_aerosol)       ::  zdvoloc        !< change of organics volume
     4384    REAL(wp), DIMENSION(nbins_aerosol)       ::  zdvolsa        !< change of sulphate volume
    43454385    REAL(wp), DIMENSION(2)                   ::  zj3n3          !< Formation massrate of molecules
    43464386                                                                !< in nucleation, (molec/m3s),
    43474387                                                                !< 1: H2SO4 and 2: organic vapor
    4348     REAL(wp), DIMENSION(nbins_aerosol)   ::  zknud          !< particle Knudsen number
     4388    REAL(wp), DIMENSION(nbins_aerosol)       ::  zknud          !< particle Knudsen number
    43494389
    43504390    TYPE(component_index), INTENT(in) :: prtcl  !< Keeps track which substances are used
     
    70957135    REAL(wp) ::  znfrac  !< number fraction to be moved to the larger bin
    70967136    REAL(wp) ::  zvfrac  !< volume fraction to be moved to the larger bin
    7097     REAL(wp) ::  zVexc   !< Volume in the grown bin which exceeds the bin upper limit
    7098     REAL(wp) ::  zVihi   !< particle volume at the high end of the bin
    7099     REAL(wp) ::  zVilo   !< particle volume at the low end of the bin
     7137    REAL(wp) ::  zvexc   !< Volume in the grown bin which exceeds the bin upper limit
     7138    REAL(wp) ::  zvihi   !< particle volume at the high end of the bin
     7139    REAL(wp) ::  zvilo   !< particle volume at the low end of the bin
    71007140    REAL(wp) ::  zvpart  !< particle volume (m3)
    7101     REAL(wp) ::  zVrat   !< volume ratio of a size bin
     7141    REAL(wp) ::  zvrat   !< volume ratio of a size bin
     7142
     7143    real(wp), dimension(nbins_aerosol) ::  dummy
    71027144
    71037145    TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero !< aerosol properties
     
    71067148    zvfrac      = 0.0_wp
    71077149    within_bins = .FALSE.
     7150
     7151    dummy = paero(:)%numc
    71087152!
    71097153!-- Check if the volume of the bin is within bin limits after update
     
    71507194!
    71517195!--          Volume ratio of the size bin
    7152              zVrat = paero(ib)%vhilim / paero(ib)%vlolim
     7196             zvrat = paero(ib)%vhilim / paero(ib)%vlolim
    71537197!
    71547198!--          Particle volume at the low end of the bin
    7155              zVilo = 2.0_wp * zvpart / ( 1.0_wp + zVrat )
     7199             zvilo = 2.0_wp * zvpart / ( 1.0_wp + zvrat )
    71567200!
    71577201!--          Particle volume at the high end of the bin
    7158              zVihi = zVrat * zVilo
     7202             zvihi = zvrat * zvilo
    71597203!
    71607204!--          Volume in the grown bin which exceeds the bin upper limit
    7161              zVexc = 0.5_wp * ( zVihi + paero(ib)%vhilim )
     7205             zvexc = 0.5_wp * ( zvihi + paero(ib)%vhilim )
    71627206!
    71637207!--          Number fraction to be moved to the larger bin
    7164              znfrac = MIN( 1.0_wp, ( zVihi - paero(ib)%vhilim) / ( zVihi - zVilo ) )
     7208             znfrac = MIN( 1.0_wp, ( zvihi - paero(ib)%vhilim) / ( zvihi - zvilo ) )
    71657209!
    71667210!--          Volume fraction to be moved to the larger bin
    7167              zvfrac = MIN( 0.99_wp, znfrac * zVexc / zvpart )
     7211             zvfrac = MIN( 0.99_wp, znfrac * zvexc / zvpart )
    71687212             IF ( zvfrac < 0.0_wp )  THEN
    71697213                message_string = 'Error: zvfrac < 0'
     
    71757219!
    71767220!--          Volume (cm3/cm3)
    7177              paero(mm)%volc(:) = paero(mm)%volc(:) + znfrac * paero(ib)%numc * zVexc *             &
    7178                                  paero(ib)%volc(:) / SUM( paero(ib)%volc(:) )
    7179              paero(ib)%volc(:) = paero(ib)%volc(:) - znfrac * paero(ib)%numc * zVexc *             &
    7180                                  paero(ib)%volc(:) / SUM( paero(ib)%volc(:) )
     7221             paero(mm)%volc(:) = paero(mm)%volc(:) + znfrac * paero(ib)%numc * zvexc *             &
     7222                                 paero(ib)%volc(:) / SUM( paero(ib)%volc(1:7) )
     7223             paero(ib)%volc(:) = paero(ib)%volc(:) - znfrac * paero(ib)%numc * zvexc *             &
     7224                                 paero(ib)%volc(:) / SUM( paero(ib)%volc(1:7) )
    71817225
    71827226!--          Number concentration (#/m3)
     
    72667310!--    Calculate total mass concentration per bin
    72677311       mcsum = 0.0_wp
    7268        DO  ic = 1, ncc
     7312       DO  ic = 1, ncomponents_mass
    72697313          icc = ( ic - 1 ) * nbins_aerosol + ib
    72707314          mcsum = mcsum + aerosol_mass(icc)%conc(:,j,i) * flag
     7315          aerosol_mass(icc)%conc(:,j,i) = MAX( mclim, aerosol_mass(icc)%conc(:,j,i) ) * flag
    72717316       ENDDO
    72727317!
    72737318!--    Check that number and mass concentration match qualitatively
    7274        IF ( ANY ( aerosol_number(ib)%conc(:,j,i) >= nclim  .AND. mcsum <= 0.0_wp ) )  THEN
     7319       IF ( ANY( aerosol_number(ib)%conc(:,j,i) > nclim  .AND. mcsum <= 0.0_wp ) )  THEN
    72757320          DO  k = nzb+1, nzt
    72767321             IF ( aerosol_number(ib)%conc(k,j,i) >= nclim  .AND. mcsum(k) <= 0.0_wp )  THEN
     
    73497394                                                 flag_zddry > 0.0_wp )
    73507395       ENDDO
    7351        aerosol_number(ib)%conc(:,j,i) = MERGE( nclim, aerosol_number(ib)%conc(:,j,i),              &
     7396       aerosol_number(ib)%conc(:,j,i) = MERGE( nclim * flag, aerosol_number(ib)%conc(:,j,i),       &
    73527397                                               flag_zddry > 0.0_wp )
    73537398       ra_dry(:,j,i,ib) = MAX( 1.0E-10_wp, 0.5_wp * zddry )
     
    83878432
    83888433    USE netcdf_data_input_mod,                                                                     &
    8389         ONLY:  check_existence, get_attribute, get_variable, inquire_num_variables,                &
    8390                inquire_variable_names, netcdf_data_input_get_dimension_length, open_read_file
     8434        ONLY:  check_existence, close_input_file, get_attribute, get_variable,                     &
     8435               inquire_num_variables, inquire_variable_names,                                      &
     8436               netcdf_data_input_get_dimension_length, open_read_file
    83918437
    83928438    USE surface_mod,                                                                               &
     
    84828528
    84838529#if defined( __netcdf )
     8530!
     8531!--       Check existence of PIDS_SALSA file
     8532          INQUIRE( FILE = TRIM( input_file_salsa ) // TRIM( coupling_char ), EXIST = netcdf_extend )
     8533          IF ( .NOT. netcdf_extend )  THEN
     8534             message_string = 'Input file '// TRIM( input_file_salsa ) //  TRIM( coupling_char )&
     8535                              // ' missing!'
     8536             CALL message( 'salsa_emission_setup', 'PA0629', 1, 2, 0, 6, 0 )
     8537          ENDIF
     8538!
     8539!--       Open file in read-only mode
     8540          CALL open_read_file( TRIM( input_file_salsa ) // TRIM( coupling_char ), id_salsa )
     8541
    84848542          IF ( init )  THEN
    8485 !
    8486 !--          Check existence of PIDS_SALSA file
    8487              INQUIRE( FILE = TRIM( input_file_salsa ) // TRIM( coupling_char ),                    &
    8488                       EXIST = netcdf_extend )
    8489              IF ( .NOT. netcdf_extend )  THEN
    8490                 message_string = 'Input file '// TRIM( input_file_salsa ) //  TRIM( coupling_char )&
    8491                                  // ' missing!'
    8492                 CALL message( 'salsa_emission_setup', 'PA0629', 1, 2, 0, 6, 0 )
    8493              ENDIF
    8494 !
    8495 !--          Open file in read-only mode
    8496              CALL open_read_file( TRIM( input_file_salsa ) // TRIM( coupling_char ), id_salsa )
    84978543!
    84988544!--          Read the index and name of chemical components
     
    88228868
    88238869
    8824              DEALLOCATE( source_array )
     8870             DEALLOCATE( nsect_emission, source_array )
    88258871!
    88268872!--       Pre-processed:
     
    88578903
    88588904          ENDIF
    8859 
     8905!
     8906!--       Close input file
     8907          CALL close_input_file( id_salsa )
    88608908#else
    88618909          message_string = 'salsa_emission_mode = "read_from_file", but preprocessor directive ' //&
     
    90739121
    90749122    USE netcdf_data_input_mod,                                                                     &
    9075         ONLY:  check_existence, chem_emis_att_type, chem_emis_val_type, get_attribute,             &
    9076                get_variable, inquire_num_variables, inquire_variable_names,                        &
     9123        ONLY:  check_existence, close_input_file, get_attribute, get_variable,                     &
     9124               inquire_num_variables, inquire_variable_names,                                      &
    90779125               netcdf_data_input_get_dimension_length, open_read_file
    90789126
     
    90989146    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_factor  !< emission time factor
    90999147
    9100     TYPE(chem_emis_att_type) ::  chem_emission_att  !< chemistry emission attributes
    9101     TYPE(chem_emis_val_type) ::  chem_emission      !< chemistry emission values
    9102 
    91039148!
    91049149!-- Reset surface fluxes
     
    91089153
    91099154#if defined( __netcdf )
     9155!
     9156!-- Check existence of PIDS_CHEM file
     9157    INQUIRE( FILE = 'PIDS_CHEM' // TRIM( coupling_char ), EXIST = netcdf_extend )
     9158    IF ( .NOT. netcdf_extend )  THEN
     9159       message_string = 'Input file PIDS_CHEM' //  TRIM( coupling_char ) // ' missing!'
     9160       CALL message( 'salsa_gas_emission_setup', 'PA0640', 1, 2, 0, 6, 0 )
     9161    ENDIF
     9162!
     9163!-- Open file in read-only mode
     9164    CALL open_read_file( 'PIDS_CHEM' // TRIM( coupling_char ), id_chem )
     9165
    91109166    IF ( init )  THEN
    9111 !
    9112 !--    Check existence of PIDS_CHEM file
    9113        INQUIRE( FILE = 'PIDS_CHEM' // TRIM( coupling_char ), EXIST = netcdf_extend )
    9114        IF ( .NOT. netcdf_extend )  THEN
    9115           message_string = 'Input file PIDS_CHEM' //  TRIM( coupling_char ) // ' missing!'
    9116           CALL message( 'salsa_gas_emission_setup', 'PA0640', 1, 2, 0, 6, 0 )
    9117        ENDIF
    9118 !
    9119 !--    Open file in read-only mode
    9120        CALL open_read_file( 'PIDS_CHEM' // TRIM( coupling_char ), id_chem )
    91219167!
    91229168!--    Read the index and name of chemical components
     
    93019347       DEALLOCATE( chem_emission%default_emission_data )
    93029348    ENDIF
     9349!
     9350!-- Close input file
     9351    CALL close_input_file( id_chem )
     9352
    93039353#else
    93049354    message_string = 'salsa_emission_mode = "read_from_file", but preprocessor directive ' //   &
     
    93479397       TYPE(surf_type), INTENT(inout) :: surface  !< respective surface type
    93489398
     9399       conv = 1.0_wp
    93499400       use_time_fac = PRESENT( time_fac )
    93509401
     
    97659816             ENDDO
    97669817
    9767           CASE ( 's_BC', 's_DU', 's_H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
     9818          CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
    97689819             IF ( is_used( prtcl, TRIM( variable(3:) ) ) )  THEN
    97699820                found_index = get_index( prtcl, TRIM( variable(3:) ) )
     
    97869837                ENDDO
    97879838             ENDIF
     9839
     9840          CASE ( 's_H2O' )
     9841             found_index = get_index( prtcl,'H2O' )
     9842             DO  i = nxlg, nxrg
     9843                DO  j = nysg, nyng
     9844                   DO  k = nzb, nzt+1
     9845                      DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
     9846                         to_be_resorted(k,j,i) = to_be_resorted(k,j,i) +                           &
     9847                                                 aerosol_mass(ic)%conc(k,j,i)
     9848                      ENDDO
     9849                   ENDDO
     9850                ENDDO
     9851             ENDDO
    97889852
    97899853          CASE DEFAULT
     
    98859949             ENDDO
    98869950
    9887           CASE ( 's_BC', 's_DU', 's_H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
     9951          CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
    98889952             IF ( is_used( prtcl, TRIM( variable(3:) ) ) )  THEN
    9889                 found_index = get_index( prtcl, TRIM( variable(3:) ) )
    98909953                IF ( TRIM( variable(3:) ) == 'BC' )   to_be_resorted => s_bc_av
    98919954                IF ( TRIM( variable(3:) ) == 'DU' )   to_be_resorted => s_du_av
     
    99049967                ENDDO
    99059968             ENDIF
     9969
     9970          CASE ( 's_H2O' )
     9971             to_be_resorted => s_h2o_av
     9972             DO  i = nxlg, nxrg
     9973                DO  j = nysg, nyng
     9974                   DO  k = nzb, nzt+1
     9975                      to_be_resorted(k,j,i) = to_be_resorted(k,j,i) /                              &
     9976                                              REAL( average_count_3d, KIND=wp )
     9977                   ENDDO
     9978                ENDDO
     9979             ENDDO
    99069980
    99079981       END SELECT
     
    1021210286          IF ( mode == 'xy' )  grid = 'zu'
    1021310287
    10214        CASE ( 's_BC', 's_DU', 's_H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
     10288       CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
    1021510289          vari = TRIM( variable( 3:LEN( TRIM( variable ) ) - 3 ) )
    1021610290          IF ( is_used( prtcl, vari ) )  THEN
     
    1025210326          IF ( mode == 'xy' )  grid = 'zu'
    1025310327
     10328       CASE ( 's_H2O' )
     10329          found_index = get_index( prtcl, 'H2O' )
     10330          IF ( av == 0 )  THEN
     10331             DO  i = nxl, nxr
     10332                DO  j = nys, nyn
     10333                   DO  k = nzb_do, nzt_do
     10334                      temp_bin = 0.0_wp
     10335                      DO  ic = ( found_index-1 ) * nbins_aerosol+1, found_index * nbins_aerosol
     10336                         temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     10337                      ENDDO
     10338                      local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     10339                                               BTEST( wall_flags_0(k,j,i), 0 ) )
     10340                   ENDDO
     10341                ENDDO
     10342             ENDDO
     10343          ELSE
     10344             to_be_resorted => s_h2o_av
     10345             DO  i = nxl, nxr
     10346                DO  j = nys, nyn
     10347                   DO  k = nzb_do, nzt_do
     10348                      local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
     10349                                               KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
     10350                   ENDDO
     10351                ENDDO
     10352             ENDDO
     10353          ENDIF
     10354
     10355          IF ( mode == 'xy' )  grid = 'zu'
     10356
    1025410357       CASE DEFAULT
    1025510358          found = .FALSE.
     
    1053410637          ENDIF
    1053510638
    10536        CASE ( 's_BC', 's_DU', 's_H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
     10639       CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
    1053710640          IF ( is_used( prtcl, TRIM( variable(3:) ) ) )  THEN
    1053810641             found_index = get_index( prtcl, TRIM( variable(3:) ) )
     
    1056710670                ENDDO
    1056810671             ENDIF
     10672          ENDIF
     10673
     10674       CASE ( 's_H2O' )
     10675          found_index = get_index( prtcl, 'H2O' )
     10676          IF ( av == 0 )  THEN
     10677             DO  i = nxl, nxr
     10678                DO  j = nys, nyn
     10679                   DO  k = nzb_do, nzt_do
     10680                      temp_bin = 0.0_wp
     10681                      DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
     10682                         temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     10683                      ENDDO
     10684                      local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),            &
     10685                                               BTEST( wall_flags_0(k,j,i), 0 ) )
     10686                   ENDDO
     10687                ENDDO
     10688             ENDDO
     10689          ELSE
     10690             to_be_resorted => s_h2o_av
     10691             DO  i = nxl, nxr
     10692                DO  j = nys, nyn
     10693                   DO  k = nzb_do, nzt_do
     10694                      local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,            &
     10695                                               KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
     10696                   ENDDO
     10697                ENDDO
     10698             ENDDO
    1056910699          ENDIF
    1057010700
     
    1096211092          ENDIF
    1096311093
     11094       CASE ( 's_H2O' )
     11095          IF ( av == 0 )  THEN
     11096             found_index = get_index( prtcl, 'H2O' )
     11097             DO  i = nxl, nxr
     11098                DO  j = nys, nyn
     11099                   DO  k = nzb, nz_do3d
     11100                      temp_bin = 0.0_wp
     11101                      DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
     11102                         temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     11103                      ENDDO
     11104                      tend(k,j,i) = temp_bin
     11105                   ENDDO
     11106                ENDDO
     11107             ENDDO
     11108             IF ( .NOT. mask_surface(mid) )  THEN
     11109                DO  i = 1, mask_size_l(mid,1)
     11110                   DO  j = 1, mask_size_l(mid,2)
     11111                      DO  k = 1, mask_size_l(mid,3)
     11112                         local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) )
     11113                      ENDDO
     11114                   ENDDO
     11115                ENDDO
     11116             ELSE
     11117                DO  i = 1, mask_size_l(mid,1)
     11118                   DO  j = 1, mask_size_l(mid,2)
     11119                      topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i),    &
     11120                                                                  grid )
     11121                      DO  k = 1, mask_size_l(mid,3)
     11122                         local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),         &
     11123                                                 mask_j(mid,j), mask_i(mid,i) )
     11124                      ENDDO
     11125                   ENDDO
     11126                ENDDO
     11127             ENDIF
     11128             resorted = .TRUE.
     11129          ELSE
     11130             to_be_resorted => s_h2o_av
     11131          ENDIF
     11132
    1096411133       CASE DEFAULT
    1096511134          found = .FALSE.
     
    1099911168 END SUBROUTINE salsa_data_output_mask
    1100011169
     11170!------------------------------------------------------------------------------!
     11171! Description:
     11172! ------------
     11173!> Creates index tables for different (aerosol) components
     11174!------------------------------------------------------------------------------!
     11175 SUBROUTINE component_index_constructor( self, ncomp, nlist, listcomp )
     11176
     11177    IMPLICIT NONE
     11178
     11179    INTEGER(iwp) ::  ii  !<
     11180    INTEGER(iwp) ::  jj  !<
     11181
     11182    INTEGER(iwp), INTENT(in) ::  nlist ! < Maximum number of components
     11183
     11184    INTEGER(iwp), INTENT(inout) ::  ncomp  !< Number of components
     11185
     11186    CHARACTER(len=3), INTENT(in) ::  listcomp(nlist)  !< List cof component names
     11187
     11188    TYPE(component_index), INTENT(inout) ::  self  !< Object containing the indices of different
     11189                                                   !< aerosol components
     11190
     11191    ncomp = 0
     11192
     11193    DO WHILE ( listcomp(ncomp+1) /= '  ' .AND. ncomp < nlist )
     11194       ncomp = ncomp + 1
     11195    ENDDO
     11196
     11197    self%ncomp = ncomp
     11198    ALLOCATE( self%ind(ncomp), self%comp(ncomp) )
     11199
     11200    DO  ii = 1, ncomp
     11201       self%ind(ii) = ii
     11202    ENDDO
     11203
     11204    jj = 1
     11205    DO  ii = 1, nlist
     11206       IF ( listcomp(ii) == '') CYCLE
     11207       self%comp(jj) = listcomp(ii)
     11208       jj = jj + 1
     11209    ENDDO
     11210
     11211 END SUBROUTINE component_index_constructor
     11212
     11213!------------------------------------------------------------------------------!
     11214! Description:
     11215! ------------
     11216!> Gives the index of a component in the component list
     11217!------------------------------------------------------------------------------!
     11218 INTEGER FUNCTION get_index( self, incomp )
     11219
     11220    IMPLICIT NONE
     11221
     11222    CHARACTER(len=*), INTENT(in) ::  incomp !< Component name
     11223
     11224    INTEGER(iwp) ::  ii  !< index
     11225
     11226    TYPE(component_index), INTENT(in) ::  self  !< Object containing the indices of different
     11227                                                !< aerosol components
     11228    IF ( ANY( self%comp == incomp ) )  THEN
     11229       ii = 1
     11230       DO WHILE ( (self%comp(ii) /= incomp) )
     11231          ii = ii + 1
     11232       ENDDO
     11233       get_index = ii
     11234    ELSEIF ( incomp == 'H2O' )  THEN
     11235       get_index = self%ncomp + 1
     11236    ELSE
     11237       WRITE( message_string, * ) 'Incorrect component name given!'
     11238       CALL message( 'get_index', 'PA0591', 1, 2, 0, 6, 0 )
     11239    ENDIF
     11240
     11241 END FUNCTION get_index
     11242
     11243!------------------------------------------------------------------------------!
     11244! Description:
     11245! ------------
     11246!> Tells if the (aerosol) component is being used in the simulation
     11247!------------------------------------------------------------------------------!
     11248 LOGICAL FUNCTION is_used( self, icomp )
     11249
     11250    IMPLICIT NONE
     11251
     11252    CHARACTER(len=*), INTENT(in) ::  icomp !< Component name
     11253
     11254    TYPE(component_index), INTENT(in) ::  self  !< Object containing the indices of different
     11255                                                !< aerosol components
     11256
     11257    IF ( ANY(self%comp == icomp) ) THEN
     11258       is_used = .TRUE.
     11259    ELSE
     11260       is_used = .FALSE.
     11261    ENDIF
     11262
     11263 END FUNCTION
     11264
    1100111265 END MODULE salsa_mod
Note: See TracChangeset for help on using the changeset viewer.