Changeset 3956 for palm/trunk/SOURCE


Ignore:
Timestamp:
May 7, 2019 12:32:52 PM (5 years ago)
Author:
monakurppa
Message:

Remove salsa calls from prognostic_equations and correct a bug in the salsa deposition for urban and land surface models

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

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

    r3931 r3956  
    2525! -----------------
    2626! $Id$
     27! - Added calls for salsa_non_advective_processes and
     28!   salsa_exchange_horiz_bounds
     29! - Moved the call for salsa_data_output_2d/3d before that of
     30!   radiation_data_output_2d/3d. radiation_data_output_2d/3d tries to read a
     31!   salsa output variable and encounters a segmentation fault for "Ntot" due
     32!   to the shortoutput name
     33!
     34! 3931 2019-04-24 16:34:28Z schwenkel
    2735! Changed non_transport_physics to non_advective_processes
    2836!
     
    158166               bcm_actions,                                                    &
    159167               bcm_non_advective_processes,                                    &
    160                bcm_exchange_horiz,                                             &               
     168               bcm_exchange_horiz,                                             &
    161169               bcm_prognostic_equations,                                       &
    162170               bcm_swap_timelevel,                                             &
     
    303311               salsa_header,                                                   &
    304312               salsa_actions,                                                  &
     313               salsa_non_advective_processes,                                  &
     314               salsa_exchange_horiz_bounds,                                    &
    305315               salsa_prognostic_equations,                                     &
    306316               salsa_swap_timelevel,                                           &
     
    963973    IF ( bulk_cloud_model    )  CALL bcm_non_advective_processes()
    964974    IF ( air_chemistry       )  CALL chem_non_advective_processes()
     975    IF ( salsa               )  CALL salsa_non_advective_processes()
    965976
    966977
     
    980991
    981992
    982     IF ( bulk_cloud_model    )  CALL bcm_non_advective_processes( i, j )   
     993    IF ( bulk_cloud_model    )  CALL bcm_non_advective_processes( i, j )
    983994    IF ( air_chemistry       )  CALL chem_non_advective_processes( i, j )
     995    IF ( salsa               )  CALL salsa_non_advective_processes( i, j )
    984996
    985997
     
    9961008    IF ( bulk_cloud_model    )  CALL bcm_exchange_horiz()
    9971009    IF ( air_chemistry       )  CALL chem_exchange_horiz_bounds()
     1010    IF ( salsa               )  CALL salsa_exchange_horiz_bounds()
    9981011
    9991012 END SUBROUTINE module_interface_exchange_horiz
     
    11491162    ENDIF
    11501163
     1164    IF ( .NOT. found  .AND.  salsa )  THEN
     1165       CALL salsa_data_output_2d(                                              &
     1166               av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do&
     1167            )
     1168    ENDIF
     1169
    11511170    IF ( .NOT. found  .AND.  radiation )  THEN
    11521171       CALL radiation_data_output_2d(                                          &
    1153                av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do&
    1154             )
    1155     ENDIF
    1156 
    1157     IF ( .NOT. found  .AND.  salsa )  THEN
    1158        CALL salsa_data_output_2d(                                              &
    11591172               av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do&
    11601173            )
     
    12261239    ENDIF
    12271240
     1241    IF ( .NOT. found  .AND.  salsa )  THEN
     1242       CALL salsa_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
     1243       resorted = .TRUE.
     1244    ENDIF
     1245
    12281246    IF ( .NOT. found  .AND.  radiation )  THEN
    12291247       CALL radiation_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
    1230        resorted = .TRUE.
    1231     ENDIF
    1232 
    1233     IF ( .NOT. found  .AND.  salsa )  THEN
    1234        CALL salsa_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
    12351248       resorted = .TRUE.
    12361249    ENDIF
  • palm/trunk/SOURCE/prognostic_equations.f90

    r3931 r3956  
    2525! -----------------
    2626! $Id$
     27! Removed salsa calls.
     28!
     29! 3931 2019-04-24 16:34:28Z schwenkel
    2730! Correct/complete module_interface introduction for chemistry model
    2831!
     
    464467#endif
    465468
    466     USE salsa_mod,                                                             &
    467         ONLY:  aerosol_mass, aerosol_number, dt_salsa, last_salsa_time,        &
    468                nbins_aerosol, ncomponents_mass, ngases_salsa,                  &
    469                salsa_boundary_conds, salsa_diagnostics, salsa_driver,          &
    470                salsa_gas, salsa_gases_from_chem, skip_time_do_salsa
    471 
    472469    USE statistics,                                                            &
    473470        ONLY:  hom
     
    518515    INTEGER(iwp) ::  i                   !<
    519516    INTEGER(iwp) ::  i_omp_start         !<
    520     INTEGER(iwp) ::  ib                  !< index for aerosol size bins (salsa)
    521     INTEGER(iwp) ::  ic                  !< index for chemical compounds (salsa)
    522     INTEGER(iwp) ::  icc                 !< additional index for chemical compounds (salsa)
    523     INTEGER(iwp) ::  ig                  !< index for gaseous compounds (salsa)
    524517    INTEGER(iwp) ::  j                   !<
    525518    INTEGER(iwp) ::  k                   !<
     
    550543!
    551544!-- Module Inferface for exchange horiz after non_advective_processes but before
    552 !-- advection. Therefore, non_advective_processes must not run for ghost points.
     545!-- advection. Therefore, non_advective_processes must not run for ghost points.
     546    !$OMP END PARALLEL
    553547    CALL module_interface_exchange_horiz()
    554     !$OMP END PARALLEL
    555 
    556 !
    557 !-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
    558 !-- step. The exchange of ghost points is required after this update of the
    559 !-- concentrations of aerosol number and mass
    560     IF ( salsa )  THEN
    561 
    562        IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
    563           IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
    564           THEN
    565              CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
    566              !$OMP PARALLEL PRIVATE (i,j)
    567 !
    568 !--          Call salsa processes
    569              !$OMP DO
    570              DO  i = nxl, nxr
    571                 DO  j = nys, nyn
    572                    CALL salsa_diagnostics( i, j )
    573                    CALL salsa_driver( i, j, 3 )
    574                    CALL salsa_diagnostics( i, j )
    575                 ENDDO
    576              ENDDO
    577              !$OMP END PARALLEL
    578 
    579              CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
    580 
    581              CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
    582 !
    583 !--          Exchange ghost points and decycle if needed.
    584              DO  ib = 1, nbins_aerosol
    585                 CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
    586                 CALL salsa_boundary_conds( aerosol_number(ib)%conc,            &
    587                                            aerosol_number(ib)%init )
    588                 DO  ic = 1, ncomponents_mass
    589                    icc = ( ic - 1 ) * nbins_aerosol + ib
    590                    CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
    591                    CALL salsa_boundary_conds( aerosol_mass(icc)%conc,          &
    592                                               aerosol_mass(icc)%init )
    593                 ENDDO
    594              ENDDO
    595 
    596              IF ( .NOT. salsa_gases_from_chem )  THEN
    597                 DO  ig = 1, ngases_salsa
    598                    CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
    599                    CALL salsa_boundary_conds( salsa_gas(ig)%conc,              &
    600                                               salsa_gas(ig)%init )
    601                 ENDDO
    602              ENDIF
    603              CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
    604              last_salsa_time = time_since_reference_point
    605 
    606           ENDIF
    607 
    608        ENDIF
    609 
    610     ENDIF
    611 
    612548!
    613549!-- Loop over all prognostic equations
     
    11531089
    11541090    INTEGER(iwp) ::  i     !<
    1155     INTEGER(iwp) ::  ib    !< index for aerosol size bins (salsa)
    1156     INTEGER(iwp) ::  ic    !< index for chemical compounds (salsa)
    1157     INTEGER(iwp) ::  icc   !< additional index for chemical compounds (salsa)
    1158     INTEGER(iwp) ::  ig    !< index for gaseous compounds (salsa)
    11591091    INTEGER(iwp) ::  j     !<
    11601092    INTEGER(iwp) ::  k     !<
     
    11681100    ENDIF
    11691101!
    1170 !-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
    1171 !-- step. The exchange of ghost points is required after this update of the
    1172 !-- concentrations of aerosol number and mass
    1173     IF ( salsa )  THEN
    1174        IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
    1175           IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
    1176           THEN
    1177              CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
    1178              !$OMP PARALLEL PRIVATE (i,j)
    1179              !$OMP DO
    1180 !
    1181 !--          Call salsa processes
    1182              DO  i = nxl, nxr
    1183                 DO  j = nys, nyn
    1184                    CALL salsa_diagnostics( i, j )
    1185                    CALL salsa_driver( i, j, 3 )
    1186                    CALL salsa_diagnostics( i, j )
    1187                 ENDDO
    1188              ENDDO
    1189              !$OMP END PARALLEL
    1190 
    1191              CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
    1192              CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
    1193 !
    1194 !--          Exchange ghost points and decycle if needed.
    1195              DO  ib = 1, nbins_aerosol
    1196                 CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
    1197                 CALL salsa_boundary_conds( aerosol_number(ib)%conc,            &
    1198                                            aerosol_number(ib)%init )
    1199                 DO  ic = 1, ncomponents_mass
    1200                    icc = ( ic - 1 ) * nbins_aerosol + ib
    1201                    CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
    1202                    CALL salsa_boundary_conds( aerosol_mass(icc)%conc,          &
    1203                                               aerosol_mass(icc)%init )
    1204                 ENDDO
    1205              ENDDO
    1206              IF ( .NOT. salsa_gases_from_chem )  THEN
    1207                 DO  ig = 1, ngases_salsa
    1208                    CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
    1209                    CALL salsa_boundary_conds( salsa_gas(ig)%conc,              &
    1210                                               salsa_gas(ig)%init )
    1211                 ENDDO
    1212              ENDIF
    1213              CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
    1214              last_salsa_time = time_since_reference_point
    1215           ENDIF
    1216        ENDIF
    1217     ENDIF
    1218 
    1219 !
    12201102!-- Calculate non advective processes for all other modules
    12211103    CALL module_interface_non_advective_processes
    12221104!
    12231105!-- Module Inferface for exchange horiz after non_advective_processes but before
    1224 !-- advection. Therefore, non_advective_processes must not run for ghost points.     
     1106!-- advection. Therefore, non_advective_processes must not run for ghost points.
    12251107    CALL module_interface_exchange_horiz()
    1226    
    12271108!
    12281109!-- u-velocity component
  • palm/trunk/SOURCE/salsa_mod.f90

    r3924 r3956  
    2626! -----------------
    2727! $Id$
     28! - Conceptual bug in depo_surf correct for urban and land surface model
     29! - Subroutine salsa_tendency_ij optimized.
     30! - Interfaces salsa_non_advective_processes and salsa_exchange_horiz_bounds
     31!   created. These are now called in module_interface.
     32!   salsa_exchange_horiz_bounds after calling salsa_driver only when needed
     33!   (i.e. every dt_salsa).
     34!
     35! 3924 2019-04-23 09:33:06Z monakurppa
    2836! Correct a bug introduced by the previous update.
    2937!
     
    173181!-- Local constants:
    174182    INTEGER(iwp), PARAMETER ::  luc_urban = 15     !< default landuse type for urban
    175     INTEGER(iwp), PARAMETER ::  ngases_salsa   = 5 !< total number of gaseous tracers:
     183    INTEGER(iwp), PARAMETER ::  ngases_salsa  = 5 !< total number of gaseous tracers:
    176184                                                   !< 1 = H2SO4, 2 = HNO3, 3 = NH3, 4 = OCNV
    177185                                                   !< (non-volatile OC), 5 = OCSV (semi-volatile)
    178     INTEGER(iwp), PARAMETER ::  nmod    = 7  !< number of modes for initialising the aerosol size
     186    INTEGER(iwp), PARAMETER ::  nmod = 7     !< number of modes for initialising the aerosol size
    179187                                             !< distribution
    180     INTEGER(iwp), PARAMETER ::  nreg    = 2  !< Number of main size subranges
     188    INTEGER(iwp), PARAMETER ::  nreg = 2     !< Number of main size subranges
    181189    INTEGER(iwp), PARAMETER ::  maxspec = 7  !< Max. number of aerosol species
    182190    INTEGER(iwp), PARAMETER ::  season = 1   !< For dry depostion by Zhang et al.: 1 = summer,
     
    208216!
    209217!-- Molar masses in kg/mol
    210     REAL(wp), PARAMETER ::  ambc   = 12.0E-3_wp     !< black carbon (BC)
    211     REAL(wp), PARAMETER ::  amdair = 28.970E-3_wp   !< dry air
    212     REAL(wp), PARAMETER ::  amdu   = 100.E-3_wp     !< mineral dust
    213     REAL(wp), PARAMETER ::  amh2o  = 18.0154E-3_wp  !< H2O
    214     REAL(wp), PARAMETER ::  amh2so4  = 98.06E-3_wp  !< H2SO4
    215     REAL(wp), PARAMETER ::  amhno3 = 63.01E-3_wp    !< HNO3
    216     REAL(wp), PARAMETER ::  amn2o  = 44.013E-3_wp   !< N2O
    217     REAL(wp), PARAMETER ::  amnh3  = 17.031E-3_wp   !< NH3
    218     REAL(wp), PARAMETER ::  amo2   = 31.9988E-3_wp  !< O2
    219     REAL(wp), PARAMETER ::  amo3   = 47.998E-3_wp   !< O3
    220     REAL(wp), PARAMETER ::  amoc   = 150.E-3_wp     !< organic carbon (OC)
    221     REAL(wp), PARAMETER ::  amss   = 58.44E-3_wp    !< sea salt (NaCl)
     218    REAL(wp), PARAMETER ::  ambc     = 12.0E-3_wp     !< black carbon (BC)
     219    REAL(wp), PARAMETER ::  amdair   = 28.970E-3_wp   !< dry air
     220    REAL(wp), PARAMETER ::  amdu     = 100.E-3_wp     !< mineral dust
     221    REAL(wp), PARAMETER ::  amh2o    = 18.0154E-3_wp  !< H2O
     222    REAL(wp), PARAMETER ::  amh2so4  = 98.06E-3_wp    !< H2SO4
     223    REAL(wp), PARAMETER ::  amhno3   = 63.01E-3_wp    !< HNO3
     224    REAL(wp), PARAMETER ::  amn2o    = 44.013E-3_wp   !< N2O
     225    REAL(wp), PARAMETER ::  amnh3    = 17.031E-3_wp   !< NH3
     226    REAL(wp), PARAMETER ::  amo2     = 31.9988E-3_wp  !< O2
     227    REAL(wp), PARAMETER ::  amo3     = 47.998E-3_wp   !< O3
     228    REAL(wp), PARAMETER ::  amoc     = 150.E-3_wp     !< organic carbon (OC)
     229    REAL(wp), PARAMETER ::  amss     = 58.44E-3_wp    !< sea salt (NaCl)
    222230!
    223231!-- Densities in kg/m3
     
    243251!-- C_B, C_IN, C_IM, beta_IM and C_IT for each land use category (15, as in Zhang et al. (2001))
    244252    REAL(wp), DIMENSION(1:15), PARAMETER :: l_p10 = &
    245         (/0.15, 4.0, 0.15, 3.0, 3.0, 0.5, 3.0, -99., 0.5, 2.0, 1.0, -99., -99., -99., 3.0/)
     253        (/0.15,  4.0,   0.15,  3.0,   3.0,   0.5,   3.0,   -99., 0.5,  2.0,  1.0,  -99., -99., -99., 3.0/)
    246254    REAL(wp), DIMENSION(1:15), PARAMETER :: c_b_p10 = &
    247         (/0.887, 1.262, 0.887, 1.262, 1.262, 0.996, 0.996, -99., 0.7, 0.93, 0.996, -99., -99., -99., 1.262/)
     255        (/0.887, 1.262, 0.887, 1.262, 1.262, 0.996, 0.996, -99., 0.7,  0.93, 0.996, -99., -99., -99., 1.262/)
    248256    REAL(wp), DIMENSION(1:15), PARAMETER :: c_in_p10 = &
    249         (/0.81, 0.216, 0.81, 0.216, 0.216, 0.191, 0.162, -99., 0.7, 0.14, 0.162, -99., -99., -99., 0.216/)
     257        (/0.81,  0.216, 0.81,  0.216, 0.216, 0.191, 0.162, -99., 0.7, 0.14, 0.162, -99., -99., -99., 0.216/)
    250258    REAL(wp), DIMENSION(1:15), PARAMETER :: c_im_p10 = &
    251         (/0.162, 0.13, 0.162, 0.13, 0.13, 0.191, 0.081, -99., 0.191, 0.086, 0.081, -99., -99., -99., 0.13/)
     259        (/0.162, 0.13,  0.162, 0.13,  0.13,  0.191, 0.081, -99., 0.191,0.086,0.081, -99., -99., -99., 0.13/)
    252260    REAL(wp), DIMENSION(1:15), PARAMETER :: beta_im_p10 = &
    253         (/0.6, 0.47, 0.6, 0.47, 0.47, 0.47, 0.47, -99., 0.6, 0.47, 0.47, -99., -99., -99., 0.47/)
     261        (/0.6,   0.47,  0.6,   0.47,  0.47,  0.47,  0.47,  -99., 0.6,  0.47, 0.47, -99., -99., -99., 0.47/)
    254262    REAL(wp), DIMENSION(1:15), PARAMETER :: c_it_p10 = &
    255         (/0.0, 0.056, 0.0, 0.056, 0.056, 0.042, 0.056, -99., 0.042, 0.014, 0.056, -99., -99., -99., 0.056/)
     263        (/0.0,   0.056, 0.0,   0.056, 0.056, 0.042, 0.056, -99., 0.042,0.014,0.056, -99., -99., -99., 0.056/)
    256264!
    257265!-- Constants for the dry deposition model by Zhang et al. (2001):
     
    307315                                   (/'SO4','   ','   ','   ','   ','   ','   '/)
    308316
     317    INTEGER(iwp) ::  depo_pcm_par_num = 1   !< parametrisation type: 1=zhang2001, 2=petroff2010
    309318    INTEGER(iwp) ::  depo_pcm_type_num = 0  !< index for the dry deposition type on the plant canopy
     319    INTEGER(iwp) ::  depo_surf_par_num = 1  !< parametrisation type: 1=zhang2001, 2=petroff2010
    310320    INTEGER(iwp) ::  dots_salsa = 0         !< starting index for salsa-timeseries
    311321    INTEGER(iwp) ::  end_subrange_1a = 1    !< last index for bin subrange 1a
     
    364374!
    365375!-- SALSA switches:
    366     LOGICAL ::  advect_particle_water = .TRUE.     !< advect water concentration of particles
    367     LOGICAL ::  decycle_lr            = .FALSE.    !< Undo cyclic boundary conditions: left and right
    368     LOGICAL ::  decycle_ns            = .FALSE.    !< north and south boundaries
    369     LOGICAL ::  include_emission      = .FALSE.    !< include or not emissions
    370     LOGICAL ::  feedback_to_palm      = .FALSE.    !< allow feedback due to condensation of H2O
    371     LOGICAL ::  nest_salsa            = .FALSE.    !< apply nesting for salsa
    372     LOGICAL ::  no_insoluble          = .FALSE.    !< Switch to exclude insoluble chemical components
    373     LOGICAL ::  read_restart_data_salsa = .FALSE.  !< read restart data for salsa
    374     LOGICAL ::  salsa_gases_from_chem = .FALSE.    !< Transfer the gaseous components to SALSA from
    375                                                    !< from chemistry model
    376     LOGICAL ::  van_der_waals_coagc   = .FALSE.    !< Enhancement of coagulation kernel by van der
     376    LOGICAL ::  advect_particle_water   = .TRUE.   !< Advect water concentration of particles
     377    LOGICAL ::  decycle_lr              = .FALSE.  !< Undo cyclic boundaries: left and right
     378    LOGICAL ::  decycle_ns              = .FALSE.  !< Undo cyclic boundaries: north and south
     379    LOGICAL ::  include_emission        = .FALSE.  !< Include or not emissions
     380    LOGICAL ::  feedback_to_palm        = .FALSE.  !< Allow feedback due to condensation of H2O
     381    LOGICAL ::  nest_salsa              = .FALSE.  !< Apply nesting for salsa
     382    LOGICAL ::  no_insoluble            = .FALSE.  !< Exclude insoluble chemical components
     383    LOGICAL ::  read_restart_data_salsa = .FALSE.  !< Read restart data for salsa
     384    LOGICAL ::  salsa_gases_from_chem   = .FALSE.  !< Transfer the gaseous components to SALSA from
     385                                                   !< the chemistry model
     386    LOGICAL ::  van_der_waals_coagc     = .FALSE.  !< Enhancement of coagulation kernel by van der
    377387                                                   !< Waals and viscous forces
    378     LOGICAL ::  write_binary_salsa    = .FALSE.    !< read binary for salsa
     388    LOGICAL ::  write_binary_salsa      = .FALSE.  !< read binary for salsa
    379389!
    380390!-- Process switches: nl* is read from the NAMELIST and is NOT changed.
     
    455465!-- SALSA derived datatypes:
    456466!
    457 !-- For matching LSM and the deposition module surface types
    458     TYPE match_lsm_depo
    459        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match
    460     END TYPE match_lsm_depo
     467!-- For matching LSM and USM surface types and the deposition module surface types
     468    TYPE match_surface
     469       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_lupg  !< index for pavement / green roofs
     470       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_luvw  !< index for vegetation / walls
     471       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_luww  !< index for water / windows
     472    END TYPE match_surface
    461473!
    462474!-- Aerosol emission data attributes
     
    598610    TYPE(t_section), DIMENSION(:), ALLOCATABLE ::  aero  !< local aerosol properties
    599611
    600     TYPE(match_lsm_depo) ::  lsm_to_depo_h
    601 
    602     TYPE(match_lsm_depo), DIMENSION(0:3) ::  lsm_to_depo_v
     612    TYPE(match_surface) ::  lsm_to_depo_h  !< to match the deposition module and horizontal LSM surfaces
     613    TYPE(match_surface) ::  usm_to_depo_h  !< to match the deposition module and horizontal USM surfaces
     614
     615    TYPE(match_surface), DIMENSION(0:3) ::  lsm_to_depo_v  !< to match the deposition mod. and vertical LSM surfaces
     616    TYPE(match_surface), DIMENSION(0:3) ::  usm_to_depo_v  !< to match the deposition mod. and vertical USM surfaces
    603617!
    604618!-- SALSA variables: as x = x(k,j,i,bin).
     
    752766    END INTERFACE salsa_actions
    753767!
     768!-- Non-advective processes (i.e. aerosol dynamic reactions) for salsa variables
     769    INTERFACE salsa_non_advective_processes
     770       MODULE PROCEDURE salsa_non_advective_processes
     771       MODULE PROCEDURE salsa_non_advective_processes_ij
     772    END INTERFACE salsa_non_advective_processes
     773!
     774!-- Exchange horiz for salsa variables
     775    INTERFACE salsa_exchange_horiz_bounds
     776       MODULE PROCEDURE salsa_exchange_horiz_bounds
     777    END INTERFACE salsa_exchange_horiz_bounds
     778!
    754779!-- Prognostics equations for salsa variables
    755780    INTERFACE salsa_prognostic_equations
     
    775800           salsa_emission_update, salsa_header, salsa_init, salsa_init_arrays, salsa_parin,        &
    776801           salsa_rrd_local, salsa_swap_timelevel, salsa_prognostic_equations, salsa_wrd_local,     &
    777            salsa_actions
     802           salsa_actions, salsa_non_advective_processes, salsa_exchange_horiz_bounds
    778803!
    779804!-- Public parameters, constants and initial values
     
    21622187
    21632188    USE surface_mod,                                                                               &
    2164         ONLY:  surf_lsm_h, surf_lsm_v
     2189        ONLY:  surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
    21652190
    21662191    IMPLICIT NONE
     
    21682193    INTEGER(iwp) ::  l  !< loop index for vertical surfaces
    21692194
     2195    LOGICAL :: match_lsm  !< flag to initilise LSM surfaces (if false, initialise USM surfaces)
     2196
     2197    IF ( depo_pcm_par == 'zhang2001' )  THEN
     2198       depo_pcm_par_num = 1
     2199    ELSEIF ( depo_pcm_par == 'petroff2010' )  THEN
     2200       depo_pcm_par_num = 2
     2201    ENDIF
     2202
     2203    IF ( depo_surf_par == 'zhang2001' )  THEN
     2204       depo_surf_par_num = 1
     2205    ELSEIF ( depo_surf_par == 'petroff2010' )  THEN
     2206       depo_surf_par_num = 2
     2207    ENDIF
     2208!
     2209!-- LSM: Pavement, vegetation and water
    21702210    IF ( nldepo_surf  .AND.  land_surface )  THEN
    2171 
    2172        ALLOCATE( lsm_to_depo_h%match(1:surf_lsm_h%ns) )
    2173        lsm_to_depo_h%match = 0
    2174        CALL match_lsm_zhang( surf_lsm_h, lsm_to_depo_h%match )
    2175 
     2211       match_lsm = .TRUE.
     2212       ALLOCATE( lsm_to_depo_h%match_lupg(1:surf_lsm_h%ns),                                         &
     2213                 lsm_to_depo_h%match_luvw(1:surf_lsm_h%ns),                                         &
     2214                 lsm_to_depo_h%match_luww(1:surf_lsm_h%ns) )
     2215       lsm_to_depo_h%match_lupg = 0
     2216       lsm_to_depo_h%match_luvw = 0
     2217       lsm_to_depo_h%match_luww = 0
     2218       CALL match_sm_zhang( surf_lsm_h, lsm_to_depo_h%match_lupg, lsm_to_depo_h%match_luvw,        &
     2219                            lsm_to_depo_h%match_luww, match_lsm )
    21762220       DO  l = 0, 3
    2177           ALLOCATE( lsm_to_depo_v(l)%match(1:surf_lsm_v(l)%ns) )
    2178           lsm_to_depo_v(l)%match = 0
    2179           CALL match_lsm_zhang( surf_lsm_v(l), lsm_to_depo_v(l)%match )
     2221          ALLOCATE( lsm_to_depo_v(l)%match_lupg(1:surf_lsm_v(l)%ns),                               &
     2222                    lsm_to_depo_v(l)%match_luvw(1:surf_lsm_v(l)%ns),                               &
     2223                    lsm_to_depo_v(l)%match_luww(1:surf_lsm_v(l)%ns) )
     2224          lsm_to_depo_v(l)%match_lupg = 0
     2225          lsm_to_depo_v(l)%match_luvw = 0
     2226          lsm_to_depo_v(l)%match_luww = 0
     2227          CALL match_sm_zhang( surf_lsm_v(l), lsm_to_depo_v(l)%match_lupg,                         &
     2228                               lsm_to_depo_v(l)%match_luvw, lsm_to_depo_v(l)%match_luww, match_lsm )
     2229       ENDDO
     2230    ENDIF
     2231!
     2232!-- USM: Green roofs/walls, wall surfaces and windows
     2233    IF ( nldepo_surf  .AND.  urban_surface )  THEN
     2234       match_lsm = .FALSE.
     2235       ALLOCATE( usm_to_depo_h%match_lupg(1:surf_usm_h%ns),                                        &
     2236                 usm_to_depo_h%match_luvw(1:surf_usm_h%ns),                                        &
     2237                 usm_to_depo_h%match_luww(1:surf_usm_h%ns) )
     2238       usm_to_depo_h%match_lupg = 0
     2239       usm_to_depo_h%match_luvw = 0
     2240       usm_to_depo_h%match_luww = 0
     2241       CALL match_sm_zhang( surf_usm_h, usm_to_depo_h%match_lupg, usm_to_depo_h%match_luvw,        &
     2242                            usm_to_depo_h%match_luww, match_lsm )
     2243       DO  l = 0, 3
     2244          ALLOCATE( usm_to_depo_v(l)%match_lupg(1:surf_usm_v(l)%ns),                               &
     2245                    usm_to_depo_v(l)%match_luvw(1:surf_usm_v(l)%ns),                               &
     2246                    usm_to_depo_v(l)%match_luww(1:surf_usm_v(l)%ns) )
     2247          usm_to_depo_v(l)%match_lupg = 0
     2248          usm_to_depo_v(l)%match_luvw = 0
     2249          usm_to_depo_v(l)%match_luww = 0
     2250          CALL match_sm_zhang( surf_usm_v(l), usm_to_depo_v(l)%match_lupg,                         &
     2251                               usm_to_depo_v(l)%match_luvw, usm_to_depo_v(l)%match_luww, match_lsm )
    21802252       ENDDO
    21812253    ENDIF
     
    22042276!> Match the surface types in PALM and Zhang et al. 2001 deposition module
    22052277!------------------------------------------------------------------------------!
    2206  SUBROUTINE match_lsm_zhang( surf, match_array )
     2278 SUBROUTINE match_sm_zhang( surf, match_pav_green, match_veg_wall, match_wat_win, match_lsm )
    22072279
    22082280    USE surface_mod,                                                           &
     
    22112283    IMPLICIT NONE
    22122284
    2213     INTEGER(iwp) ::  m                !< index for surface elements
    2214     INTEGER(iwp) ::  pav_type_palm    !< pavement type in PALM
    2215     INTEGER(iwp) ::  vege_type_palm   !< vegetation type in PALM
    2216     INTEGER(iwp) ::  water_type_palm  !< water type in PALM
    2217 
    2218     INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_array !< array matching
    2219                                                               !< the surface types
     2285    INTEGER(iwp) ::  m              !< index for surface elements
     2286    INTEGER(iwp) ::  pav_type_palm  !< pavement / green wall type in PALM
     2287    INTEGER(iwp) ::  veg_type_palm  !< vegetation / wall type in PALM
     2288    INTEGER(iwp) ::  wat_type_palm  !< water / window type in PALM
     2289
     2290    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_pav_green  !<  matching pavement/green walls
     2291    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_veg_wall   !<  matching vegetation/walls
     2292    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_wat_win    !<  matching water/windows
     2293
     2294    LOGICAL, INTENT(in) :: match_lsm  !< flag to initilise LSM surfaces (if false, initialise USM)
     2295
    22202296    TYPE(surf_type), INTENT(in) :: surf  !< respective surface type
    22212297
    22222298    DO  m = 1, surf%ns
    2223 
    2224        IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
    2225           vege_type_palm = surf%vegetation_type(m)
    2226           SELECT CASE ( vege_type_palm )
    2227              CASE ( 0 )
    2228                 message_string = 'No vegetation type defined.'
    2229                 CALL message( 'salsa_mod: init_depo_surfaces', 'PA0614', 1, 2, 0, 6, 0 )
    2230              CASE ( 1 )  ! bare soil
    2231                 match_array(m) = 6  ! grass in Z01
    2232              CASE ( 2 )  ! crops, mixed farming
    2233                 match_array(m) = 7  !  crops, mixed farming Z01
    2234              CASE ( 3 )  ! short grass
    2235                 match_array(m) = 6  ! grass in Z01
    2236              CASE ( 4 )  ! evergreen needleleaf trees
    2237                  match_array(m) = 1  ! evergreen needleleaf trees in Z01
    2238              CASE ( 5 )  ! deciduous needleleaf trees
    2239                 match_array(m) = 3  ! deciduous needleleaf trees in Z01
    2240              CASE ( 6 )  ! evergreen broadleaf trees
    2241                 match_array(m) = 2  ! evergreen broadleaf trees in Z01
    2242              CASE ( 7 )  ! deciduous broadleaf trees
    2243                 match_array(m) = 4  ! deciduous broadleaf trees in Z01
    2244              CASE ( 8 )  ! tall grass
    2245                 match_array(m) = 6  ! grass in Z01
    2246              CASE ( 9 )  ! desert
    2247                 match_array(m) = 8  ! desert in Z01
    2248              CASE ( 10 )  ! tundra
    2249                 match_array(m) = 9  ! tundra in Z01
    2250              CASE ( 11 )  ! irrigated crops
    2251                 match_array(m) = 7  !  crops, mixed farming Z01
    2252              CASE ( 12 )  ! semidesert
    2253                 match_array(m) = 8  ! desert in Z01
    2254              CASE ( 13 )  ! ice caps and glaciers
    2255                 match_array(m) = 12  ! ice cap and glacier in Z01
    2256              CASE ( 14 )  ! bogs and marshes
    2257                 match_array(m) = 11  ! wetland with plants in Z01
    2258              CASE ( 15 )  ! evergreen shrubs
    2259                 match_array(m) = 10  ! shrubs and interrupted woodlands in Z01
    2260              CASE ( 16 )  ! deciduous shrubs
    2261                 match_array(m) = 10  ! shrubs and interrupted woodlands in Z01
    2262              CASE ( 17 )  ! mixed forest/woodland
    2263                 match_array(m) = 5  ! mixed broadleaf and needleleaf trees in Z01
    2264              CASE ( 18 )  ! interrupted forest
    2265                 match_array(m) = 10  ! shrubs and interrupted woodlands in Z01
    2266           END SELECT
    2267        ENDIF
    2268 
    2269        IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
    2270           pav_type_palm = surf%pavement_type(m)
    2271           IF ( pav_type_palm == 0 )  THEN  ! error
    2272              message_string = 'No pavement type defined.'
    2273              CALL message( 'salsa_mod: match_lsm_zhang', 'PA0615', 1, 2, 0, 6, 0 )
    2274           ELSEIF ( pav_type_palm > 0  .AND.  pav_type_palm <= 15 )  THEN
    2275              match_array(m) = 15  ! urban in Z01
     2299       IF ( match_lsm )  THEN
     2300!
     2301!--       Vegetation (LSM):
     2302          IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
     2303             veg_type_palm = surf%vegetation_type(m)
     2304             SELECT CASE ( veg_type_palm )
     2305                CASE ( 0 )
     2306                   message_string = 'No vegetation type defined.'
     2307                   CALL message( 'salsa_mod: init_depo_surfaces', 'PA0614', 1, 2, 0, 6, 0 )
     2308                CASE ( 1 )  ! bare soil
     2309                   match_veg_wall(m) = 6  ! grass in Z01
     2310                CASE ( 2 )  ! crops, mixed farming
     2311                   match_veg_wall(m) = 7  !  crops, mixed farming Z01
     2312                CASE ( 3 )  ! short grass
     2313                   match_veg_wall(m) = 6  ! grass in Z01
     2314                CASE ( 4 )  ! evergreen needleleaf trees
     2315                    match_veg_wall(m) = 1  ! evergreen needleleaf trees in Z01
     2316                CASE ( 5 )  ! deciduous needleleaf trees
     2317                   match_veg_wall(m) = 3  ! deciduous needleleaf trees in Z01
     2318                CASE ( 6 )  ! evergreen broadleaf trees
     2319                   match_veg_wall(m) = 2  ! evergreen broadleaf trees in Z01
     2320                CASE ( 7 )  ! deciduous broadleaf trees
     2321                   match_veg_wall(m) = 4  ! deciduous broadleaf trees in Z01
     2322                CASE ( 8 )  ! tall grass
     2323                   match_veg_wall(m) = 6  ! grass in Z01
     2324                CASE ( 9 )  ! desert
     2325                   match_veg_wall(m) = 8  ! desert in Z01
     2326                CASE ( 10 )  ! tundra
     2327                   match_veg_wall(m) = 9  ! tundra in Z01
     2328                CASE ( 11 )  ! irrigated crops
     2329                   match_veg_wall(m) = 7  !  crops, mixed farming Z01
     2330                CASE ( 12 )  ! semidesert
     2331                   match_veg_wall(m) = 8  ! desert in Z01
     2332                CASE ( 13 )  ! ice caps and glaciers
     2333                   match_veg_wall(m) = 12  ! ice cap and glacier in Z01
     2334                CASE ( 14 )  ! bogs and marshes
     2335                   match_veg_wall(m) = 11  ! wetland with plants in Z01
     2336                CASE ( 15 )  ! evergreen shrubs
     2337                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
     2338                CASE ( 16 )  ! deciduous shrubs
     2339                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
     2340                CASE ( 17 )  ! mixed forest/woodland
     2341                   match_veg_wall(m) = 5  ! mixed broadleaf and needleleaf trees in Z01
     2342                CASE ( 18 )  ! interrupted forest
     2343                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
     2344             END SELECT
    22762345          ENDIF
    2277        ENDIF
    2278 
    2279        IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
    2280           water_type_palm = surf%water_type(m)
    2281           IF ( water_type_palm == 0 )  THEN  ! error
    2282              message_string = 'No water type defined.'
    2283              CALL message( 'salsa_mod: match_lsm_zhang', 'PA0616', 1, 2, 0, 6, 0 )
    2284           ELSEIF ( water_type_palm == 3 )  THEN
    2285              match_array(m) = 14  ! ocean in Z01
    2286           ELSEIF ( water_type_palm == 1  .OR.  water_type_palm == 2 .OR.  water_type_palm == 4     &
    2287                    .OR.  water_type_palm == 5  )  THEN
    2288              match_array(m) = 13  ! inland water in Z01
     2346!
     2347!--       Pavement (LSM):
     2348          IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
     2349             pav_type_palm = surf%pavement_type(m)
     2350             IF ( pav_type_palm == 0 )  THEN  ! error
     2351                message_string = 'No pavement type defined.'
     2352                CALL message( 'salsa_mod: match_sm_zhang', 'PA0615', 1, 2, 0, 6, 0 )
     2353             ELSE
     2354                match_pav_green(m) = 15  ! urban in Z01
     2355             ENDIF
    22892356          ENDIF
     2357!
     2358!--       Water (LSM):
     2359          IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
     2360             wat_type_palm = surf%water_type(m)
     2361             IF ( wat_type_palm == 0 )  THEN  ! error
     2362                message_string = 'No water type defined.'
     2363                CALL message( 'salsa_mod: match_sm_zhang', 'PA0616', 1, 2, 0, 6, 0 )
     2364             ELSEIF ( wat_type_palm == 3 )  THEN
     2365                match_wat_win(m) = 14  ! ocean in Z01
     2366             ELSEIF ( wat_type_palm == 1  .OR.  wat_type_palm == 2 .OR.  wat_type_palm == 4        &
     2367                      .OR.  wat_type_palm == 5  )  THEN
     2368                match_wat_win(m) = 13  ! inland water in Z01
     2369             ENDIF
     2370          ENDIF
     2371       ELSE
     2372!
     2373!--       Wall surfaces (USM):
     2374          IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
     2375             match_veg_wall(m) = 15  ! urban in Z01
     2376          ENDIF
     2377!
     2378!--       Green walls and roofs (USM):
     2379          IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
     2380             match_pav_green(m) =  6 ! (short) grass in Z01
     2381          ENDIF
     2382!
     2383!--       Windows (USM):
     2384          IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
     2385             match_wat_win(m) = 15  ! urban in Z01
     2386          ENDIF
    22902387       ENDIF
    22912388
    22922389    ENDDO
    22932390
    2294  END SUBROUTINE match_lsm_zhang
     2391 END SUBROUTINE match_sm_zhang
    22952392
    22962393!------------------------------------------------------------------------------!
     
    29313028       ENDIF
    29323029!
    2933 !--    Tendency of water vapour mixing ratio is obtained from the
    2934 !--    change in RH during SALSA run. This releases heat and changes pt.
    2935 !--    Assumes no temperature change during SALSA run.
     3030!--    Tendency of water vapour mixing ratio is obtained from the change in RH during SALSA run.
     3031!--    This releases heat and changes pt. Assumes no temperature change during SALSA run.
    29363032!--    q = r / (1+r), Euler method for integration
    29373033!
     
    29513047          CALL depo_surf( i, j, surf_def_h(0), vd, schmidt_num, kvis, in_u, .TRUE. )
    29523048          DO  l = 0, 3
    2953              CALL depo_surf( i, j, surf_def_v(l), vd, schmidt_num, kvis, in_u, .FALSE., l )
     3049             CALL depo_surf( i, j, surf_def_v(l), vd, schmidt_num, kvis, in_u, .FALSE. )
    29543050          ENDDO
    29553051       ELSE
    2956           CALL depo_surf( i, j, surf_usm_h, vd, schmidt_num, kvis, in_u, .TRUE. )
     3052          CALL depo_surf( i, j, surf_usm_h, vd, schmidt_num, kvis, in_u, .TRUE., usm_to_depo_h )
    29573053          DO  l = 0, 3
    2958              CALL depo_surf( i, j, surf_usm_v(l), vd, schmidt_num, kvis, in_u, .FALSE., l )
     3054             CALL depo_surf( i, j, surf_usm_v(l), vd, schmidt_num, kvis, in_u, .FALSE.,            &
     3055                             usm_to_depo_v(l) )
    29593056          ENDDO
    2960           CALL depo_surf( i, j, surf_lsm_h, vd, schmidt_num, kvis, in_u, .TRUE. )
     3057          CALL depo_surf( i, j, surf_lsm_h, vd, schmidt_num, kvis, in_u, .TRUE., lsm_to_depo_h )
    29613058          DO  l = 0, 3
    2962              CALL depo_surf( i, j, surf_lsm_v(l), vd, schmidt_num, kvis, in_u, .FALSE., l )
     3059             CALL depo_surf( i, j, surf_lsm_v(l), vd, schmidt_num, kvis, in_u, .FALSE.,            &
     3060                             lsm_to_depo_v(l) )
    29633061          ENDDO
    29643062       ENDIF
     
    33423440    IMPLICIT NONE
    33433441
    3344     INTEGER(iwp) ::  ib     !< loop index
    3345 
    3346     REAL(wp) ::  avis       !< molecular viscocity of air (kg/(m*s))
    3347     REAL(wp) ::  Cc         !< Cunningham slip-flow correction factor
    3348     REAL(wp) ::  Kn         !< Knudsen number
    3349     REAL(wp) ::  lambda     !< molecular mean free path (m)
    3350     REAL(wp) ::  mdiff      !< particle diffusivity coefficient
    3351     REAL(wp) ::  pdn        !< particle density (kg/m3)
    3352     REAL(wp) ::  ustar      !< friction velocity (m/s)
    3353     REAL(wp) ::  va         !< thermal speed of an air molecule (m/s)
    3354     REAL(wp) ::  zdwet      !< wet diameter (m)
     3442    INTEGER(iwp) ::  ib   !< loop index
     3443    INTEGER(iwp) ::  ic   !< loop index
     3444
     3445    REAL(wp) ::  alpha             !< parameter, Table 3 in Z01
     3446    REAL(wp) ::  avis              !< molecular viscocity of air (kg/(m*s))
     3447    REAL(wp) ::  beta_im           !< parameter for turbulent impaction
     3448    REAL(wp) ::  beta              !< Cunningham slip-flow correction factor
     3449    REAL(wp) ::  c_brownian_diff   !< coefficient for Brownian diffusion
     3450    REAL(wp) ::  c_impaction       !< coefficient for inertial impaction
     3451    REAL(wp) ::  c_interception    !< coefficient for interception
     3452    REAL(wp) ::  c_turb_impaction  !< coefficient for turbulent impaction
     3453    REAL(wp) ::  depo              !< deposition velocity (m/s)
     3454    REAL(wp) ::  gamma             !< parameter, Table 3 in Z01
     3455    REAL(wp) ::  Kn                !< Knudsen number
     3456    REAL(wp) ::  lambda            !< molecular mean free path (m)
     3457    REAL(wp) ::  mdiff             !< particle diffusivity coefficient
     3458    REAL(wp) ::  par_a             !< parameter A for the characteristic radius of collectors,
     3459                                   !< Table 3 in Z01
     3460    REAL(wp) ::  par_l             !< obstacle characteristic dimension in P10
     3461    REAL(wp) ::  pdn               !< particle density (kg/m3)
     3462    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)
    33553465
    33563466    REAL(wp), INTENT(in) ::  adn    !< air density (kg/m3)
     
    33683478!
    33693479!-- Initialise
    3370     pdn           = 1500.0_wp    ! default value
    3371     ustar         = 0.0_wp
     3480    depo  = 0.0_wp
     3481    pdn   = 1500.0_wp    ! default value
     3482    ustar = 0.0_wp
    33723483!
    33733484!-- Molecular viscosity of air (Eq. 4.54)
     
    33823493!-- Mean free path (m) (Eq. 15.24)
    33833494    lambda = 2.0_wp * avis / ( adn * va )
     3495!
     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)
    33843508
    33853509    DO  ib = 1, nbins_aerosol
     
    33923516!
    33933517!--    Cunningham slip-flow correction (Eq. 15.30)
    3394        Cc = 1.0_wp + Kn * ( 1.249_wp + 0.42_wp * EXP( -0.87_wp / Kn ) )
     3518       beta = 1.0_wp + Kn * ( 1.249_wp + 0.42_wp * EXP( -0.87_wp / Kn ) )
    33953519
    33963520!--    Particle diffusivity coefficient (Eq. 15.29)
    3397        mdiff = ( abo * tk * Cc ) / ( 3.0_wp * pi * avis * zdwet )
     3521       mdiff = ( abo * tk * beta ) / ( 3.0_wp * pi * avis * zdwet )
    33983522!
    33993523!--    Particle Schmidt number (Eq. 15.36)
     
    34013525!
    34023526!--    Critical fall speed i.e. settling velocity  (Eq. 20.4)
    3403        vc(ib) = MIN( 1.0_wp, terminal_vel( 0.5_wp * zdwet, pdn, adn, avis, Cc) )
     3527       vc(ib) = MIN( 1.0_wp, terminal_vel( 0.5_wp * zdwet, pdn, adn, avis, beta) )
    34043528!
    34053529!--    Friction velocity for deposition on vegetation. Calculated following Prandtl (1925):
    34063530       IF ( lsdepo_pcm  .AND.  plant_canopy  .AND.  lad > 0.0_wp )  THEN
    34073531          ustar = SQRT( cdc ) * mag_u
    3408           CALL depo_pcm( paero, ib, vc(ib), mag_u, ustar, kvis, schmidt_num(ib), lad )
     3532          SELECT CASE ( depo_pcm_par_num )
     3533
     3534             CASE ( 1 )   ! Zhang et al. (2001)
     3535                CALL depo_vel_Z01( vc(ib), ustar, schmidt_num(ib), paero(ib)%dwet, alpha,  gamma,  &
     3536                                   par_a, depo )
     3537             CASE ( 2 )   ! Petroff & Zhang (2010)
     3538                CALL depo_vel_P10( vc(ib), mag_u, ustar, kvis, schmidt_num(ib), paero(ib)%dwet,    &
     3539                                   par_l, c_brownian_diff, c_interception, c_impaction, beta_im,   &
     3540                                   c_turb_impaction, depo )
     3541          END SELECT
     3542!
     3543!--       Calculate the change in concentrations
     3544          paero(ib)%numc = paero(ib)%numc - depo * lad * paero(ib)%numc * dt_salsa
     3545          DO  ic = 1, maxspec+1
     3546             paero(ib)%volc(ic) = paero(ib)%volc(ic) - depo * lad * paero(ib)%volc(ic) * dt_salsa
     3547          ENDDO
    34093548       ENDIF
    34103549    ENDDO
     
    34153554! Description:
    34163555! ------------
    3417 !> Calculate change in number and volume concentrations due to deposition on
    3418 !> plant canopy.
    3419 !------------------------------------------------------------------------------!
    3420  SUBROUTINE depo_pcm( paero, ib, vc, mag_u, ustar, kvis_a, schmidt_num, lad )
     3556!> Calculate deposition velocity (m/s) based on Zhan et al. (2001, case 1).
     3557!------------------------------------------------------------------------------!
     3558
     3559 SUBROUTINE depo_vel_Z01( vc, ustar, schmidt_num, diameter, alpha, gamma, par_a, depo )
    34213560
    34223561    IMPLICIT NONE
    34233562
    3424     INTEGER(iwp) ::  ic      !< loop index
    3425 
    3426     INTEGER(iwp), INTENT(in) ::  ib  !< loop index
    3427 
    3428     REAL(wp) ::  alpha             !< parameter, Table 3 in Z01
    3429     REAL(wp) ::  beta_im           !< parameter for turbulent impaction
    3430     REAL(wp) ::  depo              !< deposition efficiency
    3431     REAL(wp) ::  c_brownian_diff   !< coefficient for Brownian diffusion
    3432     REAL(wp) ::  c_impaction       !< coefficient for inertial impaction
    3433     REAL(wp) ::  c_interception    !< coefficient for interception
    3434     REAL(wp) ::  c_turb_impaction  !< coefficient for turbulent impaction
    3435     REAL(wp) ::  gamma             !< parameter, Table 3 in Z01
    3436     REAL(wp) ::  par_a             !< parameter A for the characteristic radius of collectors,
    3437                                    !< Table 3 in Z01
    3438     REAL(wp) ::  par_l             !< obstacle characteristic dimension in P10
    34393563    REAL(wp) ::  rs                !< overall quasi-laminar resistance for particles
     3564    REAL(wp) ::  stokes_num        !< Stokes number for smooth or bluff surfaces
     3565
     3566    REAL(wp), INTENT(in) ::  alpha        !< parameter, Table 3 in Z01
     3567    REAL(wp), INTENT(in) ::  gamma        !< parameter, Table 3 in Z01
     3568    REAL(wp), INTENT(in) ::  par_a        !< parameter A for the characteristic diameter of
     3569                                          !< collectors, Table 3 in Z01
     3570    REAL(wp), INTENT(in) ::  diameter     !< particle diameter
     3571    REAL(wp), INTENT(in) ::  schmidt_num  !< particle Schmidt number
     3572    REAL(wp), INTENT(in) ::  ustar        !< friction velocity (m/s)
     3573    REAL(wp), INTENT(in) ::  vc           !< terminal velocity (m/s)
     3574
     3575    REAL(wp), INTENT(inout)  ::  depo     !< deposition efficiency (m/s)
     3576
     3577    IF ( par_a > 0.0_wp )  THEN
     3578!
     3579!--    Initialise
     3580       rs = 0.0_wp
     3581!
     3582!--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24)
     3583       stokes_num = vc * ustar / ( g * par_a )
     3584!
     3585!--    The overall quasi-laminar resistance for particles (Zhang et al., Eq. 5)
     3586       rs = MAX( EPSILON( 1.0_wp ), ( 3.0_wp * ustar * EXP( -stokes_num**0.5_wp ) *                &
     3587                 ( schmidt_num**( -gamma ) + ( stokes_num / ( alpha + stokes_num ) )**2 +          &
     3588                 0.5_wp * ( diameter / par_a )**2 ) ) )
     3589
     3590       depo = rs + vc
     3591
     3592    ELSE
     3593       depo = 0.0_wp
     3594    ENDIF
     3595
     3596 END SUBROUTINE depo_vel_Z01
     3597
     3598!------------------------------------------------------------------------------!
     3599! Description:
     3600! ------------
     3601!> Calculate deposition velocity (m/s) based on Petroff & Zhang (2010, case 2).
     3602!------------------------------------------------------------------------------!
     3603
     3604 SUBROUTINE depo_vel_P10( vc, mag_u, ustar, kvis_a, schmidt_num, diameter, par_l, c_brownian_diff, &
     3605                          c_interception, c_impaction, beta_im, c_turb_impaction, depo )
     3606
     3607    IMPLICIT NONE
     3608
    34403609    REAL(wp) ::  stokes_num        !< Stokes number for smooth or bluff surfaces
    34413610    REAL(wp) ::  tau_plus          !< dimensionless particle relaxation time
     
    34453614    REAL(wp) ::  v_it              !< deposition velocity due to turbulent impaction
    34463615
     3616    REAL(wp), INTENT(in) ::  beta_im           !< parameter for turbulent impaction
     3617    REAL(wp), INTENT(in) ::  c_brownian_diff   !< coefficient for Brownian diffusion
     3618    REAL(wp), INTENT(in) ::  c_impaction       !< coefficient for inertial impaction
     3619    REAL(wp), INTENT(in) ::  c_interception    !< coefficient for interception
     3620    REAL(wp), INTENT(in) ::  c_turb_impaction  !< coefficient for turbulent impaction
    34473621    REAL(wp), INTENT(in) ::  kvis_a       !< kinematic viscosity of air (m2/s)
    3448     REAL(wp), INTENT(in) ::  lad          !< leaf area density (m2/m3)
    34493622    REAL(wp), INTENT(in) ::  mag_u        !< wind velocity (m/s)
    3450     REAL(wp), INTENT(in) ::  schmidt_num  !< particle Schmidt number
     3623    REAL(wp), INTENT(in) ::  par_l        !< obstacle characteristic dimension in P10
     3624    REAL(wp), INTENT(in) ::  diameter       !< particle diameter
     3625    REAL(wp), INTENT(in) ::  schmidt_num  !< particle Schmidt number
    34513626    REAL(wp), INTENT(in) ::  ustar        !< friction velocity (m/s)
    34523627    REAL(wp), INTENT(in) ::  vc           !< terminal velocity (m/s)
    34533628
    3454     TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero  !< aerosol properties
    3455 !
    3456 !-- Initialise
    3457     rs       = 0.0_wp
    3458     tau_plus = 0.0_wp
    3459     v_bd     = 0.0_wp
    3460     v_im     = 0.0_wp
    3461     v_in     = 0.0_wp
    3462     v_it     = 0.0_wp
    3463 
    3464     IF ( depo_pcm_par == 'zhang2001' )  THEN
    3465 !
    3466 !--    Parameters for the land use category 'deciduous broadleaf trees'(Table 3)
    3467        alpha   = alpha_z01(depo_pcm_type_num)
    3468        gamma   = gamma_z01(depo_pcm_type_num)
    3469        par_a   = A_z01(depo_pcm_type_num, season) * 1.0E-3_wp
    3470 !
    3471 !--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24)
    3472        stokes_num = vc * ustar / ( g * par_a )
    3473 !
    3474 !--    The overall quasi-laminar resistance for particles (Zhang et al., Eq. 5)
    3475        rs = MAX( EPSILON( 1.0_wp ), ( 3.0_wp * ustar * EXP( -stokes_num**0.5_wp ) *                &
    3476                  ( schmidt_num**( -gamma ) + ( stokes_num / ( alpha + stokes_num ) )**2 +          &
    3477                  0.5_wp * ( paero(ib)%dwet / par_a )**2 ) ) )
    3478 
    3479        depo = ( rs + vc ) * lad
    3480 
    3481     ELSEIF ( depo_pcm_par == 'petroff2010' )  THEN
    3482 !
    3483 !--    vd = v_BD + v_IN + v_IM + v_IT + vc
    3484 !--    Deposition efficiencies from Table 1. Constants from Table 2.
    3485        par_l   = l_p10(depo_pcm_type_num) * 0.01_wp
    3486        c_brownian_diff     = c_b_p10(depo_pcm_type_num)
    3487        c_interception    = c_in_p10(depo_pcm_type_num)
    3488        c_impaction    = c_im_p10(depo_pcm_type_num)
    3489        beta_im = beta_im_p10(depo_pcm_type_num)
    3490        c_turb_impaction    = c_it_p10(depo_pcm_type_num)
     3629    REAL(wp), INTENT(inout)  ::  depo     !< deposition efficiency (m/s)
     3630
     3631    IF ( par_l > 0.0_wp )  THEN
     3632!
     3633!--    Initialise
     3634       tau_plus = 0.0_wp
     3635       v_bd     = 0.0_wp
     3636       v_im     = 0.0_wp
     3637       v_in     = 0.0_wp
     3638       v_it     = 0.0_wp
    34913639!
    34923640!--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24)
     
    34973645!
    34983646!--    Brownian diffusion
    3499        v_bd = mag_u * c_brownian_diff * schmidt_num**( -0.66666666_wp ) *                          &
     3647       v_bd = mag_u * c_brownian_diff * schmidt_num**( -0.66666666_wp ) *                             &
    35003648              ( mag_u * par_l / kvis_a )**( -0.5_wp )
    35013649!
    35023650!--    Interception
    3503        v_in = mag_u * c_interception * paero(ib)%dwet / par_l * ( 2.0_wp + LOG( 2.0_wp * par_l /    &
    3504               paero(ib)%dwet ) )
     3651       v_in = mag_u * c_interception * diameter / par_l * ( 2.0_wp + LOG( 2.0_wp * par_l / diameter ) )
    35053652!
    35063653!--    Impaction: Petroff (2009) Eq. 18
     
    35143661       ENDIF
    35153662
    3516        depo = ( v_bd + v_in + v_im + v_it + vc ) * lad
    3517 
     3663       depo = ( v_bd + v_in + v_im + v_it + vc )
     3664
     3665    ELSE
     3666       depo = 0.0_wp
    35183667    ENDIF
    3519 !
    3520 !-- Calculate the change in concentrations
    3521     paero(ib)%numc = paero(ib)%numc - depo * paero(ib)%numc * dt_salsa
    3522     DO  ic = 1, maxspec+1
    3523        paero(ib)%volc(ic) = paero(ib)%volc(ic) - depo * paero(ib)%volc(ic) * dt_salsa
    3524     ENDDO
    3525 
    3526  END SUBROUTINE depo_pcm
     3668
     3669 END SUBROUTINE depo_vel_P10
    35273670
    35283671!------------------------------------------------------------------------------!
     
    35343677!        high-resolution simulations)
    35353678!------------------------------------------------------------------------------!
    3536  SUBROUTINE depo_surf( i, j, surf, vc, schmidt_num, kvis, mag_u, norm, l )
    3537 
    3538     USE arrays_3d,                                                             &
     3679 SUBROUTINE depo_surf( i, j, surf, vc, schmidt_num, kvis, mag_u, norm, match_array )
     3680
     3681    USE arrays_3d,                                                                                 &
    35393682        ONLY: rho_air_zw
    35403683
    3541     USE surface_mod,                                                           &
    3542         ONLY:  surf_type
     3684    USE surface_mod,                                                                               &
     3685        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_type
    35433686
    35443687    IMPLICIT NONE
     
    35523695    INTEGER(iwp) ::  surf_s  !< Start index of surface elements at (j,i)-gridpoint
    35533696
    3554     INTEGER(iwp), INTENT(in) ::  i     !< loop index
    3555     INTEGER(iwp), INTENT(in) ::  j     !< loop index
    3556 
    3557     INTEGER(iwp), INTENT(in), OPTIONAL ::  l  !< index variable for surface facing
    3558 
    3559     LOGICAL, INTENT(in) ::  norm      !< to normalise or not
     3697    INTEGER(iwp), INTENT(in) ::  i  !< loop index
     3698    INTEGER(iwp), INTENT(in) ::  j  !< loop index
     3699
     3700    LOGICAL, INTENT(in) ::  norm   !< to normalise or not
    35603701
    35613702    REAL(wp) ::  alpha             !< parameter, Table 3 in Z01
     
    35653706    REAL(wp) ::  c_interception    !< coefficient for interception
    35663707    REAL(wp) ::  c_turb_impaction  !< coefficient for turbulent impaction
    3567     REAL(wp) ::  depo              !< deposition efficiency
    35683708    REAL(wp) ::  gamma             !< parameter, Table 3 in Z01
    35693709    REAL(wp) ::  norm_fac          !< normalisation factor (usually air density)
     
    35723712    REAL(wp) ::  par_l             !< obstacle characteristic dimension in P10
    35733713    REAL(wp) ::  rs                !< the overall quasi-laminar resistance for particles
    3574     REAL(wp) ::  stokes_num        !< Stokes number for bluff surface elements
    35753714    REAL(wp) ::  tau_plus          !< dimensionless particle relaxation time
    35763715    REAL(wp) ::  v_bd              !< deposition velocity due to Brownian diffusion
     
    35793718    REAL(wp) ::  v_it              !< deposition velocity due to turbulent impaction
    35803719
     3720    REAL(wp), DIMENSION(nbins_aerosol) ::  depo      !< deposition efficiency
     3721    REAL(wp), DIMENSION(nbins_aerosol) ::  depo_sum  !< sum of deposition efficiencies
     3722
    35813723    REAL(wp), DIMENSION(:), INTENT(in) ::  kvis   !< kinematic viscosity of air (m2/s)
    35823724    REAL(wp), DIMENSION(:), INTENT(in) ::  mag_u  !< wind velocity (m/s)
     
    35853727    REAL(wp), DIMENSION(:,:), INTENT(in) ::  vc            !< terminal velocity (m/s)
    35863728
    3587     TYPE(surf_type), INTENT(inout) :: surf  !< respective surface type
     3729    TYPE(match_surface), INTENT(in), OPTIONAL ::  match_array  !< match the deposition module and
     3730                                                               !< LSM/USM surfaces
     3731    TYPE(surf_type), INTENT(inout) :: surf                     !< respective surface type
    35883732!
    35893733!-- Initialise
     3734    depo     = 0.0_wp
     3735    depo_sum = 0.0_wp
    35903736    rs       = 0.0_wp
    35913737    surf_s   = surf%start_index(j,i)
     
    35973743    v_it     = 0.0_wp
    35983744!
    3599 !-- Model parameters for the land use category. If LSM is applied, import
     3745!-- Model parameters for the land use category. If LSM or USM is applied, import
    36003746!-- characteristics. Otherwise, apply surface type "urban".
    36013747    alpha   = alpha_z01(luc_urban)
     
    36103756    c_turb_impaction = c_it_p10(luc_urban)
    36113757
    3612     DO  m = surf_s, surf_e
    3613        k = surf%k(m)
    3614 
    3615        IF ( norm )  THEN
    3616           norm_fac = rho_air_zw(k)
    3617           IF ( land_surface )  THEN
    3618              alpha            = alpha_z01( lsm_to_depo_h%match(m) )
    3619              beta_im          = beta_im_p10( lsm_to_depo_h%match(m) )
    3620              c_brownian_diff  = c_b_p10( lsm_to_depo_h%match(m) )
    3621              c_impaction      = c_im_p10( lsm_to_depo_h%match(m) )
    3622              c_interception   = c_in_p10( lsm_to_depo_h%match(m) )
    3623              c_turb_impaction = c_it_p10( lsm_to_depo_h%match(m) )
    3624              gamma            = gamma_z01( lsm_to_depo_h%match(m) )
    3625              par_a            = A_z01( lsm_to_depo_h%match(m), season ) * 1.0E-3_wp
    3626              par_l            = l_p10( lsm_to_depo_h%match(m) ) * 0.01_wp
     3758
     3759    IF ( PRESENT( match_array ) )  THEN  ! land or urban surface model
     3760
     3761       DO  m = surf_s, surf_e
     3762
     3763          k = surf%k(m)
     3764          norm_fac = 1.0_wp
     3765
     3766          IF ( norm )  norm_fac = rho_air_zw(k)  ! normalise vertical fluxes by air density
     3767
     3768          IF ( match_array%match_lupg(m) > 0 )  THEN
     3769             alpha = alpha_z01( match_array%match_lupg(m) )
     3770             gamma = gamma_z01( match_array%match_lupg(m) )
     3771             par_a = A_z01( match_array%match_lupg(m), season ) * 1.0E-3_wp
     3772
     3773             beta_im          = beta_im_p10( match_array%match_lupg(m) )
     3774             c_brownian_diff  = c_b_p10( match_array%match_lupg(m) )
     3775             c_impaction      = c_im_p10( match_array%match_lupg(m) )
     3776             c_interception   = c_in_p10( match_array%match_lupg(m) )
     3777             c_turb_impaction = c_it_p10( match_array%match_lupg(m) )
     3778             par_l            = l_p10( match_array%match_lupg(m) ) * 0.01_wp
     3779
     3780             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
     3783
     3784                SELECT CASE ( depo_surf_par_num )
     3785
     3786                   CASE ( 1 )
     3787                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
     3788                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
     3789                   CASE ( 2 )
     3790                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
     3791                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
     3792                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
     3793                                         c_turb_impaction, depo(ib) )
     3794                END SELECT
     3795             ENDDO
     3796             depo_sum = depo_sum + surf%frac(ind_pav_green,m) * depo
    36273797          ENDIF
    3628        ELSE
    3629           norm_fac = 0.0_wp
    3630           IF ( land_surface )  THEN
    3631              alpha            = alpha_z01( lsm_to_depo_v(l)%match(m) )
    3632              beta_im          = beta_im_p10( lsm_to_depo_v(l)%match(m) )
    3633              c_brownian_diff  = c_b_p10( lsm_to_depo_v(l)%match(m) )
    3634              c_impaction      = c_im_p10( lsm_to_depo_v(l)%match(m) )
    3635              c_interception   = c_in_p10( lsm_to_depo_v(l)%match(m) )
    3636              c_turb_impaction = c_it_p10( lsm_to_depo_v(l)%match(m) )
    3637              gamma            = gamma_z01( lsm_to_depo_v(l)%match(m) )
    3638              par_a            = A_z01( lsm_to_depo_v(l)%match(m), season ) * 1.0E-3_wp
    3639              par_l            = l_p10( lsm_to_depo_v(l)%match(m) ) * 0.01_wp
     3798
     3799          IF ( match_array%match_luvw(m) > 0 )  THEN
     3800             alpha = alpha_z01( match_array%match_luvw(m) )
     3801             gamma = gamma_z01( match_array%match_luvw(m) )
     3802             par_a = A_z01( match_array%match_luvw(m), season ) * 1.0E-3_wp
     3803
     3804             beta_im          = beta_im_p10( match_array%match_luvw(m) )
     3805             c_brownian_diff  = c_b_p10( match_array%match_luvw(m) )
     3806             c_impaction      = c_im_p10( match_array%match_luvw(m) )
     3807             c_interception   = c_in_p10( match_array%match_luvw(m) )
     3808             c_turb_impaction = c_it_p10( match_array%match_luvw(m) )
     3809             par_l            = l_p10( match_array%match_luvw(m) ) * 0.01_wp
     3810
     3811             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
     3814
     3815                SELECT CASE ( depo_surf_par_num )
     3816
     3817                   CASE ( 1 )
     3818                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
     3819                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
     3820                   CASE ( 2 )
     3821                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
     3822                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
     3823                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
     3824                                         c_turb_impaction, depo(ib) )
     3825                END SELECT
     3826             ENDDO
     3827             depo_sum = depo_sum + surf%frac(ind_veg_wall,m) * depo
    36403828          ENDIF
    3641        ENDIF
    3642 
    3643        DO  ib = 1, nbins_aerosol
    3644           IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
    3645 
    3646           IF ( depo_surf_par == 'zhang2001' )  THEN
    3647 !
    3648 !--          Stokes number for smooth surfaces or surfaces with bluff roughness
    3649 !--          elements (Seinfeld and Pandis, 2nd edition (2006): Eq. 19.23)
    3650              stokes_num = MAX( 0.01_wp, vc(k+1,ib) * surf%us(m)**2 / ( g * kvis(k+1)  ) )
    3651 !
    3652 !--          The overall quasi-laminar resistance for particles (Eq. 5)
    3653              rs = MAX( EPSILON( 1.0_wp ), ( 3.0_wp * surf%us(m) * ( schmidt_num(k+1,ib)**( -gamma )&
    3654                        + ( stokes_num / ( alpha + stokes_num ) )**2 + 0.5_wp * ( ra_dry(k,j,i,ib) /&
    3655                        par_a )**2 ) * EXP( -stokes_num**0.5_wp ) ) )
    3656              depo = vc(k+1,ib) + rs
    3657 
    3658           ELSEIF ( depo_surf_par == 'petroff2010' )  THEN
    3659 !
    3660 !--          vd = v_BD + v_IN + v_IM + v_IT + vc
    3661 !
    3662 !--          Stokes number for smooth surfaces or surfaces with bluff roughness
    3663 !--          elements (Seinfeld and Pandis, 2nd edition (2006): Eq. 19.23)
    3664              stokes_num = MAX( 0.01_wp, vc(k+1,ib) * surf%us(m)**2 / ( g *  kvis(k+1) ) )
    3665 !
    3666 !--          Non-dimensional relexation time of the particle on top of canopy
    3667              tau_plus = vc(k+1,ib) * surf%us(m)**2 / ( kvis(k+1) * g )
    3668 !
    3669 !--          Brownian diffusion
    3670              v_bd = mag_u(k+1) * c_brownian_diff * schmidt_num(k+1,ib)**( -0.666666_wp ) *         &
    3671                     ( mag_u(k+1) * par_l / kvis(k+1) )**( -0.5_wp )
    3672 !
    3673 !--          Interception
    3674              v_in = mag_u(k+1) * c_interception * ra_dry(k,j,i,ib)/ par_l *                        &
    3675                     ( 2.0_wp + LOG( 2.0_wp * par_l / ra_dry(k,j,i,ib) ) )
    3676 !
    3677 !--          Impaction: Petroff (2009) Eq. 18
    3678              v_im = mag_u(k+1) * c_impaction * ( stokes_num / ( stokes_num + beta_im ) )**2
    3679 
    3680              IF ( tau_plus < 20.0_wp )  THEN
    3681                 v_it = 2.5E-3_wp * c_turb_impaction * tau_plus**2
     3829
     3830          IF ( match_array%match_luww(m) > 0 )  THEN
     3831             alpha = alpha_z01( match_array%match_luww(m) )
     3832             gamma = gamma_z01( match_array%match_luww(m) )
     3833             par_a = A_z01( match_array%match_luww(m), season ) * 1.0E-3_wp
     3834
     3835             beta_im          = beta_im_p10( match_array%match_luww(m) )
     3836             c_brownian_diff  = c_b_p10( match_array%match_luww(m) )
     3837             c_impaction      = c_im_p10( match_array%match_luww(m) )
     3838             c_interception   = c_in_p10( match_array%match_luww(m) )
     3839             c_turb_impaction = c_it_p10( match_array%match_luww(m) )
     3840             par_l            = l_p10( match_array%match_luww(m) ) * 0.01_wp
     3841
     3842             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
     3845
     3846                SELECT CASE ( depo_surf_par_num )
     3847
     3848                   CASE ( 1 )
     3849                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
     3850                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
     3851                   CASE ( 2 )
     3852                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
     3853                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
     3854                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
     3855                                         c_turb_impaction, depo(ib) )
     3856                END SELECT
     3857             ENDDO
     3858             depo_sum = depo_sum + surf%frac(ind_wat_win,m) * depo
     3859          ENDIF
     3860
     3861          DO  ib = 1, nbins_aerosol
     3862!
     3863!--          Calculate changes in surface fluxes due to dry deposition
     3864             IF ( include_emission )  THEN
     3865                surf%answs(m,ib) = aerosol_number(ib)%source(j,i) - MAX( 0.0_wp,                   &
     3866                                   depo_sum(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i) )
     3867                DO  ic = 1, ncomponents_mass
     3868                   icc = ( ic - 1 ) * nbins_aerosol + ib
     3869                   surf%amsws(m,icc) = aerosol_mass(icc)%source(j,i) - MAX( 0.0_wp,                &
     3870                                       depo_sum(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i) )
     3871                ENDDO  ! ic
    36823872             ELSE
    3683                 v_it = c_turb_impaction
     3873                surf%answs(m,ib) = -depo_sum(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i)
     3874                DO  ic = 1, ncomponents_mass
     3875                   icc = ( ic - 1 ) * nbins_aerosol + ib
     3876                   surf%amsws(m,icc) = -depo_sum(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i)
     3877                ENDDO  ! ic
    36843878             ENDIF
    3685              depo =  v_bd + v_in + v_im + v_it + vc(k+1,ib)
    3686 
    3687           ENDIF
    3688 !
    3689 !--       Calculate changes in surface fluxes due to dry deposition
    3690           IF ( include_emission )  THEN
    3691              surf%answs(m,ib) = aerosol_number(ib)%source(j,i) -                                   &
    3692                                 MAX( 0.0_wp, depo * norm_fac * aerosol_number(ib)%conc(k,j,i) )
    3693              DO  ic = 1, ncomponents_mass
    3694                 icc = ( ic - 1 ) * nbins_aerosol + ib
    3695                 surf%amsws(m,icc) = aerosol_mass(icc)%source(j,i) -                                &
    3696                                     MAX( 0.0_wp, depo *  norm_fac * aerosol_mass(icc)%conc(k,j,i) )
    3697              ENDDO  ! ic
    3698           ELSE
    3699              surf%answs(m,ib) = -depo * norm_fac * aerosol_number(ib)%conc(k,j,i)
    3700              DO  ic = 1, ncomponents_mass
    3701                 icc = ( ic - 1 ) * nbins_aerosol + ib
    3702                 surf%amsws(m,icc) = -depo *  norm_fac * aerosol_mass(icc)%conc(k,j,i)
    3703              ENDDO    ! ic
    3704           ENDIF
    3705        ENDDO    ! ib
    3706     ENDDO    ! m
     3879          ENDDO  ! ib
     3880
     3881       ENDDO
     3882
     3883    ELSE  ! default surfaces
     3884
     3885       DO  m = surf_s, surf_e
     3886
     3887          k = surf%k(m)
     3888          norm_fac = 1.0_wp
     3889
     3890          IF ( norm )  norm_fac = rho_air_zw(k)  ! normalise vertical fluxes by air density
     3891
     3892          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
     3895
     3896             SELECT CASE ( depo_surf_par_num )
     3897
     3898                CASE ( 1 )
     3899                   CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),                 &
     3900                                      ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
     3901                CASE ( 2 )
     3902                   CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),               &
     3903                                      schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,                &
     3904                                      c_brownian_diff, c_interception, c_impaction, beta_im,       &
     3905                                      c_turb_impaction, depo(ib) )
     3906             END SELECT
     3907!
     3908!--          Calculate changes in surface fluxes due to dry deposition
     3909             IF ( include_emission )  THEN
     3910                surf%answs(m,ib) = aerosol_number(ib)%source(j,i) - MAX( 0.0_wp,                   &
     3911                                   depo(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i) )
     3912                DO  ic = 1, ncomponents_mass
     3913                   icc = ( ic - 1 ) * nbins_aerosol + ib
     3914                   surf%amsws(m,icc) = aerosol_mass(icc)%source(j,i) - MAX( 0.0_wp,                &
     3915                                       depo(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i) )
     3916                ENDDO  ! ic
     3917             ELSE
     3918                surf%answs(m,ib) = -depo(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i)
     3919                DO  ic = 1, ncomponents_mass
     3920                   icc = ( ic - 1 ) * nbins_aerosol + ib
     3921                   surf%amsws(m,icc) = -depo(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i)
     3922                ENDDO  ! ic
     3923             ENDIF
     3924          ENDDO  ! ib
     3925       ENDDO
     3926
     3927    ENDIF
    37073928
    37083929 END SUBROUTINE depo_surf
     
    72007421 END SUBROUTINE salsa_actions_ij
    72017422
    7202 !
     7423!------------------------------------------------------------------------------!
     7424! Description:
     7425! ------------
     7426!> Call for all grid points
     7427!------------------------------------------------------------------------------!
     7428 SUBROUTINE salsa_non_advective_processes
     7429
     7430    USE cpulog,                                                                                    &
     7431        ONLY:  cpu_log, log_point_s
     7432
     7433    IMPLICIT NONE
     7434
     7435    INTEGER(iwp) ::  i  !<
     7436    INTEGER(iwp) ::  j  !<
     7437
     7438    IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
     7439       IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  THEN
     7440!
     7441!--       Calculate aerosol dynamic processes. salsa_driver can be run with a longer time step.
     7442          CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
     7443          DO  i = nxl, nxr
     7444             DO  j = nys, nyn
     7445                CALL salsa_diagnostics( i, j )
     7446                CALL salsa_driver( i, j, 3 )
     7447                CALL salsa_diagnostics( i, j )
     7448             ENDDO
     7449          ENDDO
     7450          CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
     7451       ENDIF
     7452    ENDIF
     7453
     7454 END SUBROUTINE salsa_non_advective_processes
     7455
     7456
     7457!------------------------------------------------------------------------------!
     7458! Description:
     7459! ------------
     7460!> Call for grid points i,j
     7461!------------------------------------------------------------------------------!
     7462 SUBROUTINE salsa_non_advective_processes_ij( i, j )
     7463
     7464    USE cpulog,                                                                &
     7465        ONLY:  cpu_log, log_point_s
     7466
     7467    IMPLICIT NONE
     7468
     7469    INTEGER(iwp), INTENT(IN) ::  i  !< grid index in x-direction
     7470    INTEGER(iwp), INTENT(IN) ::  j  !< grid index in y-direction
     7471
     7472    IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
     7473       IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  THEN
     7474!
     7475!--       Calculate aerosol dynamic processes. salsa_driver can be run with a longer time step.
     7476          CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
     7477          CALL salsa_diagnostics( i, j )
     7478          CALL salsa_driver( i, j, 3 )
     7479          CALL salsa_diagnostics( i, j )
     7480          CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
     7481       ENDIF
     7482    ENDIF
     7483
     7484 END SUBROUTINE salsa_non_advective_processes_ij
     7485
     7486!------------------------------------------------------------------------------!
     7487! Description:
     7488! ------------
     7489!> Routine for exchange horiz of salsa variables.
     7490!------------------------------------------------------------------------------!
     7491 SUBROUTINE salsa_exchange_horiz_bounds
     7492
     7493    USE cpulog,                                                                &
     7494        ONLY:  cpu_log, log_point_s
     7495
     7496    IMPLICIT NONE
     7497
     7498    INTEGER(iwp) ::  ib   !<
     7499    INTEGER(iwp) ::  ic   !<
     7500    INTEGER(iwp) ::  icc  !<
     7501    INTEGER(iwp) ::  ig   !<
     7502
     7503    IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
     7504       IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  THEN
     7505
     7506          CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
     7507!
     7508!--       Exchange ghost points and decycle if needed.
     7509          DO  ib = 1, nbins_aerosol
     7510             CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
     7511             CALL salsa_boundary_conds( aerosol_number(ib)%conc, aerosol_number(ib)%init )
     7512             DO  ic = 1, ncomponents_mass
     7513                icc = ( ic - 1 ) * nbins_aerosol + ib
     7514                CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
     7515                CALL salsa_boundary_conds( aerosol_mass(icc)%conc, aerosol_mass(icc)%init )
     7516             ENDDO
     7517          ENDDO
     7518          IF ( .NOT. salsa_gases_from_chem )  THEN
     7519             DO  ig = 1, ngases_salsa
     7520                CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
     7521                CALL salsa_boundary_conds( salsa_gas(ig)%conc, salsa_gas(ig)%init )
     7522             ENDDO
     7523          ENDIF
     7524          CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
     7525!
     7526!--       Update last_salsa_time
     7527          last_salsa_time = time_since_reference_point
     7528       ENDIF
     7529    ENDIF
     7530
     7531 END SUBROUTINE salsa_exchange_horiz_bounds
     7532
    72037533!------------------------------------------------------------------------------!
    72047534! Description:
     
    72727602! ------------
    72737603!> Calculate the prognostic equation for aerosol number and mass, and gas
    7274 !> concentrations. Cache-optimized.
     7604!> concentrations. For vector machines.
    72757605!------------------------------------------------------------------------------!
    72767606 SUBROUTINE salsa_prognostic_equations()
     
    73797709    REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  flux_l  !<
    73807710
    7381     REAL(wp), DIMENSION(:,:,:), POINTER ::  rs_p    !<
    7382     REAL(wp), DIMENSION(:,:,:), POINTER ::  rs      !<
    7383     REAL(wp), DIMENSION(:,:,:), POINTER ::  trs_m   !<
     7711    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  rs_p    !<
     7712    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  rs      !<
     7713    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  trs_m   !<
    73847714
    73857715    icc = ( ic - 1 ) * nbins_aerosol + ib
     
    74337763    END SELECT
    74347764!
    7435 !-- Sedimentation for aerosol number and mass
     7765!-- Sedimentation and prognostic equation for aerosol number and mass
    74367766    IF ( lsdepo  .AND.  do_sedimentation )  THEN
    7437        tend(nzb+1:nzt,j,i) = tend(nzb+1:nzt,j,i) - MAX( 0.0_wp,                                    &
    7438                              ( rs(nzb+2:nzt+1,j,i) * sedim_vd(nzb+2:nzt+1,j,i,ib) -                &
    7439                                rs(nzb+1:nzt,j,i) * sedim_vd(nzb+1:nzt,j,i,ib) ) * ddzu(nzb+1:nzt) )&
    7440                              * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(nzb:nzt-1,j,i), 0 ) )
     7767!DIR$ IVDEP
     7768       DO  k = nzb+1, nzt
     7769          tend(k,j,i) = tend(k,j,i) - MAX( 0.0_wp, ( rs(k+1,j,i) * sedim_vd(k+1,j,i,ib) -          &
     7770                                                     rs(k,j,i) * sedim_vd(k,j,i,ib) ) * ddzu(k) )  &
     7771                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k-1,j,i), 0 ) )
     7772          rs_p(k,j,i) = rs(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) * trs_m(k,j,i) )     &
     7773                                      - tsc(5) * rdf_sc(k) * ( rs(k,j,i) - rs_init(k) ) )          &
     7774                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     7775          IF ( rs_p(k,j,i) < 0.0_wp )  rs_p(k,j,i) = 0.1_wp * rs(k,j,i)
     7776       ENDDO
     7777    ELSE
     7778!
     7779!--    Prognostic equation
     7780!DIR$ IVDEP
     7781       DO  k = nzb+1, nzt
     7782          rs_p(k,j,i) = rs(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) * trs_m(k,j,i) )     &
     7783                                                - tsc(5) * rdf_sc(k) * ( rs(k,j,i) - rs_init(k) ) )&
     7784                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     7785          IF ( rs_p(k,j,i) < 0.0_wp )  rs_p(k,j,i) = 0.1_wp * rs(k,j,i)
     7786       ENDDO
    74417787    ENDIF
    7442 !
    7443 !-- Prognostic equation for a scalar
    7444     DO  k = nzb+1, nzt
    7445        rs_p(k,j,i) = rs(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) * trs_m(k,j,i) )        &
    7446                                            - tsc(5) * rdf_sc(k) * ( rs(k,j,i) - rs_init(k) ) )     &
    7447                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    7448        IF ( rs_p(k,j,i) < 0.0_wp )  rs_p(k,j,i) = 0.1_wp * rs(k,j,i)
    7449     ENDDO
    74507788!
    74517789!-- Calculate tendencies for the next Runge-Kutta step
     
    74687806! ------------
    74697807!> Calculate the tendencies for aerosol number and mass concentrations.
    7470 !> Vector-optimized.
     7808!> For vector machines.
    74717809!------------------------------------------------------------------------------!
    74727810 SUBROUTINE salsa_tendency( id, rs_p, rs, trs_m, ib, ic, rs_init, do_sedimentation )
     
    90849422    SELECT CASE ( TRIM( var ) )
    90859423
    9086        CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV', 'N_bin1', 'N_bin2', 'N_bin3',     &
    9087               'N_bin4',  'N_bin5',  'N_bin6', 'N_bin7', 'N_bin8', 'N_bin9', 'N_bin10', 'N_bin11',  &
    9088               'N_bin12', 'Ntot' )
     9424       CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV' )
    90899425          IF (  .NOT.  salsa )  THEN
    90909426             message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.'
    9091              CALL message( 'check_parameters', 'PA0645', 1, 2, 0, 6, 0 )
     9427             CALL message( 'check_parameters', 'PA0652', 1, 2, 0, 6, 0 )
     9428          ENDIF
     9429          IF (  salsa_gases_from_chem )  THEN
     9430             message_string = 'gases are imported from the chemistry module and thus output of "'  &
     9431                               // TRIM( var ) // '" is not allowed'
     9432             CALL message( 'check_parameters', 'PA0653', 1, 2, 0, 6, 0 )
    90929433          ENDIF
    90939434          unit = '#/m3'
     
    91089449          ENDIF
    91099450          unit = 'kg/m3'
     9451
     9452       CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4',  'N_bin5',  'N_bin6', 'N_bin7', 'N_bin8',    &
     9453              'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12', 'Ntot' )
     9454          IF (  .NOT.  salsa )  THEN
     9455             message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.'
     9456             CALL message( 'check_parameters', 'PA0645', 1, 2, 0, 6, 0 )
     9457          ENDIF
     9458          unit = '#/m3'
    91109459
    91119460       CASE DEFAULT
Note: See TracChangeset for help on using the changeset viewer.