Changeset 4416 for palm/trunk/SOURCE


Ignore:
Timestamp:
Feb 20, 2020 5:53:57 PM (5 years ago)
Author:
monakurppa
Message:

Bug fixes in restart and reformatting av data output

  • add missing arrays (averaged data output) in salsa_wrd_local and salsa_rrd_local
  • set write_binary_salsa and read_restart_data_salsa to .T. by default
  • restructure the average arrays for gases and total mass concentrations of chemical components: set to 4d arrays instead of separate arrays
  • add allocation checks for averaged data output arrays
File:
1 edited

Legend:

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

    r4399 r4416  
    2121! Current revisions:
    2222! -----------------
    23 !
     23! Bug fixes and reformatting for the restart data and averaged data output
     24! - add missing arrays (averaged data output) in salsa_wrd_local and
     25!   salsa_rrd_local
     26! - set write_binary_salsa and read_restart_data_salsa to .T. by default
     27! - restructure the average arrays for gases and total mass concentrations of
     28!   chemical components: set to 4d arrays instead of separate arrays
     29! - add allocation checks for averaged data output arrays
    2430!
    2531! Former revisions:
     
    499505    LOGICAL ::  nesting_offline_salsa   = .TRUE.   !< Apply offline nesting for salsa
    500506    LOGICAL ::  no_insoluble            = .FALSE.  !< Exclude insoluble chemical components
    501     LOGICAL ::  read_restart_data_salsa = .FALSE.  !< Read restart data for salsa
     507    LOGICAL ::  read_restart_data_salsa = .TRUE.   !< Read restart data for salsa
    502508    LOGICAL ::  salsa_gases_from_chem   = .FALSE.  !< Transfer the gaseous components to SALSA
    503509    LOGICAL ::  van_der_waals_coagc     = .FALSE.  !< Include van der Waals and viscous forces in coagulation
    504     LOGICAL ::  write_binary_salsa      = .FALSE.  !< read binary for salsa
     510    LOGICAL ::  write_binary_salsa      = .TRUE.   !< read binary for salsa
    505511!
    506512!-- Process switches: nl* is read from the NAMELIST and is NOT changed.
     
    827833!-- Data output arrays:
    828834!
    829 !-- Gases:
    830     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_h2so4_av  !< H2SO4
    831     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_hno3_av   !< HNO3
    832     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_nh3_av    !< NH3
    833     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_ocnv_av   !< non-volatile OC
    834     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_ocsv_av   !< semi-volatile OC
    835 !
    836835!-- Integrated:
    837836    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ldsa_av  !< lung-deposited surface area
     
    842841    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  pm10_av  !< PM10
    843842!
    844 !-- In the particle phase:
    845     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_bc_av   !< black carbon
    846     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_du_av   !< dust
    847     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_h2o_av  !< liquid water
    848     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_nh_av   !< ammonia
    849     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_no_av   !< nitrates
    850     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_oc_av   !< org. carbon
    851     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_so4_av  !< sulphates
    852     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_ss_av   !< sea salt
    853 !
    854843!-- Bin specific mass and number concentrations:
    855844    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  mbins_av  !< bin mas
    856845    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  nbins_av  !< bin number
     846!
     847!-- Gases:
     848    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  salsa_gases_av  !< gases
     849!
     850!-- In the particle phase:
     851    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_h2o_av  !< liquid water
     852    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  s_mass_av  !< mass components
    857853
    858854!
     
    12171213       init_gases_type = 1
    12181214    ENDIF
    1219 
     1215!
     1216!-- If the run is not a restart run, set read_restart_data to .FALSE.
     1217    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
     1218       read_restart_data_salsa = .FALSE.
     1219    ENDIF
    12201220
    12211221 END SUBROUTINE salsa_check_parameters
     
    29072907       SELECT CASE ( restart_string(1:length) )
    29082908
     2909          CASE ( 'aerosol_mass' )
     2910             DO  ic = 1, ncomponents_mass * nbins_aerosol
     2911                IF ( k == 1 )  READ ( 13 ) tmp_3d
     2912                aerosol_mass(ic)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                 &
     2913                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     2914             ENDDO
     2915             found = .TRUE.
     2916
    29092917          CASE ( 'aerosol_number' )
    29102918             DO  ib = 1, nbins_aerosol
     
    29122920                aerosol_number(ib)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =               &
    29132921                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    2914                 found = .TRUE.
    29152922             ENDDO
    2916 
    2917           CASE ( 'aerosol_mass' )
    2918              DO  ic = 1, ncomponents_mass * nbins_aerosol
     2923             found = .TRUE.
     2924
     2925          CASE( 'salsa_gases_av' )
     2926             IF ( .NOT. ALLOCATED( salsa_gases_av ) )  THEN
     2927                ALLOCATE( salsa_gases_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
     2928             ENDIF
     2929             DO  ig = 1, ngases_salsa
    29192930                IF ( k == 1 )  READ ( 13 ) tmp_3d
    2920                 aerosol_mass(ic)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                 &
     2931                salsa_gases_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp,ig) =                     &
     2932                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     2933             ENDDO
     2934             found = .TRUE.
     2935
     2936          CASE ( 'ldsa_av' )
     2937             IF ( .NOT. ALLOCATED( ldsa_av ) )  ALLOCATE( ldsa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2938             IF ( k == 1 )  READ ( 13 ) tmp_3d
     2939             ldsa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
     2940                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     2941             found = .TRUE.
     2942
     2943          CASE ( 'mbins_av' )
     2944             IF ( .NOT. ALLOCATED( mbins_av ) )  THEN
     2945                ALLOCATE( mbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
     2946             ENDIF
     2947             DO  ib = 1, nbins_aerosol
     2948                IF ( k == 1 )  READ ( 13 ) tmp_3d
     2949                mbins_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp,ib) =                           &
    29212950                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    29222951                found = .TRUE.
    29232952             ENDDO
    29242953
    2925           CASE ( 'salsa_gas' )
    2926              DO  ig = 1, ngases_salsa
     2954          CASE ( 'nbins_av' )
     2955             IF ( .NOT. ALLOCATED( nbins_av ) )  THEN
     2956                ALLOCATE( nbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
     2957             ENDIF
     2958             DO  ib = 1, nbins_aerosol
    29272959                IF ( k == 1 )  READ ( 13 ) tmp_3d
    2928                 salsa_gas(ig)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                    &
     2960                nbins_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp,ib) =                           &
    29292961                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    29302962                found = .TRUE.
    29312963             ENDDO
     2964
     2965          CASE ( 'ntot_av' )
     2966             IF ( .NOT. ALLOCATED( ntot_av ) )  ALLOCATE( ntot_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2967             IF ( k == 1 )  READ ( 13 ) tmp_3d
     2968             ntot_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
     2969                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     2970             found = .TRUE.
     2971
     2972          CASE ( 'nufp_av' )
     2973             IF ( .NOT. ALLOCATED( nufp_av ) )  ALLOCATE( nufp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2974             IF ( k == 1 )  READ ( 13 ) tmp_3d
     2975             nufp_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
     2976                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     2977             found = .TRUE.
     2978
     2979          CASE ( 'pm01_av' )
     2980             IF ( .NOT. ALLOCATED( pm01_av ) )  ALLOCATE( pm01_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2981             IF ( k == 1 )  READ ( 13 ) tmp_3d
     2982             pm01_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
     2983                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     2984             found = .TRUE.
     2985
     2986          CASE ( 'pm25_av' )
     2987             IF ( .NOT. ALLOCATED( pm25_av ) )  ALLOCATE( pm25_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2988             IF ( k == 1 )  READ ( 13 ) tmp_3d
     2989             pm25_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
     2990                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     2991             found = .TRUE.
     2992
     2993          CASE ( 'pm10_av' )
     2994             IF ( .NOT. ALLOCATED( pm10_av ) )  ALLOCATE( pm10_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2995             IF ( k == 1 )  READ ( 13 ) tmp_3d
     2996             pm10_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
     2997                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     2998             found = .TRUE.
     2999
     3000          CASE ( 's_mass_av' )
     3001             IF ( .NOT. ALLOCATED( s_mass_av ) )  THEN
     3002                ALLOCATE( s_mass_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass) )
     3003             ENDIF
     3004             DO  ic = 1, ncomponents_mass
     3005                IF ( k == 1 )  READ ( 13 ) tmp_3d
     3006                s_mass_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp,ic) =                          &
     3007                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3008             ENDDO
     3009             found = .TRUE.
     3010
     3011          CASE ( 's_h2o_av' )
     3012             IF ( .NOT. ALLOCATED( s_h2o_av ) )  ALLOCATE( s_h2o_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3013             IF ( k == 1 )  READ ( 13 ) tmp_3d
     3014             s_h2o_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                 &
     3015                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3016             found = .TRUE.
     3017
     3018          CASE ( 'salsa_gas' )
     3019             IF ( .NOT. salsa_gases_from_chem )  THEN
     3020                DO  ig = 1, ngases_salsa
     3021                   IF ( k == 1 )  READ ( 13 ) tmp_3d
     3022                   salsa_gas(ig)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                 &
     3023                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3024                ENDDO
     3025                found = .TRUE.
     3026             ENDIF
    29323027
    29333028          CASE DEFAULT
     
    29603055    IF ( write_binary  .AND.  write_binary_salsa )  THEN
    29613056
     3057       CALL wrd_write_string( 'aerosol_mass' )
     3058       DO  ic = 1, nbins_aerosol * ncomponents_mass
     3059          WRITE ( 14 )  aerosol_mass(ic)%conc
     3060       ENDDO
     3061
    29623062       CALL wrd_write_string( 'aerosol_number' )
    29633063       DO  ib = 1, nbins_aerosol
     
    29653065       ENDDO
    29663066
    2967        CALL wrd_write_string( 'aerosol_mass' )
    2968        DO  ic = 1, nbins_aerosol * ncomponents_mass
    2969           WRITE ( 14 )  aerosol_mass(ic)%conc
    2970        ENDDO
    2971 
    2972        CALL wrd_write_string( 'salsa_gas' )
    2973        DO  ig = 1, ngases_salsa
    2974           WRITE ( 14 )  salsa_gas(ig)%conc
    2975        ENDDO
     3067       IF (  .NOT. salsa_gases_from_chem )  THEN
     3068
     3069          IF ( ALLOCATED( salsa_gases_av ) )  THEN
     3070             CALL wrd_write_string( 'salsa_gases_av' )
     3071             DO  ig = 1, ngases_salsa
     3072                WRITE ( 14 )  salsa_gases_av(:,:,:,ig)
     3073             ENDDO
     3074          ENDIF
     3075       ENDIF
     3076
     3077       IF ( ALLOCATED( ldsa_av ) )  THEN
     3078       CALL wrd_write_string( 'ldsa_av' )
     3079       WRITE ( 14 )  ldsa_av
     3080       ENDIF
     3081
     3082       IF ( ALLOCATED( mbins_av ) )  THEN
     3083          CALL wrd_write_string( 'mbins_av' )
     3084          DO  ib = 1, nbins_aerosol
     3085             WRITE ( 14 )  mbins_av(:,:,:,ib)
     3086          ENDDO
     3087       ENDIF
     3088
     3089       IF ( ALLOCATED( nbins_av ) )  THEN
     3090          CALL wrd_write_string( 'nbins_av' )
     3091          DO  ib = 1, nbins_aerosol
     3092             WRITE ( 14 )  nbins_av(:,:,:,ib)
     3093          ENDDO
     3094       ENDIF
     3095
     3096       IF ( ALLOCATED( ldsa_av ) )  THEN
     3097          CALL wrd_write_string( 'ntot_av' )
     3098          WRITE ( 14 )  ntot_av
     3099       ENDIF
     3100
     3101       IF ( ALLOCATED( nufp_av ) )  THEN
     3102          CALL wrd_write_string( 'nufp_av' )
     3103          WRITE ( 14 )  nufp_av
     3104       ENDIF
     3105
     3106       IF ( ALLOCATED( pm01_av ) )  THEN
     3107          CALL wrd_write_string( 'pm01_av' )
     3108          WRITE ( 14 )  pm01_av
     3109       ENDIF
     3110
     3111       IF ( ALLOCATED( pm25_av ) )  THEN
     3112          CALL wrd_write_string( 'pm25_av' )
     3113          WRITE ( 14 )  pm25_av
     3114       ENDIF
     3115
     3116       IF ( ALLOCATED( pm10_av ) )  THEN
     3117          CALL wrd_write_string( 'pm10_av' )
     3118          WRITE ( 14 )  pm10_av
     3119       ENDIF
     3120
     3121       IF ( ALLOCATED( s_mass_av ) )  THEN
     3122          CALL wrd_write_string( 's_mass_av' )
     3123          DO  ic = 1, ncomponents_mass
     3124             WRITE ( 14 )  s_mass_av(:,:,:,ic)
     3125          ENDDO
     3126       ENDIF
     3127
     3128       IF ( ALLOCATED( s_h2o_av ) )  THEN
     3129          CALL wrd_write_string( 's_h2o_av' )
     3130          WRITE ( 14 )  s_h2o_av
     3131       ENDIF
     3132
     3133       IF ( .NOT. salsa_gases_from_chem )  THEN
     3134          CALL wrd_write_string( 'salsa_gas' )
     3135          DO  ig = 1, ngases_salsa
     3136             WRITE ( 14 )  salsa_gas(ig)%conc
     3137          ENDDO
     3138       ENDIF
    29763139
    29773140    ENDIF
     
    33783541             ic = ic+1
    33793542          ENDDO
    3380        ENDIF
    3381        IF ( prunmode == 1 )  THEN
    3382           nc_h2o = get_index( prtcl,'H2O' )
    3383           vc = 8
    3384           str = ( nc_h2o-1 ) * nbins_aerosol + 1
    3385           endi = nc_h2o * nbins_aerosol
    3386           ic = 1
    3387           DO ss = str, endi
    3388              aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), ( lo_aero(ic)%volc(vc) - &
    3389                                              aero_old(ic)%volc(vc) ) * arhoh2o )
    3390              IF ( k == nzb+1 )  THEN
    3391                 aerosol_mass(ss)%init(k-1) = aerosol_mass(ss)%init(k)
    3392              ELSEIF ( k == nzt  )  THEN
    3393                 aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k)
    3394                 aerosol_mass(ss)%conc(k+1,j,i) = aerosol_mass(ss)%init(k)
    3395              ENDIF
    3396              ic = ic+1
    3397           ENDDO
     3543          IF ( prunmode == 1 )  THEN
     3544             nc_h2o = get_index( prtcl,'H2O' )
     3545             vc = 8
     3546             str = ( nc_h2o-1 ) * nbins_aerosol + 1
     3547             endi = nc_h2o * nbins_aerosol
     3548             ic = 1
     3549             DO ss = str, endi
     3550                aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), ( lo_aero(ic)%volc(vc) - &
     3551                                                aero_old(ic)%volc(vc) ) * arhoh2o )
     3552                IF ( k == nzb+1 )  THEN
     3553                   aerosol_mass(ss)%init(k-1) = aerosol_mass(ss)%init(k)
     3554                ELSEIF ( k == nzt  )  THEN
     3555                   aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k)
     3556                   aerosol_mass(ss)%conc(k+1,j,i) = aerosol_mass(ss)%init(k)
     3557                ENDIF
     3558                ic = ic+1
     3559             ENDDO
     3560          ENDIF
    33983561       ENDIF
    33993562!
     
    98339996       ENDIF
    98349997
     9998    ELSEIF ( var(7:11) == 's_H2O' )  THEN
     9999       IF ( .NOT. advect_particle_water )  THEN
     10000          message_string = 'to output s_H2O/s_H2O_av requires that advect_particle_water = .T.'
     10001          CALL message( 'check_parameters', 'PA0707', 1, 2, 0, 6, 0 )
     10002       ENDIF
     10003
    983510004    ELSE
    983610005       SELECT CASE ( TRIM( var(7:) ) )
     
    1016910338          SELECT CASE ( TRIM( variable(7:) ) )
    1017010339
    10171              CASE ( 'g_H2SO4' )
    10172                 IF ( .NOT. ALLOCATED( g_h2so4_av ) )  THEN
    10173                    ALLOCATE( g_h2so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     10340             CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
     10341                IF ( .NOT. ALLOCATED( salsa_gases_av ) )  THEN
     10342                   ALLOCATE( salsa_gases_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
    1017410343                ENDIF
    10175                 g_h2so4_av = 0.0_wp
    10176 
    10177              CASE ( 'g_HNO3' )
    10178                 IF ( .NOT. ALLOCATED( g_hno3_av ) )  THEN
    10179                    ALLOCATE( g_hno3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10180                 ENDIF
    10181                 g_hno3_av = 0.0_wp
    10182 
    10183              CASE ( 'g_NH3' )
    10184                 IF ( .NOT. ALLOCATED( g_nh3_av ) )  THEN
    10185                    ALLOCATE( g_nh3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10186                 ENDIF
    10187                 g_nh3_av = 0.0_wp
    10188 
    10189              CASE ( 'g_OCNV' )
    10190                 IF ( .NOT. ALLOCATED( g_ocnv_av ) )  THEN
    10191                    ALLOCATE( g_ocnv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10192                 ENDIF
    10193                 g_ocnv_av = 0.0_wp
    10194 
    10195              CASE ( 'g_OCSV' )
    10196                 IF ( .NOT. ALLOCATED( g_ocsv_av ) )  THEN
    10197                    ALLOCATE( g_ocsv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10198                 ENDIF
    10199                 g_ocsv_av = 0.0_wp
     10344                salsa_gases_av = 0.0_wp
    1020010345
    1020110346             CASE ( 'LDSA' )
     
    1023510380                pm10_av = 0.0_wp
    1023610381
    10237              CASE ( 's_BC' )
    10238                 IF ( .NOT. ALLOCATED( s_bc_av ) )  THEN
    10239                    ALLOCATE( s_bc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     10382             CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
     10383                IF ( .NOT. ALLOCATED( s_mass_av ) )  THEN
     10384                   ALLOCATE( s_mass_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass+1) )
    1024010385                ENDIF
    10241                 s_bc_av = 0.0_wp
    10242 
    10243              CASE ( 's_DU' )
    10244                 IF ( .NOT. ALLOCATED( s_du_av ) )  THEN
    10245                    ALLOCATE( s_du_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10246                 ENDIF
    10247                 s_du_av = 0.0_wp
     10386                s_mass_av = 0.0_wp
    1024810387
    1024910388             CASE ( 's_H2O' )
     
    1025310392                s_h2o_av = 0.0_wp
    1025410393
    10255              CASE ( 's_NH' )
    10256                 IF ( .NOT. ALLOCATED( s_nh_av ) )  THEN
    10257                    ALLOCATE( s_nh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10258                 ENDIF
    10259                 s_nh_av = 0.0_wp
    10260 
    10261              CASE ( 's_NO' )
    10262                 IF ( .NOT. ALLOCATED( s_no_av ) )  THEN
    10263                    ALLOCATE( s_no_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10264                 ENDIF
    10265                 s_no_av = 0.0_wp
    10266 
    10267              CASE ( 's_OC' )
    10268                 IF ( .NOT. ALLOCATED( s_oc_av ) )  THEN
    10269                    ALLOCATE( s_oc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10270                 ENDIF
    10271                 s_oc_av = 0.0_wp
    10272 
    10273              CASE ( 's_SO4' )
    10274                 IF ( .NOT. ALLOCATED( s_so4_av ) )  THEN
    10275                    ALLOCATE( s_so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10276                 ENDIF
    10277                 s_so4_av = 0.0_wp
    10278 
    10279              CASE ( 's_SS' )
    10280                 IF ( .NOT. ALLOCATED( s_ss_av ) )  THEN
    10281                    ALLOCATE( s_ss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10282                 ENDIF
    10283                 s_ss_av = 0.0_wp
    10284 
    1028510394             CASE DEFAULT
    1028610395                CONTINUE
     
    1029310402
    1029410403       IF ( variable(7:11) ==  'N_bin' )  THEN
    10295           READ( variable(12:),* ) char_to_int
    10296           IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
    10297              ib = char_to_int
    10298              DO  i = nxlg, nxrg
    10299                 DO  j = nysg, nyng
    10300                    DO  k = nzb, nzt+1
    10301                       nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) + aerosol_number(ib)%conc(k,j,i)
    10302                    ENDDO
    10303                 ENDDO
    10304              ENDDO
    10305           ENDIF
    10306 
    10307        ELSEIF ( variable(7:11) ==  'm_bin' )  THEN
    10308           READ( variable(12:),* ) char_to_int
    10309           IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
    10310              ib = char_to_int
    10311              DO  i = nxlg, nxrg
    10312                 DO  j = nysg, nyng
    10313                    DO  k = nzb, nzt+1
    10314                       temp_bin = 0.0_wp
    10315                       DO  ic = ib, nbins_aerosol * ncomponents_mass, nbins_aerosol
    10316                          temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    10317                       ENDDO
    10318                       mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) + temp_bin
    10319                    ENDDO
    10320                 ENDDO
    10321              ENDDO
    10322           ENDIF
    10323        ELSE
    10324 
    10325           SELECT CASE ( TRIM( variable(7:) ) )
    10326 
    10327              CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
    10328 
    10329                 vari = TRIM( variable(9:) )  ! remove salsa_g_ from beginning
    10330 
    10331                 SELECT CASE( vari )
    10332 
    10333                    CASE( 'H2SO4' )
    10334                       found_index = 1
    10335                       to_be_resorted => g_h2so4_av
    10336 
    10337                    CASE( 'HNO3' )
    10338                       found_index = 2
    10339                       to_be_resorted => g_hno3_av
    10340 
    10341                    CASE( 'NH3' )
    10342                       found_index = 3
    10343                       to_be_resorted => g_nh3_av
    10344 
    10345                    CASE( 'OCNV' )
    10346                       found_index = 4
    10347                       to_be_resorted => g_ocnv_av
    10348 
    10349                    CASE( 'OCSV' )
    10350                       found_index = 5
    10351                       to_be_resorted => g_ocsv_av
    10352 
    10353                 END SELECT
    10354 
     10404          IF ( ALLOCATED( nbins_av ) )  THEN
     10405             READ( variable(12:),* ) char_to_int
     10406             IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     10407                ib = char_to_int
    1035510408                DO  i = nxlg, nxrg
    1035610409                   DO  j = nysg, nyng
    1035710410                      DO  k = nzb, nzt+1
    10358                          to_be_resorted(k,j,i) = to_be_resorted(k,j,i) +                           &
    10359                                                  salsa_gas(found_index)%conc(k,j,i)
     10411                         nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) + aerosol_number(ib)%conc(k,j,i)
    1036010412                      ENDDO
    1036110413                   ENDDO
    1036210414                ENDDO
    10363 
    10364              CASE ( 'LDSA' )
     10415             ENDIF
     10416          ENDIF
     10417
     10418       ELSEIF ( variable(7:11) ==  'm_bin' )  THEN
     10419          IF ( ALLOCATED( mbins_av ) )  THEN
     10420             READ( variable(12:),* ) char_to_int
     10421             IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     10422                ib = char_to_int
    1036510423                DO  i = nxlg, nxrg
    1036610424                   DO  j = nysg, nyng
    1036710425                      DO  k = nzb, nzt+1
    1036810426                         temp_bin = 0.0_wp
    10369                          DO  ib = 1, nbins_aerosol
    10370    !
    10371    !--                      Diameter in micrometres
    10372                             mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
    10373    !
    10374    !--                      Deposition factor: alveolar (use ra_dry)
    10375                             df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
    10376                                    2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
    10377                                    1.362_wp )**2 ) )
    10378    !
    10379    !--                      Lung-deposited surface area LDSA (units mum2/cm3)
    10380                             temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
    10381                                        aerosol_number(ib)%conc(k,j,i)
     10427                         DO  ic = ib, nbins_aerosol * ncomponents_mass, nbins_aerosol
     10428                            temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    1038210429                         ENDDO
    10383                          ldsa_av(k,j,i) = ldsa_av(k,j,i) + temp_bin
     10430                         mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) + temp_bin
    1038410431                      ENDDO
    1038510432                   ENDDO
    1038610433                ENDDO
    10387 
    10388              CASE ( 'N_UFP' )
    10389                 DO  i = nxlg, nxrg
    10390                    DO  j = nysg, nyng
    10391                       DO  k = nzb, nzt+1
    10392                          temp_bin = 0.0_wp
    10393                          DO  ib = 1, nbins_aerosol
    10394                             IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
    10395                                temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
    10396                             ENDIF
    10397                          ENDDO
    10398                          nufp_av(k,j,i) = nufp_av(k,j,i) + temp_bin
    10399                       ENDDO
    10400                    ENDDO
    10401                 ENDDO
    10402 
    10403              CASE ( 'Ntot' )
    10404                 DO  i = nxlg, nxrg
    10405                    DO  j = nysg, nyng
    10406                       DO  k = nzb, nzt+1
    10407                          DO  ib = 1, nbins_aerosol
    10408                             ntot_av(k,j,i) = ntot_av(k,j,i) + aerosol_number(ib)%conc(k,j,i)
     10434             ENDIF
     10435          ENDIF
     10436       ELSE
     10437
     10438          SELECT CASE ( TRIM( variable(7:) ) )
     10439
     10440             CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
     10441                IF ( ALLOCATED( salsa_gases_av ) )  THEN
     10442
     10443                   vari = TRIM( variable(9:) )  ! remove salsa_g_ from beginning
     10444
     10445                   SELECT CASE( vari )
     10446
     10447                      CASE( 'H2SO4' )
     10448                         found_index = 1
     10449                      CASE( 'HNO3' )
     10450                         found_index = 2
     10451                      CASE( 'NH3' )
     10452                         found_index = 3
     10453                      CASE( 'OCNV' )
     10454                         found_index = 4
     10455                      CASE( 'OCSV' )
     10456                         found_index = 5
     10457
     10458                   END SELECT
     10459
     10460                   DO  i = nxlg, nxrg
     10461                      DO  j = nysg, nyng
     10462                         DO  k = nzb, nzt+1
     10463                            salsa_gases_av(k,j,i,found_index) = salsa_gases_av(k,j,i,found_index)  &
     10464                                                                + salsa_gas(found_index)%conc(k,j,i)
    1040910465                         ENDDO
    1041010466                      ENDDO
    1041110467                   ENDDO
    10412                 ENDDO
    10413 
    10414              CASE ( 'PM0.1' )
    10415                 DO  i = nxlg, nxrg
    10416                    DO  j = nysg, nyng
    10417                       DO  k = nzb, nzt+1
    10418                          temp_bin = 0.0_wp
    10419                          DO  ib = 1, nbins_aerosol
    10420                             IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
    10421                                DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
    10422                                   temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    10423                                ENDDO
    10424                             ENDIF
    10425                          ENDDO
    10426                          pm01_av(k,j,i) = pm01_av(k,j,i) + temp_bin
    10427                       ENDDO
    10428                    ENDDO
    10429                 ENDDO
    10430 
    10431              CASE ( 'PM2.5' )
    10432                 DO  i = nxlg, nxrg
    10433                    DO  j = nysg, nyng
    10434                       DO  k = nzb, nzt+1
    10435                          temp_bin = 0.0_wp
    10436                          DO  ib = 1, nbins_aerosol
    10437                             IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp )  THEN
    10438                                DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
    10439                                   temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    10440                                ENDDO
    10441                             ENDIF
    10442                          ENDDO
    10443                          pm25_av(k,j,i) = pm25_av(k,j,i) + temp_bin
    10444                       ENDDO
    10445                    ENDDO
    10446                 ENDDO
    10447 
    10448              CASE ( 'PM10' )
    10449                 DO  i = nxlg, nxrg
    10450                    DO  j = nysg, nyng
    10451                       DO  k = nzb, nzt+1
    10452                          temp_bin = 0.0_wp
    10453                          DO  ib = 1, nbins_aerosol
    10454                             IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp )  THEN
    10455                                DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
    10456                                   temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    10457                                ENDDO
    10458                             ENDIF
    10459                          ENDDO
    10460                          pm10_av(k,j,i) = pm10_av(k,j,i) + temp_bin
    10461                       ENDDO
    10462                    ENDDO
    10463                 ENDDO
    10464 
    10465              CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
    10466                 IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN  ! 9: remove salsa_s_
    10467                    found_index = get_index( prtcl, TRIM( variable(9:) ) )
    10468                    IF ( TRIM( variable(9:) ) == 'BC' )   to_be_resorted => s_bc_av
    10469                    IF ( TRIM( variable(9:) ) == 'DU' )   to_be_resorted => s_du_av
    10470                    IF ( TRIM( variable(9:) ) == 'NH' )   to_be_resorted => s_nh_av
    10471                    IF ( TRIM( variable(9:) ) == 'NO' )   to_be_resorted => s_no_av
    10472                    IF ( TRIM( variable(9:) ) == 'OC' )   to_be_resorted => s_oc_av
    10473                    IF ( TRIM( variable(9:) ) == 'SO4' )  to_be_resorted => s_so4_av
    10474                    IF ( TRIM( variable(9:) ) == 'SS' )   to_be_resorted => s_ss_av
     10468                ENDIF
     10469
     10470             CASE ( 'LDSA' )
     10471                IF ( ALLOCATED( ldsa_av ) )  THEN
    1047510472                   DO  i = nxlg, nxrg
    1047610473                      DO  j = nysg, nyng
    1047710474                         DO  k = nzb, nzt+1
    10478                             DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
    10479                                to_be_resorted(k,j,i) = to_be_resorted(k,j,i) +                     &
    10480                                                        aerosol_mass(ic)%conc(k,j,i)
     10475                            temp_bin = 0.0_wp
     10476                            DO  ib = 1, nbins_aerosol
     10477!
     10478!--                            Diameter in micrometres
     10479                               mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
     10480!
     10481!--                            Deposition factor: alveolar (use ra_dry)
     10482                               df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) + &
     10483                                      2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) &
     10484                                      - 1.362_wp )**2 ) )
     10485!
     10486!--                            Lung-deposited surface area LDSA (units mum2/cm3)
     10487                               temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *             &
     10488                                          aerosol_number(ib)%conc(k,j,i)
     10489                            ENDDO
     10490                            ldsa_av(k,j,i) = ldsa_av(k,j,i) + temp_bin
     10491                         ENDDO
     10492                      ENDDO
     10493                   ENDDO
     10494                ENDIF
     10495
     10496             CASE ( 'N_UFP' )
     10497                IF ( ALLOCATED( nufp_av ) )  THEN
     10498                   DO  i = nxlg, nxrg
     10499                      DO  j = nysg, nyng
     10500                         DO  k = nzb, nzt+1
     10501                            temp_bin = 0.0_wp
     10502                            DO  ib = 1, nbins_aerosol
     10503                               IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
     10504                                  temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
     10505                               ENDIF
     10506                            ENDDO
     10507                            nufp_av(k,j,i) = nufp_av(k,j,i) + temp_bin
     10508                         ENDDO
     10509                      ENDDO
     10510                   ENDDO
     10511                ENDIF
     10512
     10513             CASE ( 'Ntot' )
     10514               IF ( ALLOCATED( ntot_av ) )  THEN
     10515                   DO  i = nxlg, nxrg
     10516                      DO  j = nysg, nyng
     10517                         DO  k = nzb, nzt+1
     10518                            DO  ib = 1, nbins_aerosol
     10519                               ntot_av(k,j,i) = ntot_av(k,j,i) + aerosol_number(ib)%conc(k,j,i)
    1048110520                            ENDDO
    1048210521                         ENDDO
     
    1048510524                ENDIF
    1048610525
     10526             CASE ( 'PM0.1' )
     10527                IF ( ALLOCATED( pm01_av ) )  THEN
     10528                   DO  i = nxlg, nxrg
     10529                      DO  j = nysg, nyng
     10530                         DO  k = nzb, nzt+1
     10531                            temp_bin = 0.0_wp
     10532                            DO  ib = 1, nbins_aerosol
     10533                               IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
     10534                                  DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     10535                                     temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     10536                                  ENDDO
     10537                               ENDIF
     10538                            ENDDO
     10539                            pm01_av(k,j,i) = pm01_av(k,j,i) + temp_bin
     10540                         ENDDO
     10541                      ENDDO
     10542                   ENDDO
     10543                ENDIF
     10544
     10545             CASE ( 'PM2.5' )
     10546                IF ( ALLOCATED( pm25_av ) )  THEN
     10547                   DO  i = nxlg, nxrg
     10548                      DO  j = nysg, nyng
     10549                         DO  k = nzb, nzt+1
     10550                            temp_bin = 0.0_wp
     10551                            DO  ib = 1, nbins_aerosol
     10552                               IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp )  THEN
     10553                                  DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     10554                                     temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     10555                                  ENDDO
     10556                               ENDIF
     10557                            ENDDO
     10558                            pm25_av(k,j,i) = pm25_av(k,j,i) + temp_bin
     10559                         ENDDO
     10560                      ENDDO
     10561                   ENDDO
     10562                ENDIF
     10563
     10564             CASE ( 'PM10' )
     10565                IF ( ALLOCATED( pm10_av ) )  THEN
     10566                   DO  i = nxlg, nxrg
     10567                      DO  j = nysg, nyng
     10568                         DO  k = nzb, nzt+1
     10569                            temp_bin = 0.0_wp
     10570                            DO  ib = 1, nbins_aerosol
     10571                               IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp )  THEN
     10572                                  DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     10573                                     temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     10574                                  ENDDO
     10575                               ENDIF
     10576                            ENDDO
     10577                            pm10_av(k,j,i) = pm10_av(k,j,i) + temp_bin
     10578                         ENDDO
     10579                      ENDDO
     10580                   ENDDO
     10581                ENDIF
     10582
     10583             CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
     10584                IF ( ALLOCATED( s_mass_av ) )  THEN
     10585                   IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN  ! 9: remove salsa_s_
     10586                      found_index = get_index( prtcl, TRIM( variable(9:) ) )
     10587                      DO  i = nxlg, nxrg
     10588                         DO  j = nysg, nyng
     10589                            DO  k = nzb, nzt+1
     10590                               DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
     10591                                  s_mass_av(k,j,i,found_index) = s_mass_av(k,j,i,found_index) +    &
     10592                                                                 aerosol_mass(ic)%conc(k,j,i)
     10593                               ENDDO
     10594                            ENDDO
     10595                         ENDDO
     10596                      ENDDO
     10597                   ENDIF
     10598                ENDIF
     10599
    1048710600             CASE ( 's_H2O' )
    10488                 found_index = get_index( prtcl,'H2O' )
    10489                 to_be_resorted => s_h2o_av
     10601                IF ( ALLOCATED( s_H2O_av ) )  THEN
     10602                   found_index = get_index( prtcl,'H2O' )
     10603                   to_be_resorted => s_h2o_av
     10604                   DO  i = nxlg, nxrg
     10605                      DO  j = nysg, nyng
     10606                         DO  k = nzb, nzt+1
     10607                            DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
     10608                               s_h2o_av(k,j,i) = s_h2o_av(k,j,i) + aerosol_mass(ic)%conc(k,j,i)
     10609                            ENDDO
     10610                         ENDDO
     10611                      ENDDO
     10612                   ENDDO
     10613                ENDIF
     10614
     10615             CASE DEFAULT
     10616                CONTINUE
     10617
     10618          END SELECT
     10619
     10620       ENDIF
     10621
     10622    ELSEIF ( mode == 'average' )  THEN
     10623
     10624       IF ( variable(7:11) ==  'N_bin' )  THEN
     10625          IF ( ALLOCATED( nbins_av ) )  THEN
     10626             READ( variable(12:),* ) char_to_int
     10627             IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     10628                ib = char_to_int
    1049010629                DO  i = nxlg, nxrg
    1049110630                   DO  j = nysg, nyng
    1049210631                      DO  k = nzb, nzt+1
    10493                          DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
    10494                             s_h2o_av(k,j,i) = s_h2o_av(k,j,i) + aerosol_mass(ic)%conc(k,j,i)
     10632                         nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp )
     10633                      ENDDO
     10634                   ENDDO
     10635                ENDDO
     10636             ENDIF
     10637          ENDIF
     10638
     10639       ELSEIF ( variable(7:11) ==  'm_bin' )  THEN
     10640          IF ( ALLOCATED( mbins_av ) )  THEN
     10641             READ( variable(12:),* ) char_to_int
     10642             IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     10643                ib = char_to_int
     10644                DO  i = nxlg, nxrg
     10645                   DO  j = nysg, nyng
     10646                      DO  k = nzb, nzt+1
     10647                         mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp)
     10648                      ENDDO
     10649                   ENDDO
     10650                ENDDO
     10651             ENDIF
     10652          ENDIF
     10653       ELSE
     10654
     10655          SELECT CASE ( TRIM( variable(7:) ) )
     10656
     10657             CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
     10658                IF ( ALLOCATED( salsa_gases_av ) )  THEN
     10659                   IF ( TRIM( variable(9:) ) == 'H2SO4' )  THEN  ! 9: remove salsa_g_ from beginning
     10660                      found_index = 1
     10661                   ELSEIF ( TRIM( variable(9:) ) == 'HNO3' )  THEN
     10662                      found_index = 2
     10663                   ELSEIF ( TRIM( variable(9:) ) == 'NH3' )  THEN
     10664                      found_index = 3
     10665                   ELSEIF ( TRIM( variable(9:) ) == 'OCNV' )  THEN
     10666                      found_index = 4
     10667                   ELSEIF ( TRIM( variable(9:) ) == 'OCSV' )  THEN
     10668                      found_index = 5
     10669                   ENDIF
     10670                   DO  i = nxlg, nxrg
     10671                      DO  j = nysg, nyng
     10672                         DO  k = nzb, nzt+1
     10673                            salsa_gases_av(k,j,i,found_index) = salsa_gases_av(k,j,i,found_index)  &
     10674                                                                / REAL( average_count_3d, KIND=wp )
    1049510675                         ENDDO
    1049610676                      ENDDO
    1049710677                   ENDDO
    10498                 ENDDO
    10499 
    10500              CASE DEFAULT
    10501                 CONTINUE
    10502 
    10503           END SELECT
    10504 
    10505        ENDIF
    10506 
    10507     ELSEIF ( mode == 'average' )  THEN
    10508 
    10509        IF ( variable(7:11) ==  'N_bin' )  THEN
    10510           READ( variable(12:),* ) char_to_int
    10511           IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
    10512              ib = char_to_int
    10513              DO  i = nxlg, nxrg
    10514                 DO  j = nysg, nyng
    10515                    DO  k = nzb, nzt+1
    10516                       nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp )
    10517                    ENDDO
    10518                 ENDDO
    10519              ENDDO
    10520           ENDIF
    10521 
    10522        ELSEIF ( variable(7:11) ==  'm_bin' )  THEN
    10523           READ( variable(12:),* ) char_to_int
    10524           IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
    10525              ib = char_to_int
    10526              DO  i = nxlg, nxrg
    10527                 DO  j = nysg, nyng
    10528                    DO  k = nzb, nzt+1
    10529                       mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp)
    10530                    ENDDO
    10531                 ENDDO
    10532              ENDDO
    10533           ENDIF
    10534        ELSE
    10535 
    10536           SELECT CASE ( TRIM( variable(7:) ) )
    10537 
    10538              CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
    10539                 IF ( TRIM( variable(9:) ) == 'H2SO4' )  THEN  ! 9: remove salsa_g_ from beginning
    10540                    found_index = 1
    10541                    to_be_resorted => g_h2so4_av
    10542                 ELSEIF ( TRIM( variable(9:) ) == 'HNO3' )  THEN
    10543                    found_index = 2
    10544                    to_be_resorted => g_hno3_av
    10545                 ELSEIF ( TRIM( variable(9:) ) == 'NH3' )  THEN
    10546                    found_index = 3
    10547                    to_be_resorted => g_nh3_av
    10548                 ELSEIF ( TRIM( variable(9:) ) == 'OCNV' )  THEN
    10549                    found_index = 4
    10550                    to_be_resorted => g_ocnv_av
    10551                 ELSEIF ( TRIM( variable(9:) ) == 'OCSV' )  THEN
    10552                    found_index = 5
    10553                    to_be_resorted => g_ocsv_av
    1055410678                ENDIF
    10555                 DO  i = nxlg, nxrg
    10556                    DO  j = nysg, nyng
    10557                       DO  k = nzb, nzt+1
    10558                          to_be_resorted(k,j,i) = to_be_resorted(k,j,i) /                           &
    10559                                                  REAL( average_count_3d, KIND=wp )
    10560                       ENDDO
    10561                    ENDDO
    10562                 ENDDO
    1056310679
    1056410680             CASE ( 'LDSA' )
    10565                 DO  i = nxlg, nxrg
    10566                    DO  j = nysg, nyng
    10567                       DO  k = nzb, nzt+1
    10568                          ldsa_av(k,j,i) = ldsa_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    10569                       ENDDO
    10570                    ENDDO
    10571                 ENDDO
    10572 
    10573              CASE ( 'N_UFP' )
    10574                 DO  i = nxlg, nxrg
    10575                    DO  j = nysg, nyng
    10576                       DO  k = nzb, nzt+1
    10577                          nufp_av(k,j,i) = nufp_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    10578                       ENDDO
    10579                    ENDDO
    10580                 ENDDO
    10581 
    10582              CASE ( 'Ntot' )
    10583                 DO  i = nxlg, nxrg
    10584                    DO  j = nysg, nyng
    10585                       DO  k = nzb, nzt+1
    10586                          ntot_av(k,j,i) = ntot_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    10587                       ENDDO
    10588                    ENDDO
    10589                 ENDDO
    10590 
    10591 
    10592              CASE ( 'PM0.1' )
    10593                 DO  i = nxlg, nxrg
    10594                    DO  j = nysg, nyng
    10595                       DO  k = nzb, nzt+1
    10596                          pm01_av(k,j,i) = pm01_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    10597                       ENDDO
    10598                    ENDDO
    10599                 ENDDO
    10600 
    10601              CASE ( 'PM2.5' )
    10602                 DO  i = nxlg, nxrg
    10603                    DO  j = nysg, nyng
    10604                       DO  k = nzb, nzt+1
    10605                          pm25_av(k,j,i) = pm25_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    10606                       ENDDO
    10607                    ENDDO
    10608                 ENDDO
    10609 
    10610              CASE ( 'PM10' )
    10611                 DO  i = nxlg, nxrg
    10612                    DO  j = nysg, nyng
    10613                       DO  k = nzb, nzt+1
    10614                          pm10_av(k,j,i) = pm10_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    10615                       ENDDO
    10616                    ENDDO
    10617                 ENDDO
    10618 
    10619              CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
    10620                 IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN  ! 9: remove salsa_s_
    10621                    IF ( TRIM( variable(9:) ) == 'BC' )   to_be_resorted => s_bc_av
    10622                    IF ( TRIM( variable(9:) ) == 'DU' )   to_be_resorted => s_du_av
    10623                    IF ( TRIM( variable(9:) ) == 'NH' )   to_be_resorted => s_nh_av
    10624                    IF ( TRIM( variable(9:) ) == 'NO' )   to_be_resorted => s_no_av
    10625                    IF ( TRIM( variable(9:) ) == 'OC' )   to_be_resorted => s_oc_av
    10626                    IF ( TRIM( variable(9:) ) == 'SO4' )  to_be_resorted => s_so4_av
    10627                    IF ( TRIM( variable(9:) ) == 'SS' )   to_be_resorted => s_ss_av
     10681                IF ( ALLOCATED( ldsa_av ) )  THEN
    1062810682                   DO  i = nxlg, nxrg
    1062910683                      DO  j = nysg, nyng
    1063010684                         DO  k = nzb, nzt+1
    10631                             to_be_resorted(k,j,i) = to_be_resorted(k,j,i) /                        &
    10632                                                     REAL( average_count_3d, KIND=wp )
     10685                            ldsa_av(k,j,i) = ldsa_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    1063310686                         ENDDO
    1063410687                      ENDDO
    1063510688                   ENDDO
     10689                ENDIF
     10690
     10691             CASE ( 'N_UFP' )
     10692                IF ( ALLOCATED( nufp_av ) )  THEN
     10693                   DO  i = nxlg, nxrg
     10694                      DO  j = nysg, nyng
     10695                         DO  k = nzb, nzt+1
     10696                            nufp_av(k,j,i) = nufp_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     10697                         ENDDO
     10698                      ENDDO
     10699                   ENDDO
     10700                ENDIF
     10701
     10702             CASE ( 'Ntot' )
     10703                IF ( ALLOCATED( ntot_av ) )  THEN
     10704                   DO  i = nxlg, nxrg
     10705                      DO  j = nysg, nyng
     10706                         DO  k = nzb, nzt+1
     10707                            ntot_av(k,j,i) = ntot_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     10708                         ENDDO
     10709                      ENDDO
     10710                   ENDDO
     10711                ENDIF
     10712
     10713             CASE ( 'PM0.1' )
     10714                IF ( ALLOCATED( pm01_av ) )  THEN
     10715                   DO  i = nxlg, nxrg
     10716                      DO  j = nysg, nyng
     10717                         DO  k = nzb, nzt+1
     10718                            pm01_av(k,j,i) = pm01_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     10719                         ENDDO
     10720                      ENDDO
     10721                   ENDDO
     10722                ENDIF
     10723
     10724             CASE ( 'PM2.5' )
     10725                IF ( ALLOCATED( pm25_av ) )  THEN
     10726                   DO  i = nxlg, nxrg
     10727                      DO  j = nysg, nyng
     10728                         DO  k = nzb, nzt+1
     10729                            pm25_av(k,j,i) = pm25_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     10730                         ENDDO
     10731                      ENDDO
     10732                   ENDDO
     10733                ENDIF
     10734
     10735             CASE ( 'PM10' )
     10736                IF ( ALLOCATED( pm10_av ) )  THEN
     10737                   DO  i = nxlg, nxrg
     10738                      DO  j = nysg, nyng
     10739                         DO  k = nzb, nzt+1
     10740                            pm10_av(k,j,i) = pm10_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     10741                         ENDDO
     10742                      ENDDO
     10743                   ENDDO
     10744                ENDIF
     10745
     10746             CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
     10747                IF ( ALLOCATED( s_mass_av ) )  THEN
     10748                   IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN  ! 9: remove salsa_s_
     10749                      found_index = get_index( prtcl, TRIM( variable(9:) ) )
     10750                      DO  i = nxlg, nxrg
     10751                         DO  j = nysg, nyng
     10752                            DO  k = nzb, nzt+1
     10753                               s_mass_av(k,j,i,found_index) = s_mass_av(k,j,i,found_index) /       &
     10754                                                              REAL( average_count_3d, KIND=wp )
     10755                            ENDDO
     10756                         ENDDO
     10757                      ENDDO
     10758                   ENDIF
    1063610759                ENDIF
    1063710760
     
    1069610819
    1069710820    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< output
    10698 
    10699     REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted           !< pointer
    1070010821!
    1070110822!-- Next statement is to avoid compiler warning about unused variable. May be removed in future.
     
    1072110842             ENDDO
    1072210843          ELSE
     10844             IF ( .NOT. ALLOCATED( nbins_av ) )  THEN
     10845                ALLOCATE( nbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
     10846                nbins_av = REAL( fill_value, KIND = wp )
     10847             ENDIF
    1072310848             DO  i = nxl, nxr
    1072410849                DO  j = nys, nyn
     
    1075310878             ENDDO
    1075410879          ELSE
     10880             IF ( .NOT. ALLOCATED( mbins_av ) )  THEN
     10881                ALLOCATE( mbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
     10882                mbins_av = REAL( fill_value, KIND = wp )
     10883             ENDIF
    1075510884             DO  i = nxl, nxr
    1075610885                DO  j = nys, nyn
     
    1077110900          CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
    1077210901             vari = TRIM( variable( 9:LEN( TRIM( variable ) ) - 3 ) )  ! 9: remove salsa_g_
     10902             IF ( vari == 'H2SO4')  found_index = 1
     10903             IF ( vari == 'HNO3')   found_index = 2
     10904             IF ( vari == 'NH3')    found_index = 3
     10905             IF ( vari == 'OCNV')   found_index = 4
     10906             IF ( vari == 'OCSV')   found_index = 5
    1077310907             IF ( av == 0 )  THEN
    10774                 IF ( vari == 'H2SO4')  found_index = 1
    10775                 IF ( vari == 'HNO3')   found_index = 2
    10776                 IF ( vari == 'NH3')    found_index = 3
    10777                 IF ( vari == 'OCNV')   found_index = 4
    10778                 IF ( vari == 'OCSV')   found_index = 5
    1077910908                DO  i = nxl, nxr
    1078010909                   DO  j = nys, nyn
     
    1078710916                ENDDO
    1078810917             ELSE
    10789                 IF ( vari == 'H2SO4' )  to_be_resorted => g_h2so4_av
    10790                 IF ( vari == 'HNO3' )   to_be_resorted => g_hno3_av
    10791                 IF ( vari == 'NH3' )    to_be_resorted => g_nh3_av
    10792                 IF ( vari == 'OCNV' )   to_be_resorted => g_ocnv_av
    10793                 IF ( vari == 'OCSV' )   to_be_resorted => g_ocsv_av
     10918                IF ( .NOT. ALLOCATED( salsa_gases_av ) )  THEN
     10919                   ALLOCATE( salsa_gases_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
     10920                   salsa_gases_av = REAL( fill_value, KIND = wp )
     10921                ENDIF
    1079410922                DO  i = nxl, nxr
    1079510923                   DO  j = nys, nyn
    1079610924                      DO  k = nzb_do, nzt_do
    10797                          local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
    10798                                                KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
     10925                         local_pf(i,j,k) = MERGE( salsa_gases_av(k,j,i,found_index),               &
     10926                                                  REAL( fill_value, KIND = wp ),                   &
     10927                                                  BTEST( wall_flags_total_0(k,j,i), 0 ) )
    1079910928                      ENDDO
    1080010929                   ENDDO
     
    1081110940                         temp_bin = 0.0_wp
    1081210941                         DO  ib = 1, nbins_aerosol
    10813    !
    10814    !--                      Diameter in micrometres
     10942!
     10943!--                         Diameter in micrometres
    1081510944                            mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
    10816    !
    10817    !--                      Deposition factor: alveolar
     10945!
     10946!--                         Deposition factor: alveolar
    1081810947                            df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
    1081910948                                   2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
    1082010949                                   1.362_wp )**2 ) )
    10821    !
    10822    !--                      Lung-deposited surface area LDSA (units mum2/cm3)
     10950!
     10951!--                         Lung-deposited surface area LDSA (units mum2/cm3)
    1082310952                            temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
    1082410953                                       aerosol_number(ib)%conc(k,j,i)
     
    1083110960                ENDDO
    1083210961             ELSE
     10962                IF ( .NOT. ALLOCATED( ldsa_av ) )  THEN
     10963                   ALLOCATE( ldsa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     10964                   ldsa_av = REAL( fill_value, KIND = wp )
     10965                ENDIF
    1083310966                DO  i = nxl, nxr
    1083410967                   DO  j = nys, nyn
     
    1086110994                ENDDO
    1086210995             ELSE
     10996                IF ( .NOT. ALLOCATED( nufp_av ) )  THEN
     10997                   ALLOCATE( nufp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     10998                   nufp_av = REAL( fill_value, KIND = wp )
     10999                ENDIF
    1086311000                DO  i = nxl, nxr
    1086411001                   DO  j = nys, nyn
     
    1088911026                ENDDO
    1089011027             ELSE
     11028                IF ( .NOT. ALLOCATED( ntot_av ) )  THEN
     11029                   ALLOCATE( ntot_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11030                   ntot_av = REAL( fill_value, KIND = wp )
     11031                ENDIF
    1089111032                DO  i = nxl, nxr
    1089211033                   DO  j = nys, nyn
     
    1092011061                ENDDO
    1092111062             ELSE
     11063                IF ( .NOT. ALLOCATED( pm01_av ) )  THEN
     11064                   ALLOCATE( pm01_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11065                   pm01_av = REAL( fill_value, KIND = wp )
     11066                ENDIF
    1092211067                DO  i = nxl, nxr
    1092311068                   DO  j = nys, nyn
     
    1095111096                ENDDO
    1095211097             ELSE
     11098                IF ( .NOT. ALLOCATED( pm25_av ) )  THEN
     11099                   ALLOCATE( pm25_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11100                   pm25_av = REAL( fill_value, KIND = wp )
     11101                ENDIF
    1095311102                DO  i = nxl, nxr
    1095411103                   DO  j = nys, nyn
     
    1098211131                ENDDO
    1098311132             ELSE
     11133                IF ( .NOT. ALLOCATED( pm10_av ) )  THEN
     11134                   ALLOCATE( pm10_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11135                   pm10_av = REAL( fill_value, KIND = wp )
     11136                ENDIF
    1098411137                DO  i = nxl, nxr
    1098511138                   DO  j = nys, nyn
     
    1101211165                   ENDDO
    1101311166                ELSE
    11014                    IF ( vari == 'BC' )   to_be_resorted => s_bc_av
    11015                    IF ( vari == 'DU' )   to_be_resorted => s_du_av
    11016                    IF ( vari == 'NH' )   to_be_resorted => s_nh_av
    11017                    IF ( vari == 'NO' )   to_be_resorted => s_no_av
    11018                    IF ( vari == 'OC' )   to_be_resorted => s_oc_av
    11019                    IF ( vari == 'SO4' )  to_be_resorted => s_so4_av
    11020                    IF ( vari == 'SS' )   to_be_resorted => s_ss_av
     11167                   IF ( .NOT. ALLOCATED( s_mass_av ) )  THEN
     11168                      ALLOCATE( s_mass_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass) )
     11169                      s_mass_av = REAL( fill_value, KIND = wp )
     11170                   ENDIF
    1102111171                   DO  i = nxl, nxr
    1102211172                      DO  j = nys, nyn
    1102311173                         DO  k = nzb_do, nzt_do
    11024                             local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,      &
    11025                                                KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
     11174                            local_pf(i,j,k) = MERGE( s_mass_av(k,j,i,found_index),                 &
     11175                                                     REAL( fill_value, KIND = wp ),                &
     11176                                                     BTEST( wall_flags_total_0(k,j,i), 0 ) )
    1102611177                         ENDDO
    1102711178                      ENDDO
     
    1105011201                ENDDO
    1105111202             ELSE
    11052                 to_be_resorted => s_h2o_av
     11203     !           to_be_resorted => s_h2o_av
     11204                IF ( .NOT. ALLOCATED( s_h2o_av ) )  THEN
     11205                   ALLOCATE( s_h2o_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11206                   s_h2o_av = REAL( fill_value, KIND = wp )
     11207                ENDIF
    1105311208                DO  i = nxl, nxr
    1105411209                   DO  j = nys, nyn
    1105511210                      DO  k = nzb_do, nzt_do
    11056                          local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
    11057                                               KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
     11211                         local_pf(i,j,k) = MERGE( s_h2o_av(k,j,i), REAL( fill_value, KIND = wp ),  &
     11212                                                  BTEST( wall_flags_total_0(k,j,i), 0 ) )
    1105811213                      ENDDO
    1105911214                   ENDDO
     
    1111011265
    1111111266    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< local
    11112 
    11113     REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< pointer
    1111411267
    1111511268    found     = .TRUE.
     
    1113111284             ENDDO
    1113211285          ELSE
     11286             IF ( .NOT. ALLOCATED( nbins_av ) )  THEN
     11287                ALLOCATE( nbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
     11288                nbins_av = REAL( fill_value, KIND = wp )
     11289             ENDIF
    1113311290             DO  i = nxl, nxr
    1113411291                DO  j = nys, nyn
     
    1116111318             ENDDO
    1116211319          ELSE
     11320             IF ( .NOT. ALLOCATED( mbins_av ) )  THEN
     11321                ALLOCATE( mbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
     11322                mbins_av = REAL( fill_value, KIND = wp )
     11323             ENDIF
    1116311324             DO  i = nxl, nxr
    1116411325                DO  j = nys, nyn
     
    1117611337
    1117711338          CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV' )
     11339             IF ( TRIM( variable(7:) ) == 'g_H2SO4')  found_index = 1
     11340             IF ( TRIM( variable(7:) ) == 'g_HNO3')   found_index = 2
     11341             IF ( TRIM( variable(7:) ) == 'g_NH3')    found_index = 3
     11342             IF ( TRIM( variable(7:) ) == 'g_OCNV')   found_index = 4
     11343             IF ( TRIM( variable(7:) ) == 'g_OCSV')   found_index = 5
     11344
    1117811345             IF ( av == 0 )  THEN
    11179                 IF ( TRIM( variable(7:) ) == 'g_H2SO4')  found_index = 1
    11180                 IF ( TRIM( variable(7:) ) == 'g_HNO3')   found_index = 2
    11181                 IF ( TRIM( variable(7:) ) == 'g_NH3')    found_index = 3
    11182                 IF ( TRIM( variable(7:) ) == 'g_OCNV')   found_index = 4
    11183                 IF ( TRIM( variable(7:) ) == 'g_OCSV')   found_index = 5
    11184 
    1118511346                DO  i = nxl, nxr
    1118611347                   DO  j = nys, nyn
     
    1119311354                ENDDO
    1119411355             ELSE
    11195 !
    11196 !--             9: remove salsa_g_ from the beginning
    11197                 IF ( TRIM( variable(9:) ) == 'H2SO4' ) to_be_resorted => g_h2so4_av
    11198                 IF ( TRIM( variable(9:) ) == 'HNO3' )  to_be_resorted => g_hno3_av
    11199                 IF ( TRIM( variable(9:) ) == 'NH3' )   to_be_resorted => g_nh3_av
    11200                 IF ( TRIM( variable(9:) ) == 'OCNV' )  to_be_resorted => g_ocnv_av
    11201                 IF ( TRIM( variable(9:) ) == 'OCSV' )  to_be_resorted => g_ocsv_av
     11356                IF ( .NOT. ALLOCATED( salsa_gases_av ) )  THEN
     11357                   ALLOCATE( salsa_gases_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
     11358                   salsa_gases_av = REAL( fill_value, KIND = wp )
     11359                ENDIF
    1120211360                DO  i = nxl, nxr
    1120311361                   DO  j = nys, nyn
    1120411362                      DO  k = nzb_do, nzt_do
    11205                          local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
    11206                                                KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
     11363!                          local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
     11364!                                                KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
     11365                         local_pf(i,j,k) = MERGE( salsa_gases_av(k,j,i,found_index),               &
     11366                                                  REAL( fill_value, KIND = wp ),                   &
     11367                                                  BTEST( wall_flags_total_0(k,j,i), 0 ) )
    1120711368                      ENDDO
    1120811369                   ENDDO
     
    1121711378                         temp_bin = 0.0_wp
    1121811379                         DO  ib = 1, nbins_aerosol
    11219    !
    11220    !--                      Diameter in micrometres
     11380!
     11381!--                         Diameter in micrometres
    1122111382                            mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
    11222    !
    11223    !--                      Deposition factor: alveolar
     11383!
     11384!--                         Deposition factor: alveolar
    1122411385                            df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
    1122511386                                   2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
    1122611387                                   1.362_wp )**2 ) )
    11227    !
    11228    !--                      Lung-deposited surface area LDSA (units mum2/cm3)
     11388!
     11389!--                         Lung-deposited surface area LDSA (units mum2/cm3)
    1122911390                            temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
    1123011391                                       aerosol_number(ib)%conc(k,j,i)
     
    1123611397                ENDDO
    1123711398             ELSE
     11399                IF ( .NOT. ALLOCATED( ldsa_av ) )  THEN
     11400                   ALLOCATE( ldsa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11401                   ldsa_av = REAL( fill_value, KIND = wp )
     11402                ENDIF
    1123811403                DO  i = nxl, nxr
    1123911404                   DO  j = nys, nyn
     
    1126311428                ENDDO
    1126411429             ELSE
     11430                IF ( .NOT. ALLOCATED( nufp_av ) )  THEN
     11431                   ALLOCATE( nufp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11432                   nufp_av = REAL( fill_value, KIND = wp )
     11433                ENDIF
    1126511434                DO  i = nxl, nxr
    1126611435                   DO  j = nys, nyn
     
    1128811457                ENDDO
    1128911458             ELSE
     11459                IF ( .NOT. ALLOCATED( ntot_av ) )  THEN
     11460                   ALLOCATE( ntot_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11461                   ntot_av = REAL( fill_value, KIND = wp )
     11462                ENDIF
    1129011463                DO  i = nxl, nxr
    1129111464                   DO  j = nys, nyn
     
    1131711490                ENDDO
    1131811491             ELSE
     11492                IF ( .NOT. ALLOCATED( pm01_av ) )  THEN
     11493                   ALLOCATE( pm01_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11494                   pm01_av = REAL( fill_value, KIND = wp )
     11495                ENDIF
    1131911496                DO  i = nxl, nxr
    1132011497                   DO  j = nys, nyn
     
    1134611523                ENDDO
    1134711524             ELSE
     11525                IF ( .NOT. ALLOCATED( pm25_av ) )  THEN
     11526                   ALLOCATE( pm25_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11527                   pm25_av = REAL( fill_value, KIND = wp )
     11528                ENDIF
    1134811529                DO  i = nxl, nxr
    1134911530                   DO  j = nys, nyn
     
    1137511556                ENDDO
    1137611557             ELSE
     11558                IF ( .NOT. ALLOCATED( pm10_av ) )  THEN
     11559                   ALLOCATE( pm10_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11560                   pm10_av = REAL( fill_value, KIND = wp )
     11561                ENDIF
    1137711562                DO  i = nxl, nxr
    1137811563                   DO  j = nys, nyn
     
    1140211587                   ENDDO
    1140311588                ELSE
    11404 !
    11405 !--                9: remove salsa_s_ from the beginning
    11406                    IF ( TRIM( variable(9:) ) == 'BC' )   to_be_resorted => s_bc_av
    11407                    IF ( TRIM( variable(9:) ) == 'DU' )   to_be_resorted => s_du_av
    11408                    IF ( TRIM( variable(9:) ) == 'NH' )   to_be_resorted => s_nh_av
    11409                    IF ( TRIM( variable(9:) ) == 'NO' )   to_be_resorted => s_no_av
    11410                    IF ( TRIM( variable(9:) ) == 'OC' )   to_be_resorted => s_oc_av
    11411                    IF ( TRIM( variable(9:) ) == 'SO4' )  to_be_resorted => s_so4_av
    11412                    IF ( TRIM( variable(9:) ) == 'SS' )   to_be_resorted => s_ss_av
     11589                   IF ( .NOT. ALLOCATED( s_mass_av ) )  THEN
     11590                      ALLOCATE( s_mass_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass) )
     11591                      s_mass_av = REAL( fill_value, KIND = wp )
     11592                   ENDIF
    1141311593                   DO  i = nxl, nxr
    1141411594                      DO  j = nys, nyn
    1141511595                         DO  k = nzb_do, nzt_do
    11416                             local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,      &
    11417                                                      KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
     11596                            local_pf(i,j,k) = MERGE( s_mass_av(k,j,i,found_index),                 &
     11597                                                     REAL( fill_value, KIND = wp ),                &
     11598                                                     BTEST( wall_flags_total_0(k,j,i), 0 ) )
    1141811599                         ENDDO
    1141911600                      ENDDO
     
    1143811619                ENDDO
    1143911620             ELSE
    11440                 to_be_resorted => s_h2o_av
     11621                IF ( .NOT. ALLOCATED( s_h2o_av ) )  THEN
     11622                   ALLOCATE( s_h2o_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     11623                   s_h2o_av = REAL( fill_value, KIND = wp )
     11624                ENDIF
    1144111625                DO  i = nxl, nxr
    1144211626                   DO  j = nys, nyn
    1144311627                      DO  k = nzb_do, nzt_do
    11444                          local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
    11445                                                KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
     11628                         local_pf(i,j,k) = MERGE( s_h2o_av(k,j,i), REAL( fill_value, KIND = wp ),  &
     11629                                                 BTEST( wall_flags_total_0(k,j,i), 0 ) )
    1144611630                      ENDDO
    1144711631                   ENDDO
     
    1162211806                IF ( vari == 'g_OCSV')   to_be_resorted => salsa_gas(5)%conc
    1162311807             ELSE
    11624                 IF ( vari == 'g_H2SO4') to_be_resorted => g_h2so4_av
    11625                 IF ( vari == 'g_HNO3')  to_be_resorted => g_hno3_av
    11626                 IF ( vari == 'g_NH3')   to_be_resorted => g_nh3_av
    11627                 IF ( vari == 'g_OCNV')  to_be_resorted => g_ocnv_av
    11628                 IF ( vari == 'g_OCSV')  to_be_resorted => g_ocsv_av
     11808                IF ( vari == 'g_H2SO4') temp_array = salsa_gases_av(:,:,:,1)
     11809                IF ( vari == 'g_HNO3')  temp_array = salsa_gases_av(:,:,:,2)
     11810                IF ( vari == 'g_NH3')   temp_array = salsa_gases_av(:,:,:,3)
     11811                IF ( vari == 'g_OCNV')  temp_array = salsa_gases_av(:,:,:,4)
     11812                IF ( vari == 'g_OCSV')  temp_array = salsa_gases_av(:,:,:,5)
     11813                to_be_resorted => temp_array
    1162911814             ENDIF
    1163011815
     
    1163611821                         temp_bin = 0.0_wp
    1163711822                         DO  ib = 1, nbins_aerosol
    11638    !
    11639    !--                      Diameter in micrometres
     11823!
     11824!--                         Diameter in micrometres
    1164011825                            mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
    11641    !
    11642    !--                      Deposition factor: alveolar
     11826!
     11827!--                         Deposition factor: alveolar
    1164311828                            df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
    1164411829                                   2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
    1164511830                                   1.362_wp )**2 ) )
    11646    !
    11647    !--                      Lung-deposited surface area LDSA (units mum2/cm3)
     11831!
     11832!--                         Lung-deposited surface area LDSA (units mum2/cm3)
    1164811833                            temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
    1164911834                                       aerosol_number(ib)%conc(k,j,i)
     
    1185412039                      ENDDO
    1185512040                   ENDDO
    11856                 ENDDO 
     12041                ENDDO
    1185712042                IF ( .NOT. mask_surface(mid) )  THEN
    1185812043                   DO  i = 1, mask_size_l(mid,1)
     
    1194312128
    1194412129          CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
    11945              IF ( av == 0 )  THEN
    11946                 IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN
    11947                    found_index = get_index( prtcl, TRIM( variable(9:) ) )
     12130             IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN
     12131                found_index = get_index( prtcl, TRIM( variable(9:) ) )
     12132                IF ( av == 0 )  THEN
    1194812133                   DO  i = nxl, nxr
    1194912134                      DO  j = nys, nyn
     
    1195712142                      ENDDO
    1195812143                   ENDDO
    11959                 ELSE
    11960                    tend = 0.0_wp
    11961                 ENDIF
    11962                 IF ( .NOT. mask_surface(mid) )  THEN
    11963                    DO  i = 1, mask_size_l(mid,1)
    11964                       DO  j = 1, mask_size_l(mid,2)
    11965                          DO  k = 1, mask_size_l(mid,3)
    11966                             local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) )
     12144                   IF ( .NOT. mask_surface(mid) )  THEN
     12145                      DO  i = 1, mask_size_l(mid,1)
     12146                         DO  j = 1, mask_size_l(mid,2)
     12147                            DO  k = 1, mask_size_l(mid,3)
     12148                               local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) )
     12149                            ENDDO
    1196712150                         ENDDO
    1196812151                      ENDDO
    11969                    ENDDO
    11970                 ELSE
    11971                    DO  i = 1, mask_size_l(mid,1)
    11972                       DO  j = 1, mask_size_l(mid,2)
    11973 !
    11974 !--                      Get k index of the highest terraing surface
    11975                          im = mask_i(mid,i)
    11976                          jm = mask_j(mid,j)
    11977                          ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
    11978                                                        DIM = 1 ) - 1
    11979                          DO  k = 1, mask_size_l(mid,3)
    11980                             kk = MIN( ktt+mask_k(mid,k), nzt+1 )
    11981 !
    11982 !--                         Set value if not in building
    11983                             IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
    11984                                local_pf(i,j,k) = fill_value
    11985                             ELSE
    11986                                local_pf(i,j,k) = tend(kk,jm,im)
    11987                             ENDIF
     12152                   ELSE
     12153                      DO  i = 1, mask_size_l(mid,1)
     12154                         DO  j = 1, mask_size_l(mid,2)
     12155   !
     12156   !--                      Get k index of the highest terraing surface
     12157                            im = mask_i(mid,i)
     12158                            jm = mask_j(mid,j)
     12159                            ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
     12160                                                          DIM = 1 ) - 1
     12161                            DO  k = 1, mask_size_l(mid,3)
     12162                               kk = MIN( ktt+mask_k(mid,k), nzt+1 )
     12163   !
     12164   !--                         Set value if not in building
     12165                               IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
     12166                                  local_pf(i,j,k) = fill_value
     12167                               ELSE
     12168                                  local_pf(i,j,k) = tend(kk,jm,im)
     12169                               ENDIF
     12170                            ENDDO
    1198812171                         ENDDO
    1198912172                      ENDDO
    11990                    ENDDO
     12173                   ENDIF
     12174                   resorted = .TRUE.
     12175                ELSE
     12176                   temp_array = s_mass_av(:,:,:,found_index)
     12177                   to_be_resorted => temp_array
    1199112178                ENDIF
    11992                 resorted = .TRUE.
    1199312179             ELSE
    11994 !
    11995 !--             9: remove salsa_s_ from the beginning
    11996                 IF ( TRIM( variable(9:) ) == 'BC' )   to_be_resorted => s_bc_av
    11997                 IF ( TRIM( variable(9:) ) == 'DU' )   to_be_resorted => s_du_av
    11998                 IF ( TRIM( variable(9:) ) == 'NH' )   to_be_resorted => s_nh_av
    11999                 IF ( TRIM( variable(9:) ) == 'NO' )   to_be_resorted => s_no_av
    12000                 IF ( TRIM( variable(9:) ) == 'OC' )   to_be_resorted => s_oc_av
    12001                 IF ( TRIM( variable(9:) ) == 'SO4' )  to_be_resorted => s_so4_av
    12002                 IF ( TRIM( variable(9:) ) == 'SS' )   to_be_resorted => s_ss_av
     12180                local_pf = fill_value
    1200312181             ENDIF
    1200412182
     
    1207412252          DO  i = 1, mask_size_l(mid,1)
    1207512253             DO  j = 1, mask_size_l(mid,2)
     12254!
    1207612255!--             Get k index of the highest terraing surface
    1207712256                im = mask_i(mid,i)
     
    1208112260                DO  k = 1, mask_size_l(mid,3)
    1208212261                   kk = MIN( ktt+mask_k(mid,k), nzt+1 )
     12262!
    1208312263!--                Set value if not in building
    1208412264                   IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
Note: See TracChangeset for help on using the changeset viewer.