Changeset 3848 for palm/trunk/SOURCE


Ignore:
Timestamp:
Apr 1, 2019 3:25:48 PM (3 years ago)
Author:
forkel
Message:

some formatting

File:
1 edited

Legend:

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

    r3833 r3848  
    2727! -----------------
    2828! $Id: chemistry_model_mod.f90 3784 2019-03-05 14:16:20Z banzhafs
     29! some formatting
     30!
     31! 3784 2019-03-05 14:16:20Z banzhafs
    2932! added cs_mech to USE chem_gasphase_mod
    3033!
     
    260263! ------------
    261264!> Chemistry model for PALM-4U
    262 !> @todo Extend chem_species type by nspec and nvar as addititional elements
     265!> @todo Extend chem_species type by nspec and nvar as addititional elements (RF)
     266!> @todo Check possibility to reduce dimension of chem_species%conc from nspec to nvar (RF)
    263267!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not
    264268!>       allowed to use the chemistry model in a precursor run and additionally
    265269!>       not using it in a main run
    266 !> @todo Update/clean-up todo list! (FK)
    267 !> @todo Add routine chem_check_parameters, add checks for inconsistent or
    268 !>       unallowed parameter settings.
    269 !>       CALL of chem_check_parameters from check_parameters. (FK)
    270270!> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
    271271!> @todo Separate boundary conditions for each chem spcs to be implemented
     
    277277!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
    278278!> @todo chemistry error messages
    279 !> @todo Format this module according to PALM coding standard (see e.g. module
    280 !>       template under http://palm.muk.uni-hannover.de/mosaik/downloads/8 or
    281 !>       D3_coding_standard.pdf under https://palm.muk.uni-hannover.de/trac/downloads/16)
    282279!
    283280!------------------------------------------------------------------------------!
     
    318315
    319316    LOGICAL ::  nest_chemistry = .TRUE.  !< flag for nesting mode of chemical species, independent on parent or not
    320     !
    321     !- Define chemical variables
    322     TYPE   species_def
    323        CHARACTER(LEN=15)                            ::  name
    324        CHARACTER(LEN=15)                            ::  unit
    325        REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc
    326        REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc_av
    327        REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc_p
    328        REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  tconc_m
    329        REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  cssws_av          
    330        REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  flux_s_cs
    331        REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  diss_s_cs
    332        REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) ::  flux_l_cs
    333        REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) ::  diss_l_cs
    334        REAL(kind=wp), ALLOCATABLE, DIMENSION(:)     ::  conc_pr_init
     317!
     318!-  Define chemical variables within chem_species
     319    TYPE species_def
     320       CHARACTER(LEN=15)                            ::  name         !< name of chemical species
     321       CHARACTER(LEN=15)                            ::  unit         !< unit (ppm for gases, kg m^-3 for aerosol tracers)
     322       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc         !< concentrations of trace gases
     323       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc_av      !< averaged concentrations
     324       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  conc_p       !< conc at prognostic time level
     325       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     ::  tconc_m      !< weighted tendency of conc for previous sub-timestep (Runge-Kutta)
     326       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  cssws_av     !< averaged fluxes of trace gases at surface
     327       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  flux_s_cs    !< 6th-order advective flux at south face of grid box of chemical species (='cs')
     328       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   ::  diss_s_cs    !< artificial numerical dissipation flux at south face of grid box of chemical species
     329       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) ::  flux_l_cs    !< 6th-order advective flux at left face of grid box of chemical species (='cs')
     330       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) ::  diss_l_cs    !< artificial numerical dissipation flux at left face of grid box of chemical species
     331       REAL(kind=wp), ALLOCATABLE, DIMENSION(:)     ::  conc_pr_init !< initial profile of chemical species
    335332    END TYPE species_def
    336 
    337 
    338     TYPE   photols_def                                                           
    339        CHARACTER(LEN=15)                            :: name
    340        CHARACTER(LEN=15)                            :: unit
    341        REAL(kind=wp), POINTER, DIMENSION(:,:,:)     :: freq
     333!
     334!-- Define photolysis frequencies in phot_frequen
     335    TYPE photols_def                                                           
     336       CHARACTER(LEN=15)                            :: name          !< name of pgotolysis frequency
     337       CHARACTER(LEN=15)                            :: unit          !< unit (1/s)
     338       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     :: freq          !< photolysis frequency
    342339    END TYPE photols_def
    343 
    344 
    345     PUBLIC  species_def
    346     PUBLIC  photols_def
    347340
    348341
     
    350343    TYPE(photols_def), ALLOCATABLE, DIMENSION(:), TARGET, PUBLIC ::  phot_frequen 
    351344
    352     REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_1
    353     REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_2
    354     REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_3
    355     REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_av       
    356     REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  freq_1
    357 
    358     INTEGER, DIMENSION(nkppctrl)                           ::  icntrl                            !< Fine tuning kpp
    359     REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl                            !< Fine tuning kpp
    360     LOGICAL ::  decycle_chem_lr    = .FALSE.
    361     LOGICAL ::  decycle_chem_ns    = .FALSE.
     345    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_1  !< pointer for swapping of timelevels for conc
     346    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_2  !< pointer for swapping of timelevels for conc
     347    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_3  !< pointer for swapping of timelevels for conc
     348    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_av !< averaged concentrations of chemical species       
     349    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  freq_1       !< pointer for phtolysis frequncies
     350                                                                            !< (only 1 timelevel required)
     351    INTEGER, DIMENSION(nkppctrl)                           ::  icntrl       !< 20 integer parameters for fine tuning KPP code
     352                                                                            !< (e.g. solver type)
     353    REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl       !< 20 real parameters for fine tuning of KPP code
     354                                                                            !< (e.g starting internal timestep of solver)
     355!
     356!-- Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
     357!-- approproate for chemicals compounds since they may accumulate too much.
     358!-- If no proper boundary conditions from a DYNAMIC input file are available,
     359!-- de-cycling applies the initial profiles at the inflow boundaries for
     360!-- Dirichlet boundary conditions
     361    LOGICAL ::  decycle_chem_lr    = .FALSE.    !< switch for setting decycling in left-right direction
     362    LOGICAL ::  decycle_chem_ns    = .FALSE.    !< switch for setting decycling in south-norht direction
    362363    CHARACTER (LEN=20), DIMENSION(4) ::  decycle_method = &
    363364         (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
     
    376377
    377378
    378     !
    379     !-- Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
     379!
     380!-- Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
    380381    !
    381382    INTEGER(iwp), PARAMETER ::  nlu_dep = 15                   !< Number of DEPAC landuse classes (lu's)
    382383    INTEGER(iwp), PARAMETER ::  ncmp = 10                      !< Number of DEPAC gas components
    383384    INTEGER(iwp), PARAMETER ::  nposp = 69                     !< Number of possible species for deposition
    384     !
    385     !-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1                             
     385!
     386!-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1                             
    386387    INTEGER(iwp) ::  ilu_grass              = 1       
    387388    INTEGER(iwp) ::  ilu_arable             = 2       
     
    400401    INTEGER(iwp) ::  ilu_semi_natural_veg   = 15
    401402
    402     !
    403     !-- NH3/SO2 ratio regimes:
     403!
     404!-- NH3/SO2 ratio regimes:
    404405    INTEGER(iwp), PARAMETER ::  iratns_low      = 1       !< low ratio NH3/SO2
    405406    INTEGER(iwp), PARAMETER ::  iratns_high     = 2       !< high ratio NH3/SO2
    406407    INTEGER(iwp), PARAMETER ::  iratns_very_low = 3       !< very low ratio NH3/SO2
    407     !
    408     !-- Default:
     408!
     409!-- Default:
    409410    INTEGER, PARAMETER ::  iratns_default = iratns_low
    410     !
    411     !-- Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) to W m-2
     411!
     412!-- Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) to W m-2
    412413    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  alpha   =(/ 0.009, 0.009, 0.009, 0.006, 0.006, -999., -999., 0.009, -999.,      &
    413414         -999., 0.009, 0.006, -999., 0.009, 0.008/)*4.57
    414     !
    415     !-- Set temperatures per land use for f_temp
     415!
     416!-- Set temperatures per land use for f_temp
    416417    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmin = (/ 12.0, 12.0,  12.0,  0.0,  0.0, -999., -999., 12.0, -999., -999.,      &
    417418         12.0,  0.0, -999., 12.0,  8.0/)
     
    420421    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmax = (/ 40.0, 40.0,  40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999.,      &
    421422         40.0, 35.0, -999., 40.0, 39.0 /)
    422     !
    423     !-- Set f_min:
     423!
     424!-- Set f_min:
    424425    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  f_min = (/ 0.01, 0.01, 0.01, 0.1, 0.1, -999., -999., 0.01, -999., -999., 0.01,  &
    425426         0.1, -999., 0.01, 0.04/)
    426427
    427     !
    428     !-- Set maximal conductance (m/s)
    429     !-- (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from  mmol O3/m2/s to m/s
     428!
     429!-- Set maximal conductance (m/s)
     430!-- (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from  mmol O3/m2/s to m/s
    430431    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  g_max = (/ 270., 300., 300., 140., 150., -999., -999., 270., -999., -999.,      &
    431          270., 150., -999., 300., 422./)/41000
    432     !
    433     !-- Set max, min for vapour pressure deficit vpd
     432         270., 150., -999., 300., 422./)/41000.
     433!
     434!-- Set max, min for vapour pressure deficit vpd
    434435    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_max = (/1.3, 0.9, 0.9, 0.5, 1.0, -999., -999., 1.3, -999., -999., 1.3,      &
    435436         1.0, -999., 0.9, 2.8/)
     
    437438         3.25, -999., 2.8, 4.5/)
    438439
    439 
    440440    PUBLIC nest_chemistry
    441     PUBLIC nspec
    442441    PUBLIC nreact
    443     PUBLIC nvar     
    444     PUBLIC spc_names
     442    PUBLIC nspec               !< number of gas phase chemical species including constant compound (e.g. N2)
     443    PUBLIC nvar                !< number of variable gas phase chemical species (nvar <= nspec)
     444    PUBLIC photols_def         !< type defining phot_frequen
     445    PUBLIC species_def         !< type defining chem_species
     446    PUBLIC spc_names           !< names of gas phase chemical species (come from KPP) (come from KPP)
    445447    PUBLIC spec_conc_2 
    446 
    447     !   
    448     !-- Interface section
     448!   
     449!-- Interface section
    449450    INTERFACE chem_3d_data_averaging
    450451       MODULE PROCEDURE chem_3d_data_averaging
     
    607608         chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo
    608609
    609 
    610 
    611610 CONTAINS
    612611
    613  !
    614  !------------------------------------------------------------------------------!
    615  !
    616  ! Description:
    617  ! ------------
    618  !> Subroutine for averaging 3D data of chemical species. Due to the fact that
    619  !> the averaged chem arrays are allocated in chem_init, no if-query concerning
    620  !> the allocation is required (in any mode). Attention: If you just specify an
    621  !> averaged output quantity in the _p3dr file during restarts the first output
    622  !> includes the time between the beginning of the restart run and the first
    623  !> output time (not necessarily the whole averaging_interval you have
    624  !> specified in your _p3d/_p3dr file )
    625  !------------------------------------------------------------------------------!
    626 
     612
     613!------------------------------------------------------------------------------!
     614! Description:
     615! ------------
     616!> Subroutine for averaging 3D data of chemical species. Due to the fact that
     617!> the averaged chem arrays are allocated in chem_init, no if-query concerning
     618!> the allocation is required (in any mode). Attention: If you just specify an
     619!> averaged output quantity in the _p3dr file during restarts the first output
     620!> includes the time between the beginning of the restart run and the first
     621!> output time (not necessarily the whole averaging_interval you have
     622!> specified in your _p3d/_p3dr file )
     623!------------------------------------------------------------------------------!
    627624 SUBROUTINE chem_3d_data_averaging( mode, variable )
    628625
     
    638635    IMPLICIT NONE
    639636
    640     CHARACTER (LEN=*) ::  mode    !<
    641     CHARACTER (LEN=*) :: variable !<
    642 
    643     LOGICAL      ::  match_def !< flag indicating natural-type surface
    644     LOGICAL      ::  match_lsm !< flag indicating natural-type surface
    645     LOGICAL      ::  match_usm !< flag indicating urban-type surface
     637    CHARACTER (LEN=*) ::  mode     !<
     638    CHARACTER (LEN=*) ::  variable !<
     639
     640    LOGICAL ::  match_def  !< flag indicating default-type surface
     641    LOGICAL ::  match_lsm !< flag indicating natural-type surface
     642    LOGICAL ::  match_usm !< flag indicating urban-type surface
    646643
    647644    INTEGER(iwp) ::  i                  !< grid index x direction
     
    649646    INTEGER(iwp) ::  k                  !< grid index z direction
    650647    INTEGER(iwp) ::  m                  !< running index surface type
    651     INTEGER(iwp) :: lsp                 !< running index for chem spcs
     648    INTEGER(iwp) ::  lsp               !< running index for chem spcs
    652649
    653650    IF ( (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_')  )  THEN
    654651
    655     IF ( mode == 'allocate' )  THEN
    656        DO lsp = 1, nspec
    657           IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
    658                TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
    659              chem_species(lsp)%conc_av = 0.0_wp
    660           ENDIF
    661        ENDDO
    662 
    663     ELSEIF ( mode == 'sum' )  THEN
    664 
    665        DO lsp = 1, nspec
    666           IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
    667                TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
    668              DO  i = nxlg, nxrg
    669                 DO  j = nysg, nyng
    670                    DO  k = nzb, nzt+1
    671                       chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) + &
    672                            chem_species(lsp)%conc(k,j,i)
     652       IF ( mode == 'allocate' )  THEN
     653
     654          DO  lsp = 1, nspec
     655             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
     656                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
     657                chem_species(lsp)%conc_av = 0.0_wp
     658             ENDIF
     659          ENDDO
     660
     661       ELSEIF ( mode == 'sum' )  THEN
     662
     663          DO  lsp = 1, nspec
     664             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
     665                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
     666                DO  i = nxlg, nxrg
     667                   DO  j = nysg, nyng
     668                      DO  k = nzb, nzt+1
     669                         chem_species(lsp)%conc_av(k,j,i) =                              &
     670                                           chem_species(lsp)%conc_av(k,j,i) +            &
     671                                           chem_species(lsp)%conc(k,j,i)
     672                      ENDDO
    673673                   ENDDO
    674674                ENDDO
    675              ENDDO
    676           ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') )        THEN
    677              DO  i = nxl, nxr
    678                 DO  j = nys, nyn
    679                    match_def = surf_def_h(0)%start_index(j,i) <=               &
    680                         surf_def_h(0)%end_index(j,i)
    681                    match_lsm = surf_lsm_h%start_index(j,i) <=                  &
    682                         surf_lsm_h%end_index(j,i)
    683                    match_usm = surf_usm_h%start_index(j,i) <=                  &
    684                         surf_usm_h%end_index(j,i)
    685 
    686                    IF ( match_def )  THEN
    687                       m = surf_def_h(0)%end_index(j,i)
    688                       chem_species(lsp)%cssws_av(j,i) =                        &
    689                            chem_species(lsp)%cssws_av(j,i) + &
    690                            surf_def_h(0)%cssws(lsp,m)
    691                    ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
    692                       m = surf_lsm_h%end_index(j,i)
    693                       chem_species(lsp)%cssws_av(j,i) =                        &
    694                            chem_species(lsp)%cssws_av(j,i) + &
    695                            surf_lsm_h%cssws(lsp,m)
    696                    ELSEIF ( match_usm )  THEN
    697                       m = surf_usm_h%end_index(j,i)
    698                       chem_species(lsp)%cssws_av(j,i) =                        &
    699                            chem_species(lsp)%cssws_av(j,i) + &
    700                            surf_usm_h%cssws(lsp,m)
    701                    ENDIF
    702                 ENDDO
    703              ENDDO
    704           ENDIF
    705        ENDDO
    706 
    707     ELSEIF ( mode == 'average' )  THEN
    708        DO lsp = 1, nspec
    709           IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
    710                TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
    711              DO  i = nxlg, nxrg
    712                 DO  j = nysg, nyng
    713                    DO  k = nzb, nzt+1
    714                       chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     675             ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) )  THEN
     676                DO  i = nxl, nxr
     677                   DO  j = nys, nyn
     678                      match_def = surf_def_h(0)%start_index(j,i) <=                      &
     679                           surf_def_h(0)%end_index(j,i)
     680                      match_lsm = surf_lsm_h%start_index(j,i) <=                         &
     681                           surf_lsm_h%end_index(j,i)
     682                      match_usm = surf_usm_h%start_index(j,i) <=                         &
     683                           surf_usm_h%end_index(j,i)
     684
     685                      IF ( match_def )  THEN
     686                         m = surf_def_h(0)%end_index(j,i)
     687                         chem_species(lsp)%cssws_av(j,i) =                               &
     688                              chem_species(lsp)%cssws_av(j,i) + &
     689                              surf_def_h(0)%cssws(lsp,m)
     690                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
     691                         m = surf_lsm_h%end_index(j,i)
     692                         chem_species(lsp)%cssws_av(j,i) =                               &
     693                              chem_species(lsp)%cssws_av(j,i) + &
     694                              surf_lsm_h%cssws(lsp,m)
     695                      ELSEIF ( match_usm )  THEN
     696                         m = surf_usm_h%end_index(j,i)
     697                         chem_species(lsp)%cssws_av(j,i) =                               &
     698                              chem_species(lsp)%cssws_av(j,i) + &
     699                              surf_usm_h%cssws(lsp,m)
     700                      ENDIF
    715701                   ENDDO
    716702                ENDDO
    717              ENDDO
    718 
    719           ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') )            THEN
    720              DO i = nxlg, nxrg
    721                 DO  j = nysg, nyng
    722                    chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
     703             ENDIF
     704          ENDDO
     705
     706       ELSEIF ( mode == 'average' )  THEN
     707
     708          DO  lsp = 1, nspec
     709             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
     710                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
     711                DO  i = nxlg, nxrg
     712                   DO  j = nysg, nyng
     713                      DO  k = nzb, nzt+1
     714                         chem_species(lsp)%conc_av(k,j,i) =                              &
     715                             chem_species(lsp)%conc_av(k,j,i) /                          &
     716                             REAL( average_count_3d, KIND=wp )
     717                      ENDDO
     718                   ENDDO
    723719                ENDDO
    724              ENDDO
    725              CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )                 
    726           ENDIF
    727        ENDDO
     720
     721             ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) )  THEN
     722                DO  i = nxlg, nxrg
     723                   DO  j = nysg, nyng
     724                      chem_species(lsp)%cssws_av(j,i) =                                  &
     725                      chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
     726                   ENDDO
     727                ENDDO
     728                CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )                 
     729             ENDIF
     730          ENDDO
     731       ENDIF
    728732
    729733    ENDIF
    730734
    731     ENDIF
    732 
    733735 END SUBROUTINE chem_3d_data_averaging
    734736
    735 !   
    736 !------------------------------------------------------------------------------!
    737 !
     737   
     738!------------------------------------------------------------------------------!
    738739! Description:
    739740! ------------
    740741!> Subroutine to initialize and set all boundary conditions for chemical species
    741742!------------------------------------------------------------------------------!
    742 
    743743 SUBROUTINE chem_boundary_conds( mode )                                           
    744744
     
    755755        ONLY:  bc_h                                                           
    756756
    757     CHARACTER (len=*), INTENT(IN) ::  mode
     757    CHARACTER (LEN=*), INTENT(IN) ::  mode
    758758    INTEGER(iwp) ::  i                            !< grid index x direction.
    759759    INTEGER(iwp) ::  j                            !< grid index y direction.
     
    768768       CASE ( 'init' )
    769769
    770           IF ( bc_cs_b == 'dirichlet' )    THEN
     770          IF ( bc_cs_b == 'dirichlet' )  THEN
    771771             ibc_cs_b = 0
    772772          ELSEIF ( bc_cs_b == 'neumann' )  THEN
     
    774774          ELSE
    775775             message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' 
    776              CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )     !<
     776             CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )
    777777          ENDIF                                                                 
    778778!
    779 !--          Set Integer flags and check for possible erroneous settings for top
    780 !--          boundary condition.
    781           IF ( bc_cs_t == 'dirichlet' )             THEN
     779!--       Set Integer flags and check for possible erroneous settings for top
     780!--       boundary condition.
     781          IF ( bc_cs_t == 'dirichlet' )  THEN
    782782             ibc_cs_t = 0
    783           ELSEIF ( bc_cs_t == 'neumann' )           THEN
     783          ELSEIF ( bc_cs_t == 'neumann' )  THEN
    784784             ibc_cs_t = 1
    785785          ELSEIF ( bc_cs_t == 'initial_gradient' )  THEN
    786786             ibc_cs_t = 2
    787           ELSEIF ( bc_cs_t == 'nested' )            THEN         
     787          ELSEIF ( bc_cs_t == 'nested' )  THEN         
    788788             ibc_cs_t = 3
    789789          ELSE
     
    792792          ENDIF
    793793
    794 
    795794       CASE ( 'set_bc_bottomtop' )                   
     795!
    796796!--          Bottom boundary condtions for chemical species     
    797           DO lsp = 1, nspec                                                     
    798              IF ( ibc_cs_b == 0 ) THEN                   
    799                 DO l = 0, 1
     797          DO  lsp = 1, nspec                                                     
     798             IF ( ibc_cs_b == 0 )  THEN                   
     799                DO  l = 0, 1
     800!
    800801!--                Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e.
    801802!--                the chem_species(nspec)%conc_p value at the topography top (k-1)
     
    805806                   kb = MERGE( -1, 1, l == 0 )
    806807                   !$OMP PARALLEL DO PRIVATE( i, j, k )
    807                    DO m = 1, bc_h(l)%ns
    808                       i = bc_h(l)%i(m)           
    809                       j = bc_h(l)%j(m)
    810                       k = bc_h(l)%k(m)
     808                   DO  m = 1, bc_h(l)%ns
     809                       i = bc_h(l)%i(m)           
     810                       j = bc_h(l)%j(m)
     811                       k = bc_h(l)%k(m)
    811812                      chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i)
    812813                   ENDDO                                       
    813814                ENDDO                                       
    814815
    815              ELSEIF ( ibc_cs_b == 1 ) THEN
     816             ELSEIF ( ibc_cs_b == 1 )  THEN
     817!
    816818!--             in boundary_conds there is som extra loop over m here for passive tracer
    817                 DO l = 0, 1
     819                DO  l = 0, 1
    818820                   kb = MERGE( -1, 1, l == 0 )
    819821                   !$OMP PARALLEL DO PRIVATE( i, j, k )                                           
     
    828830             ENDIF
    829831       ENDDO    ! end lsp loop 
    830 
     832!
    831833!--    Top boundary conditions for chemical species - Should this not be done for all species?
    832834          IF ( ibc_cs_t == 0 )  THEN                   
    833              DO lsp = 1, nspec
     835             DO  lsp = 1, nspec
    834836                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)       
    835837             ENDDO
    836838          ELSEIF ( ibc_cs_t == 1 )  THEN
    837              DO lsp = 1, nspec
     839             DO  lsp = 1, nspec
    838840                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
    839841             ENDDO
    840842          ELSEIF ( ibc_cs_t == 2 )  THEN
    841              DO lsp = 1, nspec
     843             DO  lsp = 1, nspec
    842844                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
    843845             ENDDO
    844846          ENDIF
    845 !
     847
    846848       CASE ( 'set_bc_lateral' )                       
     849!
    847850!--             Lateral boundary conditions for chem species at inflow boundary
    848851!--             are automatically set when chem_species concentration is
     
    854857
    855858          IF ( bc_radiation_s )  THEN
    856              DO lsp = 1, nspec
     859             DO  lsp = 1, nspec
    857860                chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
    858861             ENDDO
    859862          ELSEIF ( bc_radiation_n )  THEN
    860              DO lsp = 1, nspec
     863             DO  lsp = 1, nspec
    861864                chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
    862865             ENDDO
    863866          ELSEIF ( bc_radiation_l )  THEN
    864              DO lsp = 1, nspec
     867             DO  lsp = 1, nspec
    865868                chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
    866869             ENDDO
    867870          ELSEIF ( bc_radiation_r )  THEN
    868              DO lsp = 1, nspec
     871             DO  lsp = 1, nspec
    869872                chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
    870873             ENDDO
     
    875878 END SUBROUTINE chem_boundary_conds
    876879
    877 !
     880
    878881!------------------------------------------------------------------------------!
    879882! Description:
     
    881884!> Boundary conditions for prognostic variables in chemistry: decycling in the
    882885!> x-direction
    883 !-----------------------------------------------------------------------------
     886!------------------------------------------------------------------------------!
    884887 SUBROUTINE chem_boundary_conds_decycle( cs_3d, cs_pr_init )
     888!
     889!-- Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
     890!-- approproate for chemicals compounds since they may accumulate too much.
     891!-- If no proper boundary conditions from a DYNAMIC input file are available,
     892!-- de-cycling applies the initial profiles at the inflow boundaries for
     893!-- Dirichlet boundary conditions
    885894
    886895    IMPLICIT NONE
     
    897906
    898907    flag = 0.0_wp
    899 
    900 
     908!
    901909!-- Left and right boundaries
    902910    IF ( decycle_chem_lr  .AND.  bc_lr_cyc )  THEN
     
    962970       ENDDO
    963971    ENDIF
    964 
    965 
     972!
    966973!-- South and north boundaries
    967974    IF ( decycle_chem_ns  .AND.  bc_ns_cyc )  THEN
     
    10291036    ENDIF
    10301037 END SUBROUTINE chem_boundary_conds_decycle
    1031    !
    1032 !------------------------------------------------------------------------------!
    1033 !
     1038
     1039
     1040!------------------------------------------------------------------------------!
    10341041! Description:
    10351042! ------------
    10361043!> Subroutine for checking data output for chemical species
    10371044!------------------------------------------------------------------------------!
    1038 
    10391045 SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
    1040 
    10411046
    10421047    IMPLICIT NONE
     
    10501055    INTEGER(iwp) ::  k
    10511056
    1052     CHARACTER(len=16)    ::  spec_name
     1057    CHARACTER(LEN=16)    ::  spec_name
    10531058
    10541059!
     
    10581063    unit = 'illegal'
    10591064
    1060     spec_name = TRIM(var(4:))             !< var 1:3 is 'kc_' or 'em_'.
    1061 
    1062     IF ( TRIM(var(1:3)) == 'em_' )  THEN
    1063        DO lsp=1,nspec
    1064           IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
     1065    spec_name = TRIM( var(4:) )             !< var 1:3 is 'kc_' or 'em_'.
     1066
     1067    IF ( TRIM( var(1:3) ) == 'em_' )  THEN
     1068       DO  lsp=1,nspec
     1069          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
    10651070             unit = 'mol m-2 s-1'
    10661071          ENDIF
    1067           !
    1068           !-- It is possible to plant PM10 and PM25 into the gasphase chemistry code
    1069           !-- as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
    1070           !-- set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
    1071           IF (spec_name(1:2) == 'PM')   THEN
     1072!
     1073!--      It is possible to plant PM10 and PM25 into the gasphase chemistry code
     1074!--      as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
     1075!--      set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
     1076          IF (spec_name(1:2) == 'PM')  THEN
    10721077             unit = 'kg m-2 s-1'
    10731078          ENDIF
     
    10761081    ELSE
    10771082
    1078        DO lsp=1,nspec
    1079           IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
     1083       DO  lsp=1,nspec
     1084          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
    10801085             unit = 'ppm'
    10811086          ENDIF
     
    10841089!--            as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
    10851090!--            set unit to kilograms per m**3 for PM10 and PM25 (PM2.5)
    1086           IF (spec_name(1:2) == 'PM')   THEN 
     1091          IF (spec_name(1:2) == 'PM')  THEN 
    10871092            unit = 'kg m-3'
    10881093          ENDIF
    10891094       ENDDO
    10901095
    1091        DO lsp=1,nphot
    1092           IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name))   THEN
     1096       DO  lsp=1,nphot
     1097          IF (TRIM( spec_name ) == TRIM( phot_frequen(lsp)%name ) )  THEN
    10931098             unit = 'sec-1'
    10941099          ENDIF
     
    10991104    RETURN
    11001105 END SUBROUTINE chem_check_data_output
    1101 !
    1102 !------------------------------------------------------------------------------!
    1103 !
     1106
     1107
     1108!------------------------------------------------------------------------------!
    11041109! Description:
    11051110! ------------
    11061111!> Subroutine for checking data output of profiles for chemistry model
    11071112!------------------------------------------------------------------------------!
    1108 
    11091113 SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
    11101114
     
    11221126    CHARACTER (LEN=*) ::  variable !<
    11231127    CHARACTER (LEN=*) ::  dopr_unit
    1124     CHARACTER(len=16) ::  spec_name
     1128    CHARACTER (LEN=16) ::  spec_name
    11251129
    11261130    INTEGER(iwp) ::  var_count, lsp  !<
    11271131
    11281132
    1129     spec_name = TRIM(variable(4:))   
     1133    spec_name = TRIM( variable(4:) )   
    11301134
    11311135       IF (  .NOT.  air_chemistry )  THEN
     
    11361140
    11371141       ELSE
    1138           DO lsp = 1, nspec
    1139              IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN
     1142          DO  lsp = 1, nspec
     1143             IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
    11401144                cs_pr_count = cs_pr_count+1
    11411145                cs_pr_index(cs_pr_count) = lsp
    11421146                dopr_index(var_count) = pr_palm+cs_pr_count 
    11431147                dopr_unit  = 'ppm'
    1144                 IF (spec_name(1:2) == 'PM')   THEN
     1148                IF (spec_name(1:2) == 'PM')  THEN
    11451149                     dopr_unit = 'kg m-3'
    11461150                ENDIF
     
    11531157 END SUBROUTINE chem_check_data_output_pr
    11541158
    1155 !
     1159
    11561160!------------------------------------------------------------------------------!
    11571161! Description:
     
    11631167    IMPLICIT NONE
    11641168
    1165     LOGICAL       ::  found
     1169    LOGICAL  ::  found
    11661170    INTEGER (iwp) ::  lsp_usr      !< running index for user defined chem spcs
    11671171    INTEGER (iwp) ::  lsp          !< running index for chem spcs.
     
    11711175       message_string = 'Chemical reactions: ON'
    11721176       CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 )
    1173     ELSEIF ( .not. (chem_gasphase_on) ) THEN
     1177    ELSEIF ( .NOT. (chem_gasphase_on) ) THEN
    11741178       message_string = 'Chemical reactions: OFF'
    11751179       CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 )
    11761180    ENDIF
    1177 
     1181!
    11781182!-- check for chemistry time-step
    11791183    IF ( call_chem_at_all_substeps )  THEN
    11801184       message_string = 'Chemistry is calculated at all meteorology time-step'
    11811185       CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 )
    1182     ELSEIF ( .not. (call_chem_at_all_substeps) ) THEN
     1186    ELSEIF ( .not. (call_chem_at_all_substeps) )  THEN
    11831187       message_string = 'Sub-time-steps are skipped for chemistry time-steps'
    11841188       CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 )
    11851189    ENDIF
    1186 
     1190!
    11871191!-- check for photolysis scheme
    11881192    IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant')  )  THEN
     
    11901194       CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 )
    11911195    ENDIF
    1192 
     1196!
    11931197!-- check for  decycling of chem species
    1194     IF ( (.not. any(decycle_method == 'neumann') ) .AND. (.not. any(decycle_method == 'dirichlet') ) )   THEN
     1198    IF ( (.NOT. any(decycle_method == 'neumann') ) .AND. (.NOT. any(decycle_method == 'dirichlet') ) )  THEN
    11951199       message_string = 'Incorrect boundary conditions. Only neumann or '   &
    11961200                // 'dirichlet &available for decycling chemical species '
    11971201       CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
    11981202    ENDIF
     1203!
    11991204!-- check for chemical mechanism used
    12001205    CALL get_mechanism_name
    1201     IF (chem_mechanism /= trim(cs_mech) )  THEN
     1206    IF ( chem_mechanism /= TRIM( cs_mech ) )  THEN
    12021207       message_string = 'Incorrect chemistry mechanism selected, check spelling in namelist and/or chem_gasphase_mod'
    12031208       CALL message( 'chem_check_parameters', 'CM0462', 1, 2, 0, 6, 0 )
    12041209    ENDIF
    1205 
    1206 
    1207 !---------------------
     1210!
    12081211!-- chem_check_parameters is called before the array chem_species is allocated!
    12091212!-- temporary switch of this part of the check
    12101213!    RETURN                !bK commented
     1214
    12111215    CALL chem_init_internal
    1212 !---------------------
    1213 
     1216!
    12141217!-- check for initial chem species input
    12151218    lsp_usr = 1
     
    12171220    DO WHILE ( cs_name (lsp_usr) /= 'novalue')
    12181221       found = .FALSE.
    1219        DO lsp = 1, nvar
    1220           IF ( TRIM(cs_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
     1222       DO  lsp = 1, nvar
     1223          IF ( TRIM( cs_name (lsp_usr) ) == TRIM( chem_species(lsp)%name) ) THEN
    12211224             found = .TRUE.
    12221225             EXIT
    12231226          ENDIF
    12241227       ENDDO
    1225        IF ( .not. found ) THEN
    1226           message_string = 'Unused/incorrect input for initial surface value: ' // TRIM(cs_name(lsp_usr))
     1228       IF ( .NOT.  found )  THEN
     1229          message_string = 'Unused/incorrect input for initial surface value: ' //     &
     1230                            TRIM( cs_name(lsp_usr) )
    12271231          CALL message( 'chem_check_parameters', 'CM0427', 1, 2, 0, 6, 0 )
    12281232       ENDIF
    12291233       lsp_usr = lsp_usr + 1
    12301234    ENDDO
    1231 
     1235!
    12321236!-- check for surface  emission flux chem species
    1233 
    12341237    lsp_usr = 1
    12351238    lsp     = 1
    12361239    DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue')
    12371240       found = .FALSE.
    1238        DO lsp = 1, nvar
    1239           IF ( TRIM(surface_csflux_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
     1241       DO  lsp = 1, nvar
     1242          IF ( TRIM( surface_csflux_name (lsp_usr) ) == TRIM( chem_species(lsp)%name ) ) THEN
    12401243             found = .TRUE.
    12411244             EXIT
    12421245          ENDIF
    12431246       ENDDO
    1244        IF ( .not. found ) THEN
     1247       IF ( .NOT.  found ) THEN
    12451248          message_string = 'Unused/incorrect input of chemical species for surface emission fluxes: '  &
    1246                             // TRIM(surface_csflux_name(lsp_usr))
     1249                            // TRIM( surface_csflux_name(lsp_usr) )
    12471250          CALL message( 'chem_check_parameters', 'CM0428', 1, 2, 0, 6, 0 )
    12481251       ENDIF
     
    12521255 END SUBROUTINE chem_check_parameters
    12531256
    1254 !
    1255 !------------------------------------------------------------------------------!
    1256 !
     1257
     1258!------------------------------------------------------------------------------!
    12571259! Description:
    12581260! ------------
     
    12601262!> @todo: Remove "mode" from argument list, not used.
    12611263!------------------------------------------------------------------------------!
    1262    
    12631264 SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf,   &
    12641265                                  two_d, nzb_do, nzt_do, fill_value )
     
    12811282    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf 
    12821283
    1283     !
    1284     !-- local variables.
    1285     CHARACTER(len=16)    ::  spec_name
     1284!
     1285!-- local variables.
     1286    CHARACTER(LEN=16)    ::  spec_name
    12861287    INTEGER(iwp) ::  lsp
    12871288    INTEGER(iwp) ::  i               !< grid index along x-direction
     
    12891290    INTEGER(iwp) ::  k               !< grid index along z-direction
    12901291    INTEGER(iwp) ::  char_len        !< length of a character string
    1291 
    12921292!
    12931293!-- Next statement is to avoid compiler warnings about unused variables
     
    12951295
    12961296    found = .FALSE.
    1297     char_len  = LEN_TRIM(variable)
     1297    char_len  = LEN_TRIM( variable )
    12981298
    12991299    spec_name = TRIM( variable(4:char_len-3) )
    13001300
    1301        DO lsp=1,nspec
    1302           IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)    .AND.                             &
    1303                 ( (variable(char_len-2:) == '_xy')               .OR.                              &
    1304                   (variable(char_len-2:) == '_xz')               .OR.                              &
    1305                   (variable(char_len-2:) == '_yz') ) )               THEN             
     1301       DO  lsp=1,nspec
     1302          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name )  .AND.                           &
     1303                ( (variable(char_len-2:) == '_xy')  .OR.                                           &
     1304                  (variable(char_len-2:) == '_xz')  .OR.                                           &
     1305                  (variable(char_len-2:) == '_yz') ) )  THEN             
    13061306!
    13071307!--   todo: remove or replace by "CALL message" mechanism (kanani)
    1308 !                    IF(myid == 0)  WRITE(6,*) 'Output of species ' // TRIM(variable)  //               &
    1309 !                                                             TRIM(chem_species(lsp)%name)       
    1310              IF (av == 0) THEN
     1308!                    IF(myid == 0)  WRITE(6,*) 'Output of species ' // TRIM( variable )  //       &
     1309!                                                             TRIM( chem_species(lsp)%name )       
     1310             IF (av == 0)  THEN
    13111311                DO  i = nxl, nxr
    13121312                   DO  j = nys, nyn
     
    13161316                                              REAL( fill_value, KIND = wp ),                       &
    13171317                                              BTEST( wall_flags_0(k,j,i), 0 ) )
    1318 
    1319 
    13201318                      ENDDO
    13211319                   ENDDO
     
    13341332                ENDDO
    13351333             ENDIF
    1336               grid = 'zu'           
    1337               found = .TRUE.
     1334             grid = 'zu'           
     1335             found = .TRUE.
    13381336          ENDIF
    13391337       ENDDO
     
    13431341 END SUBROUTINE chem_data_output_2d     
    13441342
    1345 !
    1346 !------------------------------------------------------------------------------!
    1347 !
     1343
     1344!------------------------------------------------------------------------------!
    13481345! Description:
    13491346! ------------
    13501347!> Subroutine defining 3D output variables for chemical species
    13511348!------------------------------------------------------------------------------!
    1352 
    13531349 SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
    1354 
    13551350
    13561351    USE indices
     
    13721367    REAL(wp)             ::  fill_value   !<
    13731368    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf
    1374 
    1375 
    1376     !-- local variables
    1377     CHARACTER(len=16)    ::  spec_name
     1369!
     1370!-- local variables
     1371    CHARACTER(LEN=16)    ::  spec_name
    13781372    INTEGER(iwp)         ::  i
    13791373    INTEGER(iwp)         ::  j
     
    13851379
    13861380    found = .FALSE.
    1387     IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) ) THEN
     1381    IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) )  THEN
    13881382       RETURN
    13891383    ENDIF
    13901384
    1391     spec_name = TRIM(variable(4:))
    1392 
    1393     IF ( variable(1:3) == 'em_' ) THEN
    1394 
    1395       local_pf = 0.0_wp
    1396 
    1397       DO lsp = 1, nvar   !!! cssws - nvar species, chem_species - nspec species !!!
    1398          IF ( TRIM(spec_name) == TRIM(chem_species(lsp)%name) )   THEN
    1399             !
    1400             !-- no average for now
    1401             DO m = 1, surf_usm_h%ns
    1402                local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = &
    1403                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m)
    1404             ENDDO
    1405             DO m = 1, surf_lsm_h%ns
    1406                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = &
     1385    spec_name = TRIM( variable(4:) )
     1386
     1387    IF ( variable(1:3) == 'em_' )  THEN
     1388
     1389       local_pf = 0.0_wp
     1390
     1391       DO lsp = 1, nvar   !!! cssws - nvar species, chem_species - nspec species !!!
     1392          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
     1393!
     1394!--          no average for now
     1395             DO m = 1, surf_usm_h%ns
     1396                local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = &
     1397                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m)
     1398             ENDDO
     1399             DO m = 1, surf_lsm_h%ns
     1400                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = &
    14071401                  local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m)
    1408             ENDDO
    1409             DO l = 0, 3
    1410                DO m = 1, surf_usm_v(l)%ns
    1411                   local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = &
     1402             ENDDO
     1403             DO l = 0, 3
     1404                DO m = 1, surf_usm_v(l)%ns
     1405                   local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = &
    14121406                     local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) + surf_usm_v(l)%cssws(lsp,m)
    1413                ENDDO
    1414                DO m = 1, surf_lsm_v(l)%ns
    1415                   local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = &
    1416                      local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) + surf_lsm_v(l)%cssws(lsp,m)
    1417                ENDDO
    1418             ENDDO
    1419             found = .TRUE.
    1420          ENDIF
    1421       ENDDO
     1407                ENDDO
     1408                DO m = 1, surf_lsm_v(l)%ns
     1409                   local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = &
     1410                      local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) + surf_lsm_v(l)%cssws(lsp,m)
     1411                ENDDO
     1412             ENDDO
     1413             found = .TRUE.
     1414          ENDIF
     1415       ENDDO
    14221416    ELSE
    1423       DO lsp=1,nspec
    1424          IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
     1417      DO  lsp = 1, nspec
     1418         IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
    14251419!
    14261420!--   todo: remove or replace by "CALL message" mechanism (kanani)
    1427 !              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM(variable)  // &
    1428 !                                                           TRIM(chem_species(lsp)%name)       
    1429             IF (av == 0) THEN
     1421!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
     1422!                                                           TRIM( chem_species(lsp)%name )       
     1423            IF (av == 0)  THEN
    14301424               DO  i = nxl, nxr
    14311425                  DO  j = nys, nyn
     
    14401434
    14411435            ELSE
     1436
    14421437               DO  i = nxl, nxr
    14431438                  DO  j = nys, nyn
     
    14591454
    14601455 END SUBROUTINE chem_data_output_3d
    1461 !
    1462 !------------------------------------------------------------------------------!
    1463 !
     1456
     1457
     1458!------------------------------------------------------------------------------!
    14641459! Description:
    14651460! ------------
    14661461!> Subroutine defining mask output variables for chemical species
    14671462!------------------------------------------------------------------------------!
    1468    
    14691463 SUBROUTINE chem_data_output_mask( av, variable, found, local_pf )
    14701464
     
    14771471    IMPLICIT NONE
    14781472
    1479     CHARACTER(LEN=5) ::  grid        !< flag to distinquish between staggered grids
    1480 
    1481     CHARACTER (LEN=*)::  variable    !<
    1482     INTEGER(iwp)     ::  av          !< flag to control data output of instantaneous or time-averaged data
    1483     LOGICAL          ::  found       !<
     1473    CHARACTER(LEN=5)  ::  grid        !< flag to distinquish between staggered grids
     1474    CHARACTER(LEN=*)  ::  variable    !<
     1475    INTEGER(iwp)  ::  av              !< flag to control data output of instantaneous or time-averaged data
     1476    LOGICAL ::  found
    14841477    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
    14851478              local_pf   !<
    1486 
    1487 
    1488     !-- local variables.
    1489     CHARACTER(len=16)    ::  spec_name
     1479!
     1480!-- local variables.
     1481    CHARACTER(LEN=16)  ::  spec_name
    14901482    INTEGER(iwp) ::  lsp
    14911483    INTEGER(iwp) ::  i               !< grid index along x-direction
     
    14951487
    14961488    found = .TRUE.
    1497     grid  = 's'
     1489    grid = 's'
    14981490
    14991491    spec_name = TRIM( variable(4:) )
    15001492
    1501     DO lsp=1,nspec
    1502        IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name) )               THEN             
     1493    DO  lsp=1,nspec
     1494       IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN             
    15031495!
    15041496!-- todo: remove or replace by "CALL message" mechanism (kanani)
    1505 !              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM(variable)  // &
    1506 !                                                        TRIM(chem_species(lsp)%name)       
    1507           IF (av == 0) THEN
     1497!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
     1498!                                                        TRIM( chem_species(lsp)%name )       
     1499          IF (av == 0)  THEN
    15081500             IF ( .NOT. mask_surface(mid) )  THEN
    15091501
     
    15821574             ENDIF
    15831575
    1584 
    15851576          ENDIF
    15861577          found = .FALSE.
     
    15921583 END SUBROUTINE chem_data_output_mask     
    15931584
    1594 !
    1595 !------------------------------------------------------------------------------!
    1596 !
     1585
     1586!------------------------------------------------------------------------------!
    15971587! Description:
    15981588! ------------
     
    16121602    found  = .TRUE.
    16131603
    1614     IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )   THEN                   !< always the same grid for chemistry variables
     1604    IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )  THEN                   !< always the same grid for chemistry variables
    16151605       grid_x = 'x'
    16161606       grid_y = 'y'
     
    16251615
    16261616 END SUBROUTINE chem_define_netcdf_grid
    1627 !
    1628 !------------------------------------------------------------------------------!
    1629 !
     1617
     1618
     1619!------------------------------------------------------------------------------!
    16301620! Description:
    16311621! ------------
     
    16371627
    16381628    INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
    1639     INTEGER(iwp)             :: lsp            !< running index for chem spcs
    1640     INTEGER(iwp)             :: cs_fixed
    1641     CHARACTER (LEN=80)       :: docsflux_chr
    1642     CHARACTER (LEN=80)       :: docsinit_chr
     1629    INTEGER(iwp)  :: lsp                       !< running index for chem spcs
     1630    INTEGER(iwp)  :: cs_fixed
     1631    CHARACTER (LEN=80)  :: docsflux_chr
     1632    CHARACTER (LEN=80)  :: docsinit_chr
    16431633!
    16441634! Get name of chemical mechanism from chem_gasphase_mod
    16451635    CALL get_mechanism_name
    1646 
     1636!
    16471637!-- Write chemistry model  header
    16481638    WRITE( io, 1 )
    1649 
     1639!
    16501640!-- Gasphase reaction status
    16511641    IF ( chem_gasphase_on )  THEN
     
    16541644       WRITE( io, 3 )
    16551645    ENDIF
    1656 
    16571646!
    16581647!-- Chemistry time-step
    16591648    WRITE ( io, 4 ) cs_time_step
    1660 
     1649!
    16611650!-- Emission mode info
    16621651    IF ( mode_emis == "DEFAULT" )  THEN
     
    16671656       WRITE( io, 7 )
    16681657    ENDIF
    1669 
     1658!
    16701659!-- Photolysis scheme info
    1671     IF ( photolysis_scheme == "simple" )      THEN
     1660    IF ( photolysis_scheme == "simple" )  THEN
    16721661       WRITE( io, 8 )
    1673     ELSEIF (photolysis_scheme == "constant" ) THEN
     1662    ELSEIF (photolysis_scheme == "constant" )  THEN
    16741663       WRITE( io, 9 )
    16751664    ENDIF
    1676 
     1665!
    16771666!-- Emission flux info
    16781667    lsp = 1
     
    16811670       docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ',' 
    16821671       IF ( LEN_TRIM( docsflux_chr ) >= 75 )  THEN
    1683         WRITE ( io, 10 )  docsflux_chr
    1684         docsflux_chr = '       '
     1672          WRITE ( io, 10 )  docsflux_chr
     1673          docsflux_chr = '       '
    16851674       ENDIF
    16861675       lsp = lsp + 1
     
    16901679       WRITE ( io, 10 )  docsflux_chr
    16911680    ENDIF
    1692 
    1693 
     1681!
    16941682!-- initializatoin of Surface and profile chemical species
    16951683
     
    17081696       WRITE ( io, 11 )  docsinit_chr
    17091697    ENDIF
    1710 
     1698!
    17111699!-- number of variable and fix chemical species and number of reactions
    17121700    cs_fixed = nspec - nvar
     
    17341722 END SUBROUTINE chem_header
    17351723
    1736 !
    1737 !------------------------------------------------------------------------------!
    1738 !
     1724
     1725!------------------------------------------------------------------------------!
    17391726! Description:
    17401727! ------------
     
    17421729!------------------------------------------------------------------------------!
    17431730 SUBROUTINE chem_init_arrays
    1744 
     1731!
    17451732!-- Please use this place to allocate required arrays
    17461733
    17471734 END SUBROUTINE chem_init_arrays
    17481735
    1749 !
    1750 !------------------------------------------------------------------------------!
    1751 !
     1736
     1737!------------------------------------------------------------------------------!
    17521738! Description:
    17531739! ------------
     
    17671753    INTEGER(iwp) ::  j !< running index y dimension
    17681754    INTEGER(iwp) ::  n !< running index for chemical species
    1769 
    17701755!
    17711756!-- Next statement is to avoid compiler warning about unused variables
     
    17741759           ilu_urban ) == 0 )  CONTINUE
    17751760         
    1776 
    1777 
    17781761    IF ( emissions_anthropogenic )  CALL chem_emissions_init
    17791762!
     
    17961779 END SUBROUTINE chem_init
    17971780
    1798 !
    1799 !------------------------------------------------------------------------------!
    1800 !
     1781
     1782!------------------------------------------------------------------------------!
    18011783! Description:
    18021784! ------------
     
    18131795
    18141796    IMPLICIT NONE
    1815 
    18161797!
    18171798!-- Local variables
     
    18211802    INTEGER(iwp) ::  lpr_lev           !< running index for chem spcs profile level
    18221803
    1823     IF ( emissions_anthropogenic ) THEN
     1804    IF ( emissions_anthropogenic )  THEN
    18241805       CALL netcdf_data_input_chemistry_data( chem_emis_att, chem_emis )
    18251806    ENDIF
     
    18741855
    18751856    ENDDO
    1876 
    18771857!
    18781858!-- Initial concentration of profiles is prescribed by parameters cs_profile
    18791859!-- and cs_heights in the namelist &chemistry_parameters
     1860
    18801861    CALL chem_init_profiles
    18811862!   
     
    18921873       ENDDO
    18931874    ENDIF
    1894 
    1895 
    1896 
    18971875!
    18981876!-- Initialize model variables
    18991877    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
    19001878         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
    1901 
    1902 
     1879!
    19031880!--    First model run of a possible job queue.
    19041881!--    Initial profiles of the variables must be computed.
     
    19071884!
    19081885!--       Transfer initial profiles to the arrays of the 3D model
    1909           DO lsp = 1, nspec
     1886          DO  lsp = 1, nspec
    19101887             DO  i = nxlg, nxrg
    19111888                DO  j = nysg, nyng
     
    19211898          CALL location_message( 'initializing with constant chemistry profiles', .FALSE. )
    19221899
    1923           DO lsp = 1, nspec
     1900          DO  lsp = 1, nspec
    19241901             DO  i = nxlg, nxrg
    19251902                DO  j = nysg, nyng
     
    19301907
    19311908       ENDIF
    1932 
    19331909!
    19341910!--    If required, change the surface chem spcs at the start of the 3D run
    1935        IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN           
    1936           DO lsp = 1, nspec
     1911       IF ( cs_surface_initial_change(1) /= 0.0_wp )  THEN           
     1912          DO  lsp = 1, nspec
    19371913             chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) +  &
    19381914                                               cs_surface_initial_change(lsp)
     
    19411917!
    19421918!--    Initiale old and new time levels.
    1943        DO lsp = 1, nvar
     1919       DO  lsp = 1, nvar
    19441920          chem_species(lsp)%tconc_m = 0.0_wp                     
    19451921          chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
     
    19481924    ENDIF
    19491925
    1950 
    1951 
    1952 !--- new code add above this line
    1953     DO lsp = 1, nphot
     1926    DO  lsp = 1, nphot
    19541927       phot_frequen(lsp)%name = phot_names(lsp)
    19551928!
    19561929!-- todo: remove or replace by "CALL message" mechanism (kanani)
    1957 !         IF( myid == 0 )  THEN
    1958 !            WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM(phot_names(lsp))
    1959 !         ENDIF
    1960           phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => freq_1(:,:,:,lsp)
     1930!--       IF( myid == 0 )  THEN
     1931!--          WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM( phot_names(lsp) )
     1932!--       ENDIF
     1933          phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => freq_1(:,:,:,lsp)
    19611934    ENDDO
    19621935
     
    19671940 END SUBROUTINE chem_init_internal
    19681941
    1969 !
    1970 !------------------------------------------------------------------------------!
    1971 !
     1942
     1943!------------------------------------------------------------------------------!
    19721944! Description:
    19731945! ------------
     
    19761948!> analogue to parameters u_profile, v_profile and uv_heights)
    19771949!------------------------------------------------------------------------------!
    1978  SUBROUTINE chem_init_profiles              !< SUBROUTINE is called from chem_init in case of
    1979    !< TRIM( initializing_actions ) /= 'read_restart_data'
    1980    !< We still need to see what has to be done in case of restart run
     1950 SUBROUTINE chem_init_profiles             
     1951!
     1952!-- SUBROUTINE is called from chem_init in case of TRIM( initializing_actions ) /= 'read_restart_data'
     1953!< We still need to see what has to be done in case of restart run
     1954
    19811955    USE chem_modules
    19821956
    19831957    IMPLICIT NONE
    1984 
    1985     !-- Local variables
     1958!
     1959!-- Local variables
    19861960    INTEGER ::  lsp        !< running index for number of species in derived data type species_def
    19871961    INTEGER ::  lsp_usr     !< running index for number of species (user defined)  in cs_names, cs_profiles etc
    19881962    INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs.
    19891963    INTEGER ::  npr_lev    !< the next available profile lev
    1990 
    1991     !-----------------
    1992     !-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
    1993     !-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
    1994     !-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
    1995     !-- "cs_heights" are prescribed, their values will!override the constant profile given by
    1996     !-- "cs_surface".
    1997     !
     1964!
     1965!-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
     1966!-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
     1967!-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
     1968!-- "cs_heights" are prescribed, their values will!override the constant profile given by
     1969!-- "cs_surface".
    19981970    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    19991971       lsp_usr = 1
    20001972       DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
    20011973          DO  lsp = 1, nspec                                !
    2002              !-- create initial profile (conc_pr_init) for each chemical species
     1974!
     1975!--          create initial profile (conc_pr_init) for each chemical species
    20031976             IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
    2004                 IF ( cs_profile(lsp_usr,1) == 9999999.9_wp ) THEN
    2005                    !-- set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
     1977                IF ( cs_profile(lsp_usr,1) == 9999999.9_wp )  THEN
     1978!
     1979!--               set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
    20061980                   DO lpr_lev = 0, nzt+1
    20071981                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
     
    20161990
    20171991                   npr_lev = 1
    2018                    !                     chem_species(lsp)%conc_pr_init(0) = 0.0_wp
     1992!                   chem_species(lsp)%conc_pr_init(0) = 0.0_wp
    20191993                   DO  lpr_lev = 1, nz+1
    20201994                      IF ( npr_lev < 100 )  THEN
     
    20282002                         ENDDO
    20292003                      ENDIF
    2030                       IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr, npr_lev + 1) /= 9999999.9_wp )  THEN
     2004                      IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp )  THEN
    20312005                         chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +       &
    20322006                              ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
     
    20382012                   ENDDO
    20392013                ENDIF
    2040                 !-- If a profile is prescribed explicity using cs_profiles and cs_heights, then 
    2041                 !-- chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
    2042                 !-- on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
     2014!
     2015!--          If a profile is prescribed explicity using cs_profiles and cs_heights, then 
     2016!--          chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
     2017!--          on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
    20432018             ENDIF
    20442019          ENDDO
     
    20482023
    20492024 END SUBROUTINE chem_init_profiles
    2050  !
    2051  !------------------------------------------------------------------------------!
    2052  !
    2053  ! Description:
    2054  ! ------------
    2055  !> Subroutine to integrate chemical species in the given chemical mechanism
    2056  !------------------------------------------------------------------------------!
    2057 
     2025
     2026 
     2027!------------------------------------------------------------------------------!
     2028! Description:
     2029! ------------
     2030!> Subroutine to integrate chemical species in the given chemical mechanism
     2031!------------------------------------------------------------------------------!
    20582032 SUBROUTINE chem_integrate_ij( i, j )
    20592033
     
    20682042    INTEGER,INTENT(IN)       :: i
    20692043    INTEGER,INTENT(IN)       :: j
    2070 
    2071     !--   local variables
     2044!
     2045!--   local variables
    20722046    INTEGER(iwp) ::  lsp                                                     !< running index for chem spcs.
    20732047    INTEGER(iwp) ::  lph                                                     !< running index for photolysis frequencies
     
    20962070    REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
    20972071
    2098 
    20992072    REAL(kind=wp)  :: dt_chem                                             
    2100 
    21012073!
    21022074!-- Set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
    2103     IF (chem_gasphase_on) THEN
     2075    IF (chem_gasphase_on)  THEN
    21042076       nacc = 0
    21052077       nrej = 0
    21062078
    21072079       tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
    2108        !
    2109        !--       ppm to molecules/cm**3
    2110        !--       tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp * hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp ) 
     2080!
     2081!--    convert ppm to molecules/cm**3
     2082!--    tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp *
     2083!--               hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp ) 
    21112084       conv = ppm2fr * xna / vmolcm
    21122085       tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
    21132086       tmp_fact_i = 1.0_wp/tmp_fact
    21142087
    2115        IF ( humidity ) THEN
     2088       IF ( humidity )  THEN
    21162089          IF ( bulk_cloud_model )  THEN
    21172090             tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) *  &
     
    21242097       ENDIF
    21252098
    2126        DO lsp = 1,nspec
     2099       DO  lsp = 1,nspec
    21272100          tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:)
    21282101       ENDDO
     
    21312104          tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)               
    21322105       ENDDO
    2133        !
    2134        !--   todo: remove (kanani)
    2135        !           IF(myid == 0 .AND. chem_debug0 ) THEN
    2136        !              IF (i == 10 .AND. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d
    2137        !           ENDIF
    2138 
    2139        !--       Compute length of time step
     2106!
     2107!--    Compute length of time step
    21402108       IF ( call_chem_at_all_substeps )  THEN
    21412109          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
     
    21462114       cs_time_step = dt_chem
    21472115
    2148        IF(maxval(rcntrl) > 0.0)   THEN    ! Only if rcntrl is set
     2116       IF(maxval(rcntrl) > 0.0)  THEN    ! Only if rcntrl is set
    21492117          IF( time_since_reference_point <= 2*dt_3d)  THEN
    21502118             rcntrl_local = 0
     
    21592127            icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus )
    21602128
    2161        DO lsp = 1,nspec
     2129       DO  lsp = 1,nspec
    21622130          chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)
    21632131       ENDDO
     
    21682136    RETURN
    21692137 END SUBROUTINE chem_integrate_ij
    2170  !
    2171  !------------------------------------------------------------------------------!
    2172  !
    2173  ! Description:
    2174  ! ------------
    2175  !> Subroutine defining parin for &chemistry_parameters for chemistry model
    2176  !------------------------------------------------------------------------------!
     2138
     2139
     2140!------------------------------------------------------------------------------!
     2141! Description:
     2142! ------------
     2143!> Subroutine defining parin for &chemistry_parameters for chemistry model
     2144!------------------------------------------------------------------------------!
    21772145 SUBROUTINE chem_parin
    21782146
     
    22302198         surface_csflux_name,              &
    22312199         time_fac_type
    2232 
    2233     !-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
    2234     !-- so this way we could prescribe a specific flux value for each species
     2200!
     2201!-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
     2202!-- so this way we could prescribe a specific flux value for each species
    22352203    !>  chemistry_parameters for initial profiles
    22362204    !>  cs_names = 'O3', 'NO2', 'NO', ...   to set initial profiles)
     
    22402208    !>  If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
    22412209    !>  then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
    2242 
    2243     !
    2244     !--   Read chem namelist   
     2210!
     2211!--   Read chem namelist   
    22452212
    22462213    CHARACTER(LEN=8)    :: solver_type
     
    22522219    atol = 1.0_wp
    22532220    rtol = 0.01_wp
    2254     !
    2255     !--   Try to find chemistry package
     2221!
     2222!--   Try to find chemistry package
    22562223    REWIND ( 11 )
    22572224    line = ' '
     
    22602227    ENDDO
    22612228    BACKSPACE ( 11 )
    2262     !
    2263     !--   Read chemistry namelist
     2229!
     2230!--   Read chemistry namelist
    22642231    READ ( 11, chemistry_parameters, ERR = 10, END = 20 )     
    2265     !
    2266     !--   Enable chemistry model
     2232!
     2233!--   Enable chemistry model
    22672234    air_chemistry = .TRUE.                   
    22682235    GOTO 20
     
    22732240
    22742241 20 CONTINUE
    2275 
    2276     !
    2277     !--    check for  emission mode for chem species
     2242!
     2243!--    check for emission mode for chem species
    22782244    IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
    22792245       message_string = 'Incorrect mode_emiss  option select. Please check spelling'
     
    22822248
    22832249    t_steps = my_steps         
    2284 
    2285     !--    Determine the number of user-defined profiles and append them to the
    2286     !--    standard data output (data_output_pr)
     2250!
     2251!--    Determine the number of user-defined profiles and append them to the
     2252!--    standard data output (data_output_pr)
    22872253    max_pr_cs_tmp = 0
    22882254    i = 1
    22892255
    22902256    DO  WHILE ( data_output_pr(i)  /= ' '  .AND.  i <= 100 )
    2291        IF ( TRIM( data_output_pr(i)(1:3)) == 'kc_' ) THEN
     2257       IF ( TRIM( data_output_pr(i)(1:3) ) == 'kc_' ) THEN
    22922258          max_pr_cs_tmp = max_pr_cs_tmp+1
    22932259       ENDIF
     
    22952261    ENDDO
    22962262
    2297     IF ( max_pr_cs_tmp > 0 ) THEN
     2263    IF ( max_pr_cs_tmp > 0 )  THEN
    22982264       cs_pr_namelist_found = .TRUE.
    22992265       max_pr_cs = max_pr_cs_tmp
     
    23012267
    23022268    !      Set Solver Type
    2303     IF(icntrl(3) == 0)   THEN
     2269    IF(icntrl(3) == 0)  THEN
    23042270       solver_type = 'rodas3'           !Default
    2305     ELSE IF(icntrl(3) == 1)   THEN
     2271    ELSE IF(icntrl(3) == 1)  THEN
    23062272       solver_type = 'ros2'
    2307     ELSE IF(icntrl(3) == 2)   THEN
     2273    ELSE IF(icntrl(3) == 2)  THEN
    23082274       solver_type = 'ros3'
    2309     ELSE IF(icntrl(3) == 3)   THEN
     2275    ELSE IF(icntrl(3) == 3)  THEN
    23102276       solver_type = 'ro4'
    2311     ELSE IF(icntrl(3) == 4)   THEN
     2277    ELSE IF(icntrl(3) == 4)  THEN
    23122278       solver_type = 'rodas3'
    2313     ELSE IF(icntrl(3) == 5)   THEN
     2279    ELSE IF(icntrl(3) == 5)  THEN
    23142280       solver_type = 'rodas4'
    2315     ELSE IF(icntrl(3) == 6)   THEN
     2281    ELSE IF(icntrl(3) == 6)  THEN
    23162282       solver_type = 'Rang3'
    23172283    ELSE
     
    23202286    END IF
    23212287
    2322     !
    2323     !--   todo: remove or replace by "CALL message" mechanism (kanani)
    2324     !       write(text,*) 'gas_phase chemistry: solver_type = ',TRIM(solver_type)
    2325     !kk    Has to be changed to right calling sequence
    2326     !kk       CALL location_message( TRIM(text), .FALSE. )
    2327     !        IF(myid == 0)   THEN
    2328     !           write(9,*) ' '
    2329     !           write(9,*) 'kpp setup '
    2330     !           write(9,*) ' '
    2331     !           write(9,*) '    gas_phase chemistry: solver_type = ',TRIM(solver_type)
    2332     !           write(9,*) ' '
    2333     !           write(9,*) '    Hstart  = ',rcntrl(3)
    2334     !           write(9,*) '    FacMin  = ',rcntrl(4)
    2335     !           write(9,*) '    FacMax  = ',rcntrl(5)
    2336     !           write(9,*) ' '
    2337     !           IF(vl_dim > 1)    THEN
    2338     !              write(9,*) '    Vector mode                   vektor length = ',vl_dim
    2339     !           ELSE
    2340     !              write(9,*) '    Scalar mode'
    2341     !           ENDIF
    2342     !           write(9,*) ' '
    2343     !        END IF
     2288!
     2289!--   todo: remove or replace by "CALL message" mechanism (kanani)
     2290!       write(text,*) 'gas_phase chemistry: solver_type = ',TRIM( solver_type )
     2291!kk    Has to be changed to right calling sequence
     2292!kk       CALL location_message( TRIM( text ), .FALSE. )
     2293!        IF(myid == 0)  THEN
     2294!           write(9,*) ' '
     2295!           write(9,*) 'kpp setup '
     2296!           write(9,*) ' '
     2297!           write(9,*) '    gas_phase chemistry: solver_type = ',TRIM( solver_type )
     2298!           write(9,*) ' '
     2299!           write(9,*) '    Hstart  = ',rcntrl(3)
     2300!           write(9,*) '    FacMin  = ',rcntrl(4)
     2301!           write(9,*) '    FacMax  = ',rcntrl(5)
     2302!           write(9,*) ' '
     2303!           IF(vl_dim > 1)  THEN
     2304!              write(9,*) '    Vector mode                   vektor length = ',vl_dim
     2305!           ELSE
     2306!              write(9,*) '    Scalar mode'
     2307!           ENDIF
     2308!           write(9,*) ' '
     2309!        END IF
    23442310
    23452311    RETURN
     
    23472313 END SUBROUTINE chem_parin
    23482314
    2349  !
    2350  !------------------------------------------------------------------------------!
    2351  !
    2352  ! Description:
    2353  ! ------------
    2354  !> Subroutine calculating prognostic equations for chemical species
    2355  !> (vector-optimized).
    2356  !> Routine is called separately for each chemical species over a loop from
    2357  !> prognostic_equations.
    2358  !------------------------------------------------------------------------------!
     2315 
     2316!------------------------------------------------------------------------------!
     2317! Description:
     2318! ------------
     2319!> Subroutine calculating prognostic equations for chemical species
     2320!> (vector-optimized).
     2321!> Routine is called separately for each chemical species over a loop from
     2322!> prognostic_equations.
     2323!------------------------------------------------------------------------------!
    23592324 SUBROUTINE chem_prognostic_equations( cs_scalar_p, cs_scalar, tcs_scalar_m,   &
    23602325      pr_init_cs, ilsp )
     
    23902355
    23912356
    2392     !
    2393     !-- Tendency terms for chemical species
     2357!
     2358!-- Tendency terms for chemical species
    23942359    tend = 0.0_wp
    2395     !   
    2396     !-- Advection terms
     2360!   
     2361!-- Advection terms
    23972362    IF ( timestep_scheme(1:5) == 'runge' )  THEN
    23982363       IF ( ws_scheme_sca )  THEN
     
    24042369       CALL advec_s_up( cs_scalar )
    24052370    ENDIF
    2406     !
    2407     !-- Diffusion terms  (the last three arguments are zero)
     2371!
     2372!-- Diffusion terms  (the last three arguments are zero)
    24082373    CALL diffusion_s( cs_scalar,                                               &
    24092374         surf_def_h(0)%cssws(ilsp,:),                             &
     
    24242389         surf_usm_v(2)%cssws(ilsp,:),                             &
    24252390         surf_usm_v(3)%cssws(ilsp,:) )
    2426     !   
    2427     !-- Prognostic equation for chemical species
     2391!   
     2392!-- Prognostic equation for chemical species
    24282393    DO  i = nxl, nxr
    24292394       DO  j = nys, nyn   
     
    24432408       ENDDO
    24442409    ENDDO
    2445     !
    2446     !-- Calculate tendencies for the next Runge-Kutta step
     2410!
     2411!-- Calculate tendencies for the next Runge-Kutta step
    24472412    IF ( timestep_scheme(1:5) == 'runge' )  THEN
    24482413       IF ( intermediate_timestep_count == 1 )  THEN
     
    24692434 END SUBROUTINE chem_prognostic_equations
    24702435
    2471  !------------------------------------------------------------------------------!
    2472  !
    2473  ! Description:
    2474  ! ------------
    2475  !> Subroutine calculating prognostic equations for chemical species
    2476  !> (cache-optimized).
    2477  !> Routine is called separately for each chemical species over a loop from
    2478  !> prognostic_equations.
    2479  !------------------------------------------------------------------------------!
    2480  SUBROUTINE chem_prognostic_equations_ij( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs,     &
    2481       i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs,                  &
     2436
     2437!------------------------------------------------------------------------------!
     2438! Description:
     2439! ------------
     2440!> Subroutine calculating prognostic equations for chemical species
     2441!> (cache-optimized).
     2442!> Routine is called separately for each chemical species over a loop from
     2443!> prognostic_equations.
     2444!------------------------------------------------------------------------------!
     2445 SUBROUTINE chem_prognostic_equations_ij( cs_scalar_p, cs_scalar, tcs_scalar_m,                 &
     2446      pr_init_cs, i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs,                            &
    24822447      flux_l_cs, diss_l_cs )
    24832448    USE pegrid         
     
    25072472    REAL(wp), DIMENSION(0:nz+1)                                  :: pr_init_cs  !<
    25082473
    2509     !
    2510     !-- local variables
     2474!
     2475!-- local variables
    25112476
    25122477    INTEGER :: k
    2513     !
    2514     !--    Tendency-terms for chem spcs.
     2478!
     2479!--    Tendency-terms for chem spcs.
    25152480    tend(:,j,i) = 0.0_wp
    2516     !   
    2517     !-- Advection terms
     2481!   
     2482!-- Advection terms
    25182483    IF ( timestep_scheme(1:5) == 'runge' )  THEN
    25192484       IF ( ws_scheme_sca )  THEN
     
    25262491       CALL advec_s_up( i, j, cs_scalar )
    25272492    ENDIF
    2528 
    2529     !
    2530 
    2531     !-- Diffusion terms (the last three arguments are zero)
     2493!
     2494!-- Diffusion terms (the last three arguments are zero)
    25322495
    25332496    CALL diffusion_s( i, j, cs_scalar,                                                 &
     
    25412504         surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:),        &
    25422505         surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
    2543 
    2544     !   
    2545     !-- Prognostic equation for chem spcs
    2546     DO k = nzb+1, nzt
     2506!   
     2507!-- Prognostic equation for chem spcs
     2508    DO  k = nzb+1, nzt
    25472509       cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d  *                       &
    25482510            ( tsc(2) * tend(k,j,i) +         &
     
    25572519       IF ( cs_scalar_p(k,j,i) < 0.0_wp )  cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)    !FKS6
    25582520    ENDDO
    2559 
    2560     !
    2561     !-- Calculate tendencies for the next Runge-Kutta step
     2521!
     2522!-- Calculate tendencies for the next Runge-Kutta step
    25622523    IF ( timestep_scheme(1:5) == 'runge' )  THEN
    25632524       IF ( intermediate_timestep_count == 1 )  THEN
     
    25762537 END SUBROUTINE chem_prognostic_equations_ij
    25772538
    2578  !
    2579  !------------------------------------------------------------------------------!
    2580  !
    2581  ! Description:
    2582  ! ------------
    2583  !> Subroutine to read restart data of chemical species
    2584  !------------------------------------------------------------------------------!
    2585 
     2539
     2540!------------------------------------------------------------------------------!
     2541! Description:
     2542! ------------
     2543!> Subroutine to read restart data of chemical species
     2544!------------------------------------------------------------------------------!
    25862545 SUBROUTINE chem_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,             &
    25872546                            nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,   &
     
    26232582    IF ( ALLOCATED(chem_species) )  THEN
    26242583
    2625        DO lsp = 1, nspec
     2584       DO  lsp = 1, nspec
    26262585
    26272586          !< for time-averaged chemical conc.
    2628           spc_name_av  =  TRIM(chem_species(lsp)%name)//'_av'
    2629 
    2630           IF ( restart_string(1:length) == TRIM(chem_species(lsp)%name) )    &
    2631                THEN
     2587          spc_name_av  =  TRIM( chem_species(lsp)%name )//'_av'
     2588
     2589          IF ( restart_string(1:length) == TRIM( chem_species(lsp)%name) )    &
     2590             THEN
    26322591             !< read data into tmp_3d
    26332592             IF ( k == 1 )  READ ( 13 )  tmp_3d 
    26342593             !< fill ..%conc in the restart run   
    2635              chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,                  &
     2594             chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,                    &
    26362595                  nxlc-nbgp:nxrc+nbgp) =                  &
    26372596                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     
    26392598          ELSEIF (restart_string(1:length) == spc_name_av )  THEN
    26402599             IF ( k == 1 )  READ ( 13 )  tmp_3d
    2641              chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,               &
     2600             chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,                 &
    26422601                  nxlc-nbgp:nxrc+nbgp) =               &
    26432602                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     
    26512610
    26522611 END SUBROUTINE chem_rrd_local
    2653  !
    2654  !-------------------------------------------------------------------------------!
    2655  !> Description:
    2656  !> Calculation of horizontally averaged profiles
    2657  !> This routine is called for every statistic region (sr) defined by the user,
    2658  !> but at least for the region "total domain" (sr=0).
    2659  !> quantities.
    2660  !------------------------------------------------------------------------------!
     2612
     2613
     2614!-------------------------------------------------------------------------------!
     2615!> Description:
     2616!> Calculation of horizontally averaged profiles
     2617!> This routine is called for every statistic region (sr) defined by the user,
     2618!> but at least for the region "total domain" (sr=0).
     2619!> quantities.
     2620!-------------------------------------------------------------------------------!
    26612621 SUBROUTINE chem_statistics( mode, sr, tn )
    26622622
     
    26812641    IF ( mode == 'profiles' )  THEN
    26822642       !
    2683        !--    Sample on how to calculate horizontally averaged profiles of user-
    2684        !--    defined quantities. Each quantity is identified by the index
    2685        !--    "pr_palm+#" where "#" is an integer starting from 1. These
    2686        !--    user-profile-numbers must also be assigned to the respective strings
    2687        !--    given by data_output_pr_cs in routine user_check_data_output_pr.
    2688        !--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
    2689        !                       w*pt*), dim-4 = statistical region.
    2690 
    2691        !$OMP DO
     2643!
     2644!--    Sample on how to calculate horizontally averaged profiles of user-
     2645!--    defined quantities. Each quantity is identified by the index
     2646!--    "pr_palm+#" where "#" is an integer starting from 1. These
     2647!--    user-profile-numbers must also be assigned to the respective strings
     2648!--    given by data_output_pr_cs in routine user_check_data_output_pr.
     2649!--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
     2650!--                     w*pt*), dim-4 = statistical region.
     2651
     2652!$OMP DO
    26922653       DO  i = nxl, nxr
    26932654          DO  j = nys, nyn
     
    27102671 END SUBROUTINE chem_statistics
    27112672
    2712  !
    2713  !------------------------------------------------------------------------------!
    2714  !
    2715  ! Description:
    2716  ! ------------
    2717  !> Subroutine for swapping of timelevels for chemical species
    2718  !> called out from subroutine swap_timelevel
    2719  !------------------------------------------------------------------------------!
     2673
     2674!------------------------------------------------------------------------------!
     2675! Description:
     2676! ------------
     2677!> Subroutine for swapping of timelevels for chemical species
     2678!> called out from subroutine swap_timelevel
     2679!------------------------------------------------------------------------------!
     2680
    27202681
    27212682 SUBROUTINE chem_swap_timelevel( level )
     
    27242685
    27252686    INTEGER(iwp), INTENT(IN) ::  level
    2726     !
    2727     !-- local variables
     2687!
     2688!-- local variables
    27282689    INTEGER(iwp)             ::  lsp
    27292690
    27302691
    2731     IF ( level == 0 ) THEN
    2732        DO lsp=1, nvar                                       
     2692    IF ( level == 0 )  THEN
     2693       DO  lsp=1, nvar                                       
    27332694          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
    27342695          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
    27352696       ENDDO
    27362697    ELSE
    2737        DO lsp=1, nvar                                       
     2698       DO  lsp=1, nvar                                       
    27382699          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_2(:,:,:,lsp)
    27392700          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_1(:,:,:,lsp)
     
    27432704    RETURN
    27442705 END SUBROUTINE chem_swap_timelevel
    2745  !
    2746  !------------------------------------------------------------------------------!
    2747  !
    2748  ! Description:
    2749  ! ------------
    2750  !> Subroutine to write restart data for chemistry model
    2751  !------------------------------------------------------------------------------!
     2706
     2707
     2708!------------------------------------------------------------------------------!
     2709! Description:
     2710! ------------
     2711!> Subroutine to write restart data for chemistry model
     2712!------------------------------------------------------------------------------!
    27522713 SUBROUTINE chem_wrd_local
     2714
    27532715
    27542716    IMPLICIT NONE
     
    27662728
    27672729
    2768 
    2769  !-------------------------------------------------------------------------------!
    2770  !
    2771  ! Description:
    2772  ! ------------
    2773  !> Subroutine to calculate the deposition of gases and PMs. For now deposition
    2774  !> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
    2775  !> considered. The deposition of particles is derived following Zhang et al.,
    2776  !> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
    2777  !>     
    2778  !> @TODO: Consider deposition on vertical surfaces   
    2779  !> @TODO: Consider overlaying horizontal surfaces
    2780  !> @TODO: Consider resolved vegetation   
    2781  !> @TODO: Check error messages
    2782  !-------------------------------------------------------------------------------!
    2783 
     2730!-------------------------------------------------------------------------------!
     2731! Description:
     2732! ------------
     2733!> Subroutine to calculate the deposition of gases and PMs. For now deposition
     2734!> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
     2735!> considered. The deposition of particles is derived following Zhang et al.,
     2736!> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
     2737!>     
     2738!> @TODO: Consider deposition on vertical surfaces   
     2739!> @TODO: Consider overlaying horizontal surfaces
     2740!> @TODO: Consider resolved vegetation   
     2741!> @TODO: Check error messages
     2742!-------------------------------------------------------------------------------!
    27842743 SUBROUTINE chem_depo( i, j )
    27852744
     
    28012760
    28022761
    2803 
    28042762    IMPLICIT NONE
    28052763
    28062764    INTEGER(iwp), INTENT(IN) ::  i
    28072765    INTEGER(iwp), INTENT(IN) ::  j
    2808 
    28092766    INTEGER(iwp) ::  k                             !< matching k to surface m at i,j
    28102767    INTEGER(iwp) ::  lsp                           !< running index for chem spcs.
    2811 !    INTEGER(iwp) ::  lu                            !< running index for landuse classes
     2768!    INTEGER(iwp) ::  lu                           !< running index for landuse classes
    28122769    INTEGER(iwp) ::  lu_palm                       !< index of PALM LSM vegetation_type at current surface element
    28132770    INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current surface element
     
    28262783    INTEGER(iwp) ::  pspec                         !< running index
    28272784    INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
    2828 
    2829     !
    2830     !-- Vegetation                                              !<Match to DEPAC
     2785!
     2786!-- Vegetation                                              !<Match to DEPAC
    28312787    INTEGER(iwp) ::  ind_lu_user = 0                         !< --> ERROR 
    28322788    INTEGER(iwp) ::  ind_lu_b_soil = 1                       !< --> ilu_desert
     
    28482804    INTEGER(iwp) ::  ind_lu_mixed_forest = 17                !< --> ilu_coniferous_forest (ave(decid+conif))
    28492805    INTEGER(iwp) ::  ind_lu_intrup_forest = 18               !< --> ilu_other (ave(other+decid))
    2850 
    2851     !
    2852     !-- Water
     2806!
     2807!-- Water
    28532808    INTEGER(iwp) ::  ind_luw_user = 0                        !< --> ERROR
    28542809    INTEGER(iwp) ::  ind_luw_lake = 1                        !< --> ilu_water_inland
     
    28572812    INTEGER(iwp) ::  ind_luw_pond = 4                        !< --> ilu_water_inland
    28582813    INTEGER(iwp) ::  ind_luw_fountain = 5                    !< --> ilu_water_inland
    2859 
    2860     !
    2861     !-- Pavement
     2814!
     2815!-- Pavement
    28622816    INTEGER(iwp) ::  ind_lup_user = 0                        !< --> ERROR
    28632817    INTEGER(iwp) ::  ind_lup_asph_conc = 1                   !< --> ilu_desert
     
    28762830    INTEGER(iwp) ::  ind_lup_art_turf = 14                   !< --> ilu_desert
    28772831    INTEGER(iwp) ::  ind_lup_clay = 15                       !< --> ilu_desert
    2878 
    2879 
    2880     !
    2881     !-- Particle parameters according to the respective aerosol classes (PM25, PM10)
     2832!
     2833!-- Particle parameters according to the respective aerosol classes (PM25, PM10)
    28822834    INTEGER(iwp) ::  ind_p_size = 1     !< index for partsize in particle_pars
    28832835    INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
     
    29292881    REAL(wp) ::  ts             !< surface temperatur in degrees celsius
    29302882    REAL(wp) ::  qv_tmp         !< surface mixing ratio at current surface element
    2931 
    2932     !
    2933     !-- Particle parameters (PM10 (1), PM25 (2))
    2934     !-- partsize (diameter in m), rhopart (density in kg/m3), slipcor
    2935     !-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
     2883!
     2884!-- Particle parameters (PM10 (1), PM25 (2))
     2885!-- partsize (diameter in m), rhopart (density in kg/m3), slipcor
     2886!-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
    29362887    REAL(wp), DIMENSION(1:3,1:2), PARAMETER ::  particle_pars = RESHAPE( (/ &
    29372888         8.0e-6_wp, 1.14e3_wp, 1.016_wp, &  !<  1
     
    29392890         /), (/ 3, 2 /) )
    29402891
    2941 
    29422892    LOGICAL ::  match_lsm     !< flag indicating natural-type surface
    29432893    LOGICAL ::  match_usm     !< flag indicating urban-type surface
    2944 
    2945 
    2946     !
    2947     !-- List of names of possible tracers
    2948     CHARACTER(len=*), PARAMETER ::  pspecnames(nposp) = (/ &
     2894!
     2895!-- List of names of possible tracers
     2896    CHARACTER(LEN=*), PARAMETER ::  pspecnames(nposp) = (/ &
    29492897         'NO2           ', &    !< NO2
    29502898         'NO            ', &    !< NO
     
    30162964         'tpm25_biascorr', &    !< tpm25_biascorr
    30172965         'O3_biascorr   ' /)    !< o3_biascorr
    3018 
    3019 
    3020     !
    3021     !-- tracer mole mass:
     2966!
     2967!-- tracer mole mass:
    30222968    REAL(wp), PARAMETER ::  specmolm(nposp) = (/ &
    30232969         xm_O * 2 + xm_N, &                         !< NO2
     
    30903036         xm_dummy, &                                !< tpm25_biascorr
    30913037         xm_O * 3 /)                                !< o3_biascorr
    3092 
    3093 
    3094     !
    3095     !-- Initialize surface element m
     3038!
     3039!-- Initialize surface element m
    30963040    m = 0
    30973041    k = 0
    3098     !
    3099     !-- LSM or USM surface present at i,j:
    3100     !-- Default surfaces are NOT considered for deposition
     3042!
     3043!-- LSM or USM surface present at i,j:
     3044!-- Default surfaces are NOT considered for deposition
    31013045    match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
    31023046    match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
    3103 
    3104 
    3105     !
    3106     !--For LSM surfaces
     3047!
     3048!--For LSM surfaces
    31073049
    31083050    IF ( match_lsm )  THEN
    3109 
    3110        !
    3111        !-- Get surface element information at i,j:
     3051!
     3052!--    Get surface element information at i,j:
    31123053       m = surf_lsm_h%start_index(j,i)
    31133054       k = surf_lsm_h%k(m)
    3114 
    3115        !
    3116        !-- Get needed variables for surface element m
     3055!
     3056!--    Get needed variables for surface element m
    31173057       ustar_lu  = surf_lsm_h%us(m)
    31183058       z0h_lu    = surf_lsm_h%z0h(m)
     
    31213061       lai = surf_lsm_h%lai(m)
    31223062       sai = lai + 1
    3123 
    3124        !
    3125        !-- For small grid spacing neglect R_a
     3063!
     3064!--    For small grid spacing neglect R_a
    31263065       IF ( dzw(k) <= 1.0 )  THEN
    31273066          r_aero_lu = 0.0_wp
    31283067       ENDIF
    3129 
    3130        !
    3131        !--Initialize lu's
     3068!
     3069!--    Initialize lu's
    31323070       lu_palm = 0
    31333071       lu_dep = 0
     
    31363074       luw_palm = 0
    31373075       luw_dep = 0
    3138 
    3139        !
    3140        !--Initialize budgets
     3076!
     3077!--    Initialize budgets
    31413078       bud_now_lu  = 0.0_wp
    31423079       bud_now_lup = 0.0_wp
    31433080       bud_now_luw = 0.0_wp
    3144 
    3145 
    3146        !
    3147        !-- Get land use for i,j and assign to DEPAC lu
     3081!
     3082!--    Get land use for i,j and assign to DEPAC lu
    31483083       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
    31493084          lu_palm = surf_lsm_h%vegetation_type(m)
     
    32453180          ENDIF
    32463181       ENDIF
    3247 
    3248 
    3249        !
    3250        !-- Set wetness indicator to dry or wet for lsm vegetation or pavement
     3182!
     3183!--    Set wetness indicator to dry or wet for lsm vegetation or pavement
    32513184       IF ( surf_lsm_h%c_liq(m) > 0 )  THEN
    32523185          nwet = 1
     
    32543187          nwet = 0
    32553188       ENDIF
    3256 
    3257        !
    3258        !--Compute length of time step
     3189!
     3190!--    Compute length of time step
    32593191       IF ( call_chem_at_all_substeps )  THEN
    32603192          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
     
    32633195       ENDIF
    32643196
    3265 
    32663197       dh = dzw(k)
    32673198       inv_dh = 1.0_wp / dh
    32683199       dt_dh = dt_chem/dh
    3269 
    3270        !
    3271        !-- Concentration at i,j,k
     3200!
     3201!--    Concentration at i,j,k
    32723202       DO  lsp = 1, nspec
    32733203          cc(lsp) = chem_species(lsp)%conc(k,j,i)
    32743204       ENDDO
    32753205
    3276 
    3277        !
    3278        !-- Temperature at i,j,k
     3206!--    Temperature at i,j,k
    32793207       ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
    32803208
    32813209       ts       = ttemp - 273.15  !< in degrees celcius
    3282 
    3283        !
    3284        !-- Viscosity of air
     3210!
     3211!--    Viscosity of air
    32853212       visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
    3286 
    3287        !
    3288        !-- Air density at k
     3213!
     3214!--    Air density at k
    32893215       dens = rho_air_zw(k)
    3290 
    3291        !
    3292        !-- Calculate relative humidity from specific humidity for DEPAC
     3216!
     3217!--    Calculate relative humidity from specific humidity for DEPAC
    32933218       qv_tmp = MAX(q(k,j,i),0.0_wp)
    32943219       relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(k) )
    3295 
    3296 
    3297 
    3298        !
    3299        !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
    3300        !-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
    3301 
    3302        !
    3303        !-- Vegetation
     3220!
     3221!-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
     3222!-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
     3223!
     3224!--    Vegetation
    33043225       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
    33053226
    3306 
    33073227          slinnfac = 1.0_wp
    3308 
    3309           !
    3310           !-- Get deposition velocity vd
     3228!
     3229!--       Get deposition velocity vd
    33113230          DO  lsp = 1, nvar
    3312              !
    3313              !-- Initialize
     3231!
     3232!--          Initialize
    33143233             vs = 0.0_wp
    33153234             vd_lu = 0.0_wp
     
    33193238             IF ( spc_names(lsp) == 'PM10' )  THEN
    33203239                part_type = 1
    3321                 !
    3322                 !-- Sedimentation velocity
     3240!
     3241!--            Sedimentation velocity
    33233242                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    33243243                     particle_pars(ind_p_size, part_type), &
     
    33403259             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
    33413260                part_type = 2
    3342                 !
    3343                 !-- Sedimentation velocity
     3261!
     3262!--            Sedimentation velocity
    33443263                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    33453264                     particle_pars(ind_p_size, part_type), &
    33463265                     particle_pars(ind_p_slip, part_type), &
    33473266                     visc)
    3348 
    33493267
    33503268                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     
    33593277                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    33603278
    3361 
    33623279             ELSE !< GASES
    3363 
    3364                 !
    3365                 !-- Read spc_name of current species for gas parameter
    3366 
     3280!
     3281!--             Read spc_name of current species for gas parameter
    33673282                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
    33683283                   i_pspec = 0
     
    33743289
    33753290                ELSE
    3376                    !
    3377                    !-- For now species not deposited
     3291!
     3292!--            For now species not deposited
    33783293                   CYCLE
    33793294                ENDIF
    3380 
    3381                 !
    3382                 !-- Factor used for conversion from ppb to ug/m3 :
    3383                 !-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
    3384                 !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
    3385                 !   c           1e-9              xm_tracer         1e9       /       xm_air            dens
    3386                 !-- thus:
    3387                 !   c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
    3388                 !-- Use density at k:
     3295!
     3296!--             Factor used for conversion from ppb to ug/m3 :
     3297!--             ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3298!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3299!--                 c           1e-9              xm_tracer         1e9       /       xm_air            dens
     3300!--             thus:
     3301!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3302!--             Use density at k:
    33893303
    33903304                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    3391 
    3392                 !
    3393                 !-- Atmospheric concentration in DEPAC is requested in ug/m3:
     3305!
     3306!--             Atmospheric concentration in DEPAC is requested in ug/m3:
    33943307                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    33953308                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
    3396 
    3397                 !
    3398                 !-- Diffusivity for DEPAC relevant gases
    3399                 !-- Use default value
     3309!
     3310!--             Diffusivity for DEPAC relevant gases
     3311!--             Use default value
    34003312                diffc            = 0.11e-4
    3401                 !
    3402                 !-- overwrite with known coefficients of diffusivity from Massman (1998)
     3313!
     3314!--            overwrite with known coefficients of diffusivity from Massman (1998)
    34033315                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    34043316                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     
    34083320                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    34093321                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    3410 
    3411 
    3412                 !
    3413                 !-- Get quasi-laminar boundary layer resistance Rb:
     3322!
     3323!--             Get quasi-laminar boundary layer resistance Rb:
    34143324                CALL get_rb_cell( (lu_dep == ilu_water_sea) .OR. (lu_dep == ilu_water_inland), &
    34153325                     z0h_lu, ustar_lu, diffc, &
    34163326                     Rb )
    3417 
    3418                 !
    3419                 !-- Get Rc_tot
     3327!
     3328!--             Get Rc_tot
    34203329                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    34213330                     relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    34223331                     r_aero_lu , Rb )
    3423 
    3424 
    3425                 !
    3426                 !-- Calculate budget
     3332!
     3333!--             Calculate budget
    34273334                IF ( Rc_tot <= 0.0 )  THEN
    34283335
     
    34403347          ENDDO
    34413348       ENDIF
    3442 
    3443 
    3444        !
    3445        !-- Pavement
     3349!
     3350!--    Pavement
    34463351       IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
    3447 
    3448 
    3449           !
    3450           !-- No vegetation on pavements:
     3352!
     3353!--       No vegetation on pavements:
    34513354          lai = 0.0_wp
    34523355          sai = 0.0_wp
    34533356
    34543357          slinnfac = 1.0_wp
    3455 
    3456           !
    3457           !-- Get vd
     3358!
     3359!--       Get vd
    34583360          DO  lsp = 1, nvar
    3459              !
    3460              !-- Initialize
     3361!
     3362!--      Initialize
    34613363             vs = 0.0_wp
    34623364             vd_lu = 0.0_wp
     
    34663368             IF ( spc_names(lsp) == 'PM10' )  THEN
    34673369                part_type = 1
    3468                 !
    3469                 !-- Sedimentation velocity:
     3370!
     3371!--            Sedimentation velocity:
    34703372                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    34713373                     particle_pars(ind_p_size, part_type), &
     
    34873389             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
    34883390                part_type = 2
    3489                 !
    3490                 !-- Sedimentation velocity:
     3391!
     3392!--            Sedimentation velocity:
    34913393                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    34923394                     particle_pars(ind_p_size, part_type), &
    34933395                     particle_pars(ind_p_slip, part_type), &
    34943396                     visc)
    3495 
    34963397
    34973398                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     
    35063407                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    35073408
    3508 
    35093409             ELSE  !<GASES
    3510 
    3511                 !
    3512                 !-- Read spc_name of current species for gas parameter
     3410!
     3411!--             Read spc_name of current species for gas parameter
    35133412
    35143413                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
     
    35213420
    35223421                ELSE
    3523                    !
    3524                    !-- For now species not deposited
     3422!
     3423!--                For now species not deposited
    35253424                   CYCLE
    35263425                ENDIF
    3527 
    3528                 !-- Factor used for conversion from ppb to ug/m3 :
    3529                 !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
    3530                 !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
    3531                 !   c           1e-9               xm_tracer         1e9       /       xm_air            dens
    3532                 !-- thus:
    3533                 !   c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
    3534                 !-- Use density at lowest layer:
     3426!
     3427!--            Factor used for conversion from ppb to ug/m3 :
     3428!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3429!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3430!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
     3431!--            thus:
     3432!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3433!--            Use density at lowest layer:
    35353434
    35363435                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    3537 
    3538                 !
    3539                 !-- Atmospheric concentration in DEPAC is requested in ug/m3:
     3436!
     3437!--             Atmospheric concentration in DEPAC is requested in ug/m3:
    35403438                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    35413439                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
    3542 
    3543                 !
    3544                 !-- Diffusivity for DEPAC relevant gases
    3545                 !-- Use default value
     3440!
     3441!--             Diffusivity for DEPAC relevant gases
     3442!--             Use default value
    35463443                diffc            = 0.11e-4
    3547                 !
    3548                 !-- overwrite with known coefficients of diffusivity from Massman (1998)
     3444!
     3445!--            overwrite with known coefficients of diffusivity from Massman (1998)
    35493446                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    35503447                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     
    35543451                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    35553452                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    3556 
    3557 
    3558                 !
    3559                 !-- Get quasi-laminar boundary layer resistance Rb:
    3560                 CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland), &
    3561                      z0h_lu, ustar_lu, diffc, &
    3562                      Rb )
    3563 
    3564 
    3565                 !
    3566                 !-- Get Rc_tot
    3567                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    3568                      relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    3569                      r_aero_lu , Rb )
    3570 
    3571 
    3572                 !
    3573                 !-- Calculate budget
     3453!
     3454!--             Get quasi-laminar boundary layer resistance Rb:
     3455                CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland),   &
     3456                     z0h_lu, ustar_lu, diffc, Rb )
     3457!
     3458!--             Get Rc_tot
     3459                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,      &
     3460                                         glrad, cos_zenith, relh, lai, sai, nwet, lup_dep, 2,      &
     3461                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,             &
     3462                                         r_aero_lu , Rb )
     3463!
     3464!--             Calculate budget
    35743465                IF ( Rc_tot <= 0.0 )  THEN
    3575 
    35763466                   bud_now_lup(lsp) = 0.0_wp
    3577 
    35783467                ELSE
    3579 
    35803468                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
    3581 
    35823469                   bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
    35833470                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    35843471                ENDIF
    35853472
    3586 
    35873473             ENDIF
    35883474          ENDDO
    35893475       ENDIF
    3590 
    3591 
    3592        !
    3593        !-- Water
     3476!
     3477!--    Water
    35943478       IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
    3595 
    3596 
    3597           !
    3598           !-- No vegetation on water:
     3479!
     3480!--       No vegetation on water:
    35993481          lai = 0.0_wp
    36003482          sai = 0.0_wp
    3601 
    36023483          slinnfac = 1.0_wp
    3603 
    3604           !
    3605           !-- Get vd
     3484!
     3485!--       Get vd
    36063486          DO  lsp = 1, nvar
    3607              !
    3608              !-- Initialize
     3487!
     3488!--          Initialize
    36093489             vs = 0.0_wp
    36103490             vd_lu = 0.0_wp
     
    36143494             IF ( spc_names(lsp) == 'PM10' )  THEN
    36153495                part_type = 1
    3616                 !
    3617                 !-- Sedimentation velocity:
     3496!
     3497!--            Sedimentation velocity:
    36183498                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    36193499                     particle_pars(ind_p_size, part_type), &
     
    36323512                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    36333513
    3634 
    36353514             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
    36363515                part_type = 2
    3637                 !
    3638                 !-- Sedimentation velocity:
     3516!
     3517!--            Sedimentation velocity:
    36393518                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    36403519                     particle_pars(ind_p_size, part_type), &
    36413520                     particle_pars(ind_p_slip, part_type), &
    36423521                     visc)
    3643 
    36443522
    36453523                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     
    36543532                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    36553533
    3656 
    36573534             ELSE  !<GASES
    3658 
    3659                 !
    3660                 !-- Read spc_name of current species for gas PARAMETER
     3535!
     3536!--             Read spc_name of current species for gas PARAMETER
    36613537
    36623538                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
     
    36693545
    36703546                ELSE
    3671                    !
    3672                    !-- For now species not deposited
     3547!
     3548!--                For now species not deposited
    36733549                   CYCLE
    36743550                ENDIF
    3675 
    3676                 !-- Factor used for conversion from ppb to ug/m3 :
    3677                 !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
    3678                 !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
    3679                 !   c           1e-9               xm_tracer         1e9       /       xm_air            dens
    3680                 !-- thus:
    3681                 !   c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
    3682                 !-- Use density at lowest layer:
    3683 
    3684                 ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    3685 
    3686                 !
    3687                 !-- Atmospheric concentration in DEPAC is requested in ug/m3:
    3688                 !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    3689                 catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
    3690 
    3691                 !
    3692                 !-- Diffusivity for DEPAC relevant gases
    3693                 !-- Use default value
     3551!
     3552!--             Factor used for conversion from ppb to ug/m3 :
     3553!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3554!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3555!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
     3556!--             thus:
     3557!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3558!--             Use density at lowest layer:
     3559
     3560                ppm_to_ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
     3561!
     3562!--             Atmospheric concentration in DEPAC is requested in ug/m3:
     3563!--                 ug/m3        ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     3564                catm = cc(lsp) * ppm_to_ugm3 *  specmolm(i_pspec)  ! in ug/m3
     3565!
     3566!--             Diffusivity for DEPAC relevant gases
     3567!--             Use default value
    36943568                diffc            = 0.11e-4
    3695                 !
    3696                 !-- overwrite with known coefficients of diffusivity from Massman (1998)
     3569!
     3570!--            overwrite with known coefficients of diffusivity from Massman (1998)
    36973571                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    36983572                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     
    37023576                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    37033577                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    3704 
    3705 
    3706                 !
    3707                 !-- Get quasi-laminar boundary layer resistance Rb:
    3708                 CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland), &
    3709                      z0h_lu, ustar_lu, diffc, &
    3710                      Rb )
    3711 
    3712                 !
    3713                 !-- Get Rc_tot
    3714                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    3715                      relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    3716                      r_aero_lu , Rb )
    3717 
    3718 
    3719                 ! Calculate budget
     3578!
     3579!--             Get quasi-laminar boundary layer resistance Rb:
     3580                CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland),  &
     3581                     z0h_lu, ustar_lu, diffc, rb )
     3582
     3583!--             Get Rc_tot
     3584                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,      &
     3585                                         glrad, cos_zenith, relh, lai, sai, nwet, luw_dep, 2,      &
     3586                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,            &
     3587                                          r_aero_lu , Rb )
     3588!
     3589!--             Calculate budget
    37203590                IF ( Rc_tot <= 0.0 )  THEN
    37213591
     
    37363606
    37373607       bud_now = 0.0_wp
    3738        !
    3739        !-- Calculate overall budget for surface m and adapt concentration
     3608!
     3609!--    Calculate overall budget for surface m and adapt concentration
    37403610       DO  lsp = 1, nspec
    3741 
    37423611
    37433612          bud_now(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_now_lu(lsp) + &
    37443613               surf_lsm_h%frac(ind_pav_green,m) * bud_now_lup(lsp) + &
    37453614               surf_lsm_h%frac(ind_wat_win,m) * bud_now_luw(lsp)
    3746 
    3747           !
    3748           !-- Compute new concentration:
     3615!
     3616!--       Compute new concentration:
    37493617          cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh
    37503618
     
    37543622
    37553623    ENDIF
    3756 
    3757 
    3758 
    3759 
    3760     !
    3761     !-- For USM surfaces   
     3624!
     3625!-- For USM surfaces   
    37623626
    37633627    IF ( match_usm )  THEN
    3764 
    3765        !
    3766        !-- Get surface element information at i,j:
     3628!
     3629!--    Get surface element information at i,j:
    37673630       m = surf_usm_h%start_index(j,i)
    37683631       k = surf_usm_h%k(m)
    3769 
    3770        !
    3771        !-- Get needed variables for surface element m
     3632!
     3633!--    Get needed variables for surface element m
    37723634       ustar_lu  = surf_usm_h%us(m)
    37733635       z0h_lu    = surf_usm_h%z0h(m)
     
    37763638       lai = surf_usm_h%lai(m)
    37773639       sai = lai + 1
    3778 
    3779        !
    3780        !-- For small grid spacing neglect R_a
     3640!
     3641!--    For small grid spacing neglect R_a
    37813642       IF ( dzw(k) <= 1.0 )  THEN
    37823643          r_aero_lu = 0.0_wp
    37833644       ENDIF
    3784 
    3785        !
    3786        !-- Initialize lu's
     3645!
     3646!--    Initialize lu's
    37873647       luu_palm = 0
    37883648       luu_dep = 0
     
    37913651       lud_palm = 0
    37923652       lud_dep = 0
    3793 
    3794        !
    3795        !-- Initialize budgets
     3653!
     3654!--    Initialize budgets
    37963655       bud_now_luu  = 0.0_wp
    37973656       bud_now_lug = 0.0_wp
    37983657       bud_now_lud = 0.0_wp
    3799 
    3800 
    3801        !
    3802        !-- Get land use for i,j and assign to DEPAC lu
     3658!
     3659!--    Get land use for i,j and assign to DEPAC lu
    38033660       IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
    3804           !
    3805           !-- For green urban surfaces (e.g. green roofs
    3806           !-- assume LU short grass
     3661!
     3662!--      For green urban surfaces (e.g. green roofs
     3663!--      assume LU short grass
    38073664          lug_palm = ind_lu_s_grass
    38083665          IF ( lug_palm == ind_lu_user )  THEN
     
    38493706
    38503707       IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 )  THEN
    3851           !
    3852           !-- For walls in USM assume concrete walls/roofs,
    3853           !-- assumed LU class desert as also assumed for
    3854           !-- pavements in LSM
     3708!
     3709!--      For walls in USM assume concrete walls/roofs,
     3710!--      assumed LU class desert as also assumed for
     3711!--      pavements in LSM
    38553712          luu_palm = ind_lup_conc
    38563713          IF ( luu_palm == ind_lup_user )  THEN
     
    38913748
    38923749       IF ( surf_usm_h%frac(ind_wat_win,m) > 0 )  THEN
    3893           !
    3894           !-- For windows in USM assume metal as this is
    3895           !-- as close as we get, assumed LU class desert
    3896           !-- as also assumed for pavements in LSM
     3750!
     3751!--      For windows in USM assume metal as this is
     3752!--      as close as we get, assumed LU class desert
     3753!--      as also assumed for pavements in LSM
    38973754          lud_palm = ind_lup_metal     
    38983755          IF ( lud_palm == ind_lup_user )  THEN
     
    39313788          ENDIF
    39323789       ENDIF
    3933 
    3934 
    3935        !
    3936        !-- @TODO: Activate these lines as soon as new ebsolver branch is merged:
    3937        !-- Set wetness indicator to dry or wet for usm vegetation or pavement
     3790!
     3791!--    @TODO: Activate these lines as soon as new ebsolver branch is merged:
     3792!--    Set wetness indicator to dry or wet for usm vegetation or pavement
    39383793       !IF ( surf_usm_h%c_liq(m) > 0 )  THEN
    39393794       !   nwet = 1
     
    39413796       nwet = 0
    39423797       !ENDIF
    3943 
    3944        !
    3945        !-- Compute length of time step
     3798!
     3799!--    Compute length of time step
    39463800       IF ( call_chem_at_all_substeps )  THEN
    39473801          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
     
    39503804       ENDIF
    39513805
    3952 
    39533806       dh = dzw(k)
    39543807       inv_dh = 1.0_wp / dh
    39553808       dt_dh = dt_chem/dh
    3956 
    3957        !
    3958        !-- Concentration at i,j,k
     3809!
     3810!--    Concentration at i,j,k
    39593811       DO  lsp = 1, nspec
    39603812          cc(lsp) = chem_species(lsp)%conc(k,j,i)
    39613813       ENDDO
    3962 
    3963        !
    3964        !-- Temperature at i,j,k
     3814!
     3815!--    Temperature at i,j,k
    39653816       ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
    39663817
    39673818       ts       = ttemp - 273.15  !< in degrees celcius
    3968 
    3969        !
    3970        !-- Viscosity of air
     3819!
     3820!--    Viscosity of air
    39713821       visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
    3972 
    3973        !
    3974        !-- Air density at k
     3822!
     3823!--    Air density at k
    39753824       dens = rho_air_zw(k)
    3976 
    3977        !
    3978        !-- Calculate relative humidity from specific humidity for DEPAC
     3825!
     3826!--    Calculate relative humidity from specific humidity for DEPAC
    39793827       qv_tmp = MAX(q(k,j,i),0.0_wp)
    39803828       relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(nzb) )
    3981 
    3982 
    3983 
    3984        !
    3985        !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
    3986        !-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
    3987 
    3988        !
    3989        !-- Walls/roofs
     3829!
     3830!--    Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
     3831!--    for each surface fraction. Then derive overall budget taking into account the surface fractions.
     3832!
     3833!--    Walls/roofs
    39903834       IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 )  THEN
    3991 
    3992 
    3993           !
    3994           !-- No vegetation on non-green walls:
     3835!
     3836!--       No vegetation on non-green walls:
    39953837          lai = 0.0_wp
    39963838          sai = 0.0_wp
    39973839
    39983840          slinnfac = 1.0_wp
    3999 
    4000           !
    4001           !-- Get vd
     3841!
     3842!--       Get vd
    40023843          DO  lsp = 1, nvar
    4003              !
    4004              !-- Initialize
     3844!
     3845!--          Initialize
    40053846             vs = 0.0_wp
    40063847             vd_lu = 0.0_wp
     
    40103851             IF (spc_names(lsp) == 'PM10' )  THEN
    40113852                part_type = 1
    4012                 !
    4013                 !-- Sedimentation velocity
     3853!
     3854!--            Sedimentation velocity
    40143855                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    40153856                     particle_pars(ind_p_size, part_type), &
     
    40283869                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    40293870
    4030 
    40313871             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
    40323872                part_type = 2
    4033                 !
    4034                 !-- Sedimentation velocity
     3873!
     3874!--            Sedimentation velocity
    40353875                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    40363876                     particle_pars(ind_p_size, part_type), &
    40373877                     particle_pars(ind_p_slip, part_type), &
    40383878                     visc)
    4039 
    40403879
    40413880                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     
    40503889                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    40513890
    4052 
    40533891             ELSE  !< GASES
    4054 
    4055                 !
    4056                 !-- Read spc_name of current species for gas parameter
     3892!
     3893!--             Read spc_name of current species for gas parameter
    40573894
    40583895                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
     
    40643901                   ENDDO
    40653902                ELSE
    4066                    !
    4067                    !-- For now species not deposited
     3903!
     3904!--                For now species not deposited
    40683905                   CYCLE
    40693906                ENDIF
     3907!
     3908!--             Factor used for conversion from ppb to ug/m3 :
     3909!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3910!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3911!--                 c           1e-9              xm_tracer         1e9       /       xm_air            dens
     3912!--             thus:
     3913!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3914!--             Use density at k:
     3915
     3916                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    40703917
    40713918                !
    4072                 !-- Factor used for conversion from ppb to ug/m3 :
    4073                 !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
    4074                 !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
    4075                 !   c           1e-9              xm_tracer         1e9       /       xm_air            dens
    4076                 !-- thus:
    4077                 !   c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
    4078                 !-- Use density at k:
    4079 
    4080                 ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    4081 
    4082                 !
    4083                 !-- Atmospheric concentration in DEPAC is requested in ug/m3:
    4084                 !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    4085                 catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
    4086 
    4087                 !
    4088                 !-- Diffusivity for DEPAC relevant gases
    4089                 !-- Use default value
     3919!--             Atmospheric concentration in DEPAC is requested in ug/m3:
     3920!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     3921!--             catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3922!
     3923!--             Diffusivity for DEPAC relevant gases
     3924!--             Use default value
    40903925                diffc            = 0.11e-4
    4091                 !
    4092                 !-- overwrite with known coefficients of diffusivity from Massman (1998)
     3926!
     3927!--            overwrite with known coefficients of diffusivity from Massman (1998)
    40933928                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    40943929                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     
    40983933                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    40993934                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    4100 
    4101 
    4102                 !
    4103                 !-- Get quasi-laminar boundary layer resistance Rb:
    4104                 CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland), &
     3935!
     3936!--             Get quasi-laminar boundary layer resistance Rb:
     3937                CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland),   &
    41053938                     z0h_lu, ustar_lu, diffc, &
    41063939                     Rb )
    4107 
    4108                 !
    4109                 !-- Get Rc_tot
    4110                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    4111                      relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    4112                      r_aero_lu , Rb )
    4113 
    4114 
    4115                 !
    4116                 !-- Calculate budget
     3940!
     3941!--             Get Rc_tot
     3942                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,       &
     3943                                         glrad, cos_zenith, relh, lai, sai, nwet, luu_dep, 2,       &
     3944                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,             &
     3945                                         r_aero_lu, rb )
     3946!
     3947!--             Calculate budget
    41173948                IF ( Rc_tot <= 0.0 )  THEN
    41183949
     
    41303961          ENDDO
    41313962       ENDIF
    4132 
    4133 
    4134        !
    4135        !-- Green usm surfaces
     3963!
     3964!--    Green usm surfaces
    41363965       IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
    41373966
    4138 
    41393967          slinnfac = 1.0_wp
    4140 
    4141           !
    4142           !-- Get vd
     3968!
     3969!--       Get vd
    41433970          DO  lsp = 1, nvar
    4144              !
    4145              !-- Initialize
     3971!
     3972!--          Initialize
    41463973             vs = 0.0_wp
    41473974             vd_lu = 0.0_wp
     
    41513978             IF ( spc_names(lsp) == 'PM10' )  THEN
    41523979                part_type = 1
    4153                 ! Sedimentation velocity
     3980!
     3981!--             Sedimentation velocity
    41543982                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    41553983                     particle_pars(ind_p_size, part_type), &
     
    41683996                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    41693997
    4170 
    41713998             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
    41723999                part_type = 2
    4173                 ! Sedimentation velocity
     4000!
     4001!--             Sedimentation velocity
    41744002                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    41754003                     particle_pars(ind_p_size, part_type), &
    41764004                     particle_pars(ind_p_slip, part_type), &
    41774005                     visc)
    4178 
    41794006
    41804007                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     
    41894016                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    41904017
    4191 
    41924018             ELSE  !< GASES
    4193 
    4194                 !
    4195                 !-- Read spc_name of current species for gas parameter
     4019!
     4020!--             Read spc_name of current species for gas parameter
    41964021
    41974022                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
     
    42034028                   ENDDO
    42044029                ELSE
    4205                    !
    4206                    !-- For now species not deposited
     4030!
     4031!--                For now species not deposited
    42074032                   CYCLE
    42084033                ENDIF
    4209 
    4210                 !
    4211                 !-- Factor used for conversion from ppb to ug/m3 :
    4212                 !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
    4213                 !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
    4214                 !   c           1e-9               xm_tracer         1e9       /       xm_air            dens
    4215                 !-- thus:
    4216                 !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
    4217                 !-- Use density at k:
     4034!
     4035!--             Factor used for conversion from ppb to ug/m3 :
     4036!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     4037!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     4038!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
     4039!--             thus:
     4040!--                  c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     4041!--             Use density at k:
    42184042
    42194043                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
    4220 
    4221                 !
    4222                 !-- Atmospheric concentration in DEPAC is requested in ug/m3:
     4044!
     4045!--             Atmospheric concentration in DEPAC is requested in ug/m3:
    42234046                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    42244047                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
    4225 
    4226                 !
    4227                 !-- Diffusivity for DEPAC relevant gases
    4228                 !-- Use default value
     4048!
     4049!--             Diffusivity for DEPAC relevant gases
     4050!--             Use default value
    42294051                diffc            = 0.11e-4
    4230                 !
    4231                 !-- overwrite with known coefficients of diffusivity from Massman (1998)
     4052!
     4053!--            overwrite with known coefficients of diffusivity from Massman (1998)
    42324054                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    42334055                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     
    42374059                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    42384060                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    4239 
    4240 
    4241                 !
    4242                 !-- Get quasi-laminar boundary layer resistance Rb:
    4243                 CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland), &
     4061!
     4062!--             Get quasi-laminar boundary layer resistance Rb:
     4063                CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland),    &
    42444064                     z0h_lu, ustar_lu, diffc, &
    42454065                     Rb )
    4246 
    4247 
    4248                 !
    4249                 !-- Get Rc_tot
    4250                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    4251                      relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    4252                      r_aero_lu , Rb )
    4253 
    4254 
    4255                 !
    4256                 !-- Calculate budget
     4066!
     4067!--             Get Rc_tot
     4068                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,        &
     4069                                         glrad, cos_zenith, relh, lai, sai, nwet, lug_dep, 2,        &
     4070                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,              &
     4071                                         r_aero_lu , Rb )
     4072!
     4073!--             Calculate budget
    42574074                IF ( Rc_tot <= 0.0 )  THEN
    42584075
     
    42674084                ENDIF
    42684085
    4269 
    42704086             ENDIF
    42714087          ENDDO
    42724088       ENDIF
    4273 
    4274 
    4275        !
    4276        !-- Windows
     4089!
     4090!--    Windows
    42774091       IF ( surf_usm_h%frac(ind_wat_win,m) > 0 )  THEN
    4278 
    4279 
    4280           !
    4281           !-- No vegetation on windows:
     4092!
     4093!--       No vegetation on windows:
    42824094          lai = 0.0_wp
    42834095          sai = 0.0_wp
    42844096
    42854097          slinnfac = 1.0_wp
    4286 
    4287           !
    4288           !-- Get vd
     4098!
     4099!--       Get vd
    42894100          DO  lsp = 1, nvar
    4290              !
    4291              !-- Initialize
     4101!
     4102!--          Initialize
    42924103             vs = 0.0_wp
    42934104             vd_lu = 0.0_wp
     
    42974108             IF ( spc_names(lsp) == 'PM10' )  THEN
    42984109                part_type = 1
    4299                 !
    4300                 !-- Sedimentation velocity
     4110!
     4111!--            Sedimentation velocity
    43014112                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    43024113                     particle_pars(ind_p_size, part_type), &
     
    43044115                     visc)
    43054116
    4306                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
    4307                      vs, &
     4117                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
     4118                     particle_pars(ind_p_size, part_type), &
     4119                     particle_pars(ind_p_slip, part_type), &
     4120                     nwet, ttemp, dens, visc,              &
     4121                     lud_dep, r_aero_lu, ustar_lu )
     4122
     4123                bud_now_lud(lsp) = - cc(lsp) * &
     4124                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     4125
     4126             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
     4127                part_type = 2
     4128!
     4129!--             Sedimentation velocity
     4130                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     4131                     particle_pars(ind_p_size, part_type), &
     4132                     particle_pars(ind_p_slip, part_type), &
     4133                     visc)
     4134
     4135                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
    43084136                     particle_pars(ind_p_size, part_type), &
    43094137                     particle_pars(ind_p_slip, part_type), &
     
    43154143                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    43164144
    4317 
    4318              ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
    4319                 part_type = 2
    4320                 !
    4321                 !-- Sedimentation velocity
    4322                 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
    4323                      particle_pars(ind_p_size, part_type), &
    4324                      particle_pars(ind_p_slip, part_type), &
    4325                      visc)
    4326 
    4327 
    4328                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
    4329                      vs, &
    4330                      particle_pars(ind_p_size, part_type), &
    4331                      particle_pars(ind_p_slip, part_type), &
    4332                      nwet, ttemp, dens, visc, &
    4333                      lud_dep, &
    4334                      r_aero_lu, ustar_lu )
    4335 
    4336                 bud_now_lud(lsp) = - cc(lsp) * &
    4337                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    4338 
    4339 
    43404145             ELSE  !< GASES
    4341 
    4342                 !
    4343                 !-- Read spc_name of current species for gas PARAMETER
     4146!
     4147!--             Read spc_name of current species for gas PARAMETER
    43444148
    43454149                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
     
    43514155                   ENDDO
    43524156                ELSE
    4353                    ! For now species not deposited
     4157!
     4158!--                For now species not deposited
    43544159                   CYCLE
    43554160                ENDIF
    4356 
    4357                 !
    4358                 !-- Factor used for conversion from ppb to ug/m3 :
    4359                 !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
    4360                 !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
    4361                 !   c           1e-9               xm_tracer         1e9       /       xm_air            dens
    4362                 !-- thus:
    4363                 !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
    4364                 !-- Use density at k:
     4161!
     4162!--             Factor used for conversion from ppb to ug/m3 :
     4163!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     4164!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     4165!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
     4166!--             thus:
     4167!--                  c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     4168!--             Use density at k:
    43654169
    43664170                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
    4367 
    4368                 !
    4369                 !-- Atmospheric concentration in DEPAC is requested in ug/m3:
    4370                 !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     4171!
     4172!--             Atmospheric concentration in DEPAC is requested in ug/m3:
     4173!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    43714174                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
    4372 
    4373                 !
    4374                 !-- Diffusivity for DEPAC relevant gases
    4375                 !-- Use default value
    4376                 diffc            = 0.11e-4
    4377                 !
    4378                 !-- overwrite with known coefficients of diffusivity from Massman (1998)
     4175!
     4176!--             Diffusivity for DEPAC relevant gases
     4177!--             Use default value
     4178                diffc = 0.11e-4
     4179!
     4180!--             overwrite with known coefficients of diffusivity from Massman (1998)
    43794181                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    43804182                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     
    43844186                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    43854187                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    4386 
    4387 
    4388                 !
    4389                 !-- Get quasi-laminar boundary layer resistance Rb:
    4390                 CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland), &
    4391                      z0h_lu, ustar_lu, diffc, &
    4392                      Rb )
    4393 
    4394                 !
    4395                 !-- Get Rc_tot
    4396                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    4397                      relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    4398                      r_aero_lu , Rb )
    4399 
    4400 
    4401                 !
    4402                 !-- Calculate budget
     4188!
     4189!--             Get quasi-laminar boundary layer resistance Rb:
     4190                CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland),   &
     4191                     z0h_lu, ustar_lu, diffc, rb )
     4192!
     4193!--             Get Rc_tot
     4194                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,      &
     4195                                         glrad, cos_zenith, relh, lai, sai, nwet, lud_dep, 2,      &
     4196                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,            &
     4197                                         r_aero_lu , rb )
     4198!
     4199!--             Calculate budget
    44034200                IF ( Rc_tot <= 0.0 )  THEN
    44044201
     
    44194216
    44204217       bud_now = 0.0_wp
    4421        !
    4422        !-- Calculate overall budget for surface m and adapt concentration
     4218!
     4219!--    Calculate overall budget for surface m and adapt concentration
    44234220       DO  lsp = 1, nspec
    44244221
     
    44274224               surf_usm_h%frac(ind_pav_green,m) * bud_now_lug(lsp) + &
    44284225               surf_usm_h%frac(ind_wat_win,m) * bud_now_lud(lsp)
    4429 
    4430           !
    4431           !-- Compute new concentration
     4226!
     4227!--       Compute new concentration
    44324228          cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh
    44334229
     
    44394235
    44404236
    4441 
    44424237 END SUBROUTINE chem_depo
    44434238
    44444239
    4445 
    4446  !----------------------------------------------------------------------------------
    4447  !
    4448  !> DEPAC:
    4449  !> Code of the DEPAC routine and corresponding subroutines below from the DEPAC
    4450  !> module of the LOTOS-EUROS model (Manders et al., 2017)
    4451  !   
    4452  !> Original DEPAC routines by RIVM and TNO (2015), for Documentation see
    4453  !> van Zanten et al., 2010.
    4454  !---------------------------------------------------------------------------------
    4455  !   
    4456  !----------------------------------------------------------------------------------   
    4457  !
    4458  !> depac: compute total canopy (or surface) resistance Rc for gases
    4459  !----------------------------------------------------------------------------------
    4460 
     4240!------------------------------------------------------------------------------!
     4241! Description:
     4242! ------------
     4243!> Subroutine to compute total canopy (or surface) resistance Rc for gases
     4244!>
     4245!> DEPAC:
     4246!> Code of the DEPAC routine and corresponding subroutines below from the DEPAC
     4247!> module of the LOTOS-EUROS model (Manders et al., 2017)
     4248!>
     4249!> Original DEPAC routines by RIVM and TNO (2015), for Documentation see
     4250!> van Zanten et al., 2010.
     4251!------------------------------------------------------------------------------!
    44614252 SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, glrad, sinphi,  &
    44624253      rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, catm, diffc,           &
    44634254      ra, rb ) 
    4464 
    4465 
    4466    !
    4467    !--Some of depac arguments are OPTIONAL:
    4468    !
    4469    !-- A. compute Rc_tot without compensation points (ccomp_tot will be zero):
    4470    !--     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi])
    4471    !
    4472    !-- B. compute Rc_tot with compensation points (used for LOTOS-EUROS):
    4473    !--     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
    4474    !                  c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
    4475    !
    4476    !-- C. compute effective Rc based on compensation points (used for OPS):
    4477    !--     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
    4478    !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
    4479    !                 ra, rb, rc_eff)
    4480    !-- X1. Extra (OPTIONAL) output variables:
    4481    !--     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
    4482    !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
    4483    !                 ra, rb, rc_eff, &
    4484    !                 gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out)
    4485    !-- X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves):
    4486    !--     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
    4487    !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
    4488    !                 ra, rb, rc_eff, &
    4489    !                 gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, &
    4490    !                 calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
    4491 
     4255!
     4256!--   Some of depac arguments are OPTIONAL:
     4257!--    A. compute Rc_tot without compensation points (ccomp_tot will be zero):
     4258!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi])
     4259!--    B. compute Rc_tot with compensation points (used for LOTOS-EUROS):
     4260!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
     4261!--                c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
     4262!--
     4263!--    C. compute effective Rc based on compensation points (used for OPS):
     4264!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
     4265!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
     4266!--               ra, rb, rc_eff)
     4267!--    X1. Extra (OPTIONAL) output variables:
     4268!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
     4269!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
     4270!--               ra, rb, rc_eff, &
     4271!--               gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out)
     4272!--    X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves):
     4273!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
     4274!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
     4275!--               ra, rb, rc_eff, &
     4276!--               gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, &
     4277!--               calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
    44924278
    44934279    IMPLICIT NONE
    44944280
    4495     CHARACTER(len=*), INTENT(IN) ::  compnam         !< component name
     4281    CHARACTER(LEN=*), INTENT(IN) ::  compnam         !< component name
    44964282                                                     !< 'HNO3','NO','NO2','O3','SO2','NH3'
    44974283    INTEGER(iwp), INTENT(IN) ::  day_of_year         !< day of year, 1 ... 365 (366)
     
    45184304    REAL(wp), INTENT(OUT) ::  rc_tot                 !< total canopy resistance Rc (s/m)
    45194305    REAL(wp), INTENT(OUT) ::  ccomp_tot              !< total compensation point (ug/m3)
    4520                                                      !< [= 0 for species that don't have a compensation
    4521     !
    4522     !-- Local variables:
    4523     !
    4524     !-- Component number taken from component name, paramteres matched with include files
     4306!                                                     !< [= 0 for species that don't have a compensation
     4307!-- Local variables:
     4308!
     4309!-- Component number taken from component name, paramteres matched with include files
    45254310    INTEGER(iwp) ::  icmp
    4526 
    4527     !
    4528     !-- Component numbers:
     4311!
     4312!-- Component numbers:
    45294313    INTEGER(iwp), PARAMETER ::  icmp_o3   = 1
    45304314    INTEGER(iwp), PARAMETER ::  icmp_so2  = 2
     
    45414325    !< = 1 -> constant Rc
    45424326    !< = 2 -> temperature dependent Rc
    4543     !
    4544     !-- Vegetation indicators:
     4327!
     4328!-- Vegetation indicators:
    45454329    LOGICAL ::  lai_present                          !< leaves are present for current land use type
    45464330    LOGICAL ::  sai_present                          !< vegetation is present for current land use type
     
    45544338    REAL(wp) ::  cstom                               !< stomatal compensation point (ug/m3)
    45554339    REAL(wp) ::  csoil                               !< soil compensation point (ug/m3)
    4556 
    4557 
    45584340!
    45594341!-- Next statement is just to avoid compiler warning about unused variable
    45604342    IF ( day_of_year == 0  .OR.  ( catm + lat + ra + rb ) > 0.0_wp )  CONTINUE
    4561 
    4562 
    4563     !
    4564     !-- Define component number
     4343!
     4344!-- Define component number
    45654345    SELECT CASE ( TRIM( compnam ) )
    45664346
     
    45964376
    45974377    CASE default
    4598        !
    4599        !-- Component not part of DEPAC --> not deposited
     4378!
     4379!--    Component not part of DEPAC --> not deposited
    46004380       RETURN
    46014381
    46024382    END SELECT
    46034383
    4604     !
    4605     !-- Inititalize
     4384!
     4385!-- Inititalize
    46064386    gw        = 0.0_wp
    46074387    gstom     = 0.0_wp
     
    46114391    cstom     = 0.0_wp
    46124392    csoil     = 0.0_wp
    4613 
    4614 
    4615     !
    4616     !-- Check whether vegetation is present:
     4393!
     4394!-- Check whether vegetation is present:
    46174395    lai_present = ( lai > 0.0 )
    46184396    sai_present = ( sai > 0.0 )
    4619 
    4620     !
    4621     !-- Set Rc (i.e. rc_tot) in special cases:
     4397!
     4398!-- Set Rc (i.e. rc_tot) in special cases:
    46224399    CALL rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot )
    4623 
    4624     !
    4625     !-- If Rc is not set:
     4400!
     4401!-- If Rc is not set:
    46264402    IF ( .NOT. ready ) then
    4627 
    4628        !
    4629        !-- External conductance:
     4403!
     4404!--    External conductance:
    46304405       CALL rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai,gw )         
    4631 
    4632        !
    4633        !-- Stomatal conductance:
     4406!
     4407!--    Stomatal conductance:
    46344408       CALL rc_gstom( icmp, compnam, lu, lai_present, lai, glrad, sinphi, t, rh, diffc, gstom, p )
    4635        !
    4636        !-- Effective soil conductance:
     4409!
     4410!--    Effective soil conductance:
    46374411       CALL rc_gsoil_eff( icmp, lu, sai, ust, nwet, t, gsoil_eff )
    4638 
    4639        !
    4640        !-- Total canopy conductance (gc_tot) and resistance Rc (rc_tot):
     4412!
     4413!--    Total canopy conductance (gc_tot) and resistance Rc (rc_tot):
    46414414       CALL rc_rctot( gstom, gsoil_eff, gw, gc_tot, rc_tot )
    4642 
    4643        !
    4644        !-- Needed to include compensation point for NH3
    4645        !-- Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3):
    4646        !-- CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,&
    4647        !     lai_present, sai_present, &
    4648        !     ccomp_tot, &
    4649        !     catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, &
    4650        !     c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, &
    4651        !     tsea=tsea,cw=cw,cstom=cstom,csoil=csoil )
    4652        !
    4653        !-- Effective Rc based on compensation points:
    4654        !   IF ( present(rc_eff) ) then
    4655        !     check on required arguments:
    4656        !      IF ( (.not. present(catm)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
    4657        !         stop 'output argument rc_eff requires input arguments catm, ra and rb'
    4658        !      END IF
    4659        !--    Compute rc_eff :
     4415!
     4416!--    Needed to include compensation point for NH3
     4417!--    Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3):
     4418!--    CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,&
     4419!--          lai_present, sai_present, &
     4420!--          ccomp_tot, &
     4421!--          catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, &
     4422!--          c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, &
     4423!--          tsea=tsea,cw=cw,cstom=cstom,csoil=csoil )
     4424!
     4425!--    Effective Rc based on compensation points:
     4426!--        IF ( present(rc_eff) ) then
     4427!--          check on required arguments:
     4428!--           IF ( (.not. present(catm)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
     4429!--              stop 'output argument rc_eff requires input arguments catm, ra and rb'
     4430!--           END IF
     4431!
     4432!--       Compute rc_eff :
    46604433       !      CALL rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff)
    46614434       !   ENDIF
     
    46654438
    46664439
    4667 
    4668  !-------------------------------------------------------------------
    4669  !> rc_special: compute total canopy resistance in special CASEs
    4670  !-------------------------------------------------------------------
     4440!------------------------------------------------------------------------------!
     4441! Description:
     4442! ------------
     4443!> Subroutine to compute total canopy resistance in special cases
     4444!------------------------------------------------------------------------------!
    46714445 SUBROUTINE rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot )
    46724446
    46734447    IMPLICIT NONE
    46744448   
    4675     CHARACTER(len=*), INTENT(IN) ::  compnam     !< component name
    4676 
    4677     INTEGER(iwp), INTENT(IN) ::  icmp            !< component index     
    4678     INTEGER(iwp), INTENT(IN) ::  lu              !< land use type, lu = 1,...,nlu
    4679     INTEGER(iwp), INTENT(IN) ::  nwet            !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
    4680 
    4681     REAL(wp), INTENT(IN) ::  t                   !< temperature (C)
     4449    CHARACTER(LEN=*), INTENT(IN) ::  compnam     !< component name
     4450
     4451    INTEGER(iwp), INTENT(IN)  ::  icmp            !< component index     
     4452    INTEGER(iwp), INTENT(IN)  ::  lu              !< land use type, lu = 1,...,nlu
     4453    INTEGER(iwp), INTENT(IN)  ::  nwet            !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     4454
     4455    REAL(wp), INTENT(IN)  ::  t                   !< temperature (C)
    46824456
    46834457    REAL(wp), INTENT(OUT) ::  rc_tot             !< total canopy resistance Rc (s/m)
     
    46864460    LOGICAL, INTENT(OUT) ::  ready               !< Rc has been set
    46874461                                                 !< = 1 -> constant Rc
    4688                                                  !< = 2 -> temperature dependent Rc
    4689 
    46904462!
    46914463!-- Next line is to avoid compiler warning about unused variable
    46924464    IF ( icmp == 0 )  CONTINUE
    4693 
    4694     !
    4695     !-- rc_tot is not yet set:
     4465!
     4466!-- rc_tot is not yet set:
    46964467    ready = .FALSE.
    4697 
    4698     !
    4699     !-- Default compensation point in special CASEs = 0:
     4468!
     4469!-- Default compensation point in special CASEs = 0:
    47004470    ccomp_tot = 0.0_wp
    47014471
    47024472    SELECT CASE( TRIM( compnam ) )
    47034473    CASE( 'HNO3', 'N2O5', 'NO3', 'H2O2' )
    4704        !
    4705        !-- No separate resistances for HNO3; just one total canopy resistance:
     4474!
     4475!--    No separate resistances for HNO3; just one total canopy resistance:
    47064476       IF ( t < -5.0_wp .AND. nwet == 9 )  THEN
    4707           !
    4708           !-- T < 5 C and snow:
     4477!
     4478!--      T < 5 C and snow:
    47094479          rc_tot = 50.0_wp
    47104480       ELSE
    4711           !
    4712           !-- all other circumstances:
     4481!
     4482!--      all other circumstances:
    47134483          rc_tot = 10.0_wp
    47144484       ENDIF
     
    47244494       ENDIF
    47254495    CASE( 'NO2', 'O3', 'SO2', 'NH3' )
    4726        !
    4727        !-- snow surface:
     4496!
     4497!--    snow surface:
    47284498       IF ( nwet == 9 )  THEN
    4729           !
    4730           !-- To be activated when snow is implemented
     4499!
     4500!--      To be activated when snow is implemented
    47314501          !CALL rc_snow(ipar_snow(icmp),t,rc_tot)
    47324502          ready = .TRUE.
     
    47404510
    47414511
    4742 
    4743  !-------------------------------------------------------------------
    4744  !> rc_gw: compute external conductance
    4745  !-------------------------------------------------------------------
     4512!------------------------------------------------------------------------------!
     4513! Description:
     4514! ------------
     4515!> Subroutine to compute external conductance
     4516!------------------------------------------------------------------------------!
    47464517 SUBROUTINE rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai, gw )
    47474518
    47484519    IMPLICIT NONE
    4749 
    4750     !
    4751     !-- Input/output variables:
    4752     CHARACTER(len=*), INTENT(IN) ::  compnam      !< component name
     4520!
     4521!-- Input/output variables:
     4522    CHARACTER(LEN=*), INTENT(IN) ::  compnam      !< component name
    47534523
    47544524    INTEGER(iwp), INTENT(IN) ::  nwet             !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     
    47574527                                                  !< iratns = 2: high NH3/SO2
    47584528                                                  !< iratns = 3: very low NH3/SO2
    4759 
    47604529    LOGICAL, INTENT(IN) ::  sai_present
    47614530
     
    47664535    REAL(wp), INTENT(OUT) ::  gw                  !< external leaf conductance (m/s)
    47674536
    4768 
    47694537    SELECT CASE( TRIM( compnam ) )
    47704538
     
    47834551    CASE( 'NH3' )
    47844552       CALL rw_nh3_sutton( t, rh, sai_present, gw )
    4785 
    4786        !
    4787        !-- conversion from leaf resistance to canopy resistance by multiplying with sai:
     4553!
     4554!--    conversion from leaf resistance to canopy resistance by multiplying with sai:
    47884555       gw = sai * gw
    47894556
    47904557    CASE default
    4791        message_string = 'Component '// TRIM(compnam) // ' not supported'
     4558       message_string = 'Component '// TRIM( compnam ) // ' not supported'
    47924559       CALL message( 'rc_gw', 'CM0458', 1, 2, 0, 6, 0 )
    47934560    END SELECT
     
    47964563
    47974564
    4798 
    4799  !-------------------------------------------------------------------
    4800  !> rw_so2: compute external leaf conductance for SO2
    4801  !-------------------------------------------------------------------
     4565!------------------------------------------------------------------------------!
     4566! Description:
     4567! ------------
     4568!> Subroutine to compute external leaf conductance for SO2
     4569!------------------------------------------------------------------------------!
    48024570 SUBROUTINE rw_so2( t, nwet, rh, iratns, sai_present, gw )
    48034571
    48044572    IMPLICIT NONE
    4805 
    4806     !
    4807     !-- Input/output variables:
     4573!
     4574!-- Input/output variables:
    48084575    INTEGER(iwp), INTENT(IN) ::  nwet        !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
    48094576    INTEGER(iwp), INTENT(IN) ::  iratns      !< index for NH3/SO2 ratio:
     
    48174584
    48184585    REAL(wp), INTENT(OUT) ::  gw             !< external leaf conductance (m/s)
    4819 
    4820     !
    4821     !-- Local variables:
     4586!
     4587!-- Local variables:
    48224588    REAL(wp) ::  rw                          !< external leaf resistance (s/m)
    4823 
    4824 
    4825     !
    4826     !-- Check if vegetation present:
     4589!
     4590!-- Check if vegetation present:
    48274591    IF ( sai_present )  THEN
    48284592
    48294593       IF ( nwet == 0 )  THEN
    4830           !--------------------------
    4831           ! dry surface
    4832           !--------------------------
    4833           ! T > -1 C
     4594!
     4595!--   ------------------------
     4596!--         dry surface
     4597!--   ------------------------
     4598!--         T > -1 C
    48344599          IF ( t > -1.0_wp )  THEN
    48354600             IF ( rh < 81.3_wp )  THEN
     
    48484613          ENDIF
    48494614       ELSE
    4850           !--------------------------
    4851           ! wet surface
    4852           !--------------------------
     4615!
     4616!--   ------------------------
     4617!--         wet surface
     4618!--   ------------------------
    48534619          rw = 10.0_wp !see Table 5, Erisman et al, 1994 Atm. Environment, 0 is impl. as 10
    48544620       ENDIF
    4855 
    4856        ! very low NH3/SO2 ratio:
     4621!
     4622!--    very low NH3/SO2 ratio:
    48574623       IF ( iratns == iratns_very_low ) rw = rw + 50.0_wp
    4858 
    4859        ! Conductance:
     4624!
     4625!--      Conductance:
     4626       gw = 1.0_wp / rw
     4627    ELSE
     4628!
     4629!--      no vegetation:
     4630       gw = 0.0_wp
     4631    ENDIF
     4632
     4633 END SUBROUTINE rw_so2
     4634
     4635
     4636!------------------------------------------------------------------------------!
     4637! Description:
     4638! ------------
     4639!> Subroutine to compute external leaf conductance for NH3,
     4640!>                  following Sutton & Fowler, 1993
     4641!------------------------------------------------------------------------------!
     4642 SUBROUTINE rw_nh3_sutton( tsurf, rh,sai_present, gw )
     4643
     4644    IMPLICIT NONE
     4645!
     4646!-- Input/output variables:
     4647    LOGICAL, INTENT(IN) ::  sai_present
     4648
     4649    REAL(wp), INTENT(IN) ::  tsurf          !< surface temperature (C)
     4650    REAL(wp), INTENT(IN) ::  rh             !< relative humidity (%)
     4651
     4652    REAL(wp), INTENT(OUT) ::  gw            !< external leaf conductance (m/s)
     4653!
     4654!-- Local variables:
     4655    REAL(wp) ::  rw                         !< external leaf resistance (s/m)
     4656    REAL(wp) ::  sai_grass_haarweg          !< surface area index at experimental site Haarweg
     4657!
     4658!-- Fix sai_grass at value valid for Haarweg data for which gamma_w parametrization is derived
     4659    sai_grass_haarweg = 3.5_wp
     4660!
     4661!-- Calculation rw:
     4662!--                    100 - rh
     4663!--    rw = 2.0 * exp(----------)
     4664!--                      12
     4665
     4666    IF ( sai_present )  THEN
     4667!
     4668!--    External resistance according to Sutton & Fowler, 1993
     4669       rw = 2.0_wp * exp( ( 100.0_wp - rh ) / 12.0_wp )
     4670       rw = sai_grass_haarweg * rw
     4671!
     4672!--    Frozen soil (from Depac v1):
     4673       IF ( tsurf < 0.0_wp ) rw = 200.0_wp
     4674!
     4675!--    Conductance:
    48604676       gw = 1.0_wp / rw
    48614677    ELSE
     
    48644680    ENDIF
    48654681
    4866  END SUBROUTINE rw_so2
    4867 
    4868 
    4869 
    4870  !-------------------------------------------------------------------
    4871  !> rw_nh3_sutton: compute external leaf conductance for NH3,
    4872  !>                  following Sutton & Fowler, 1993
    4873  !-------------------------------------------------------------------
    4874  SUBROUTINE rw_nh3_sutton( tsurf, rh,sai_present, gw )
     4682 END SUBROUTINE rw_nh3_sutton
     4683
     4684
     4685!------------------------------------------------------------------------------!
     4686! Description:
     4687! ------------
     4688!> Subroutine to compute external leaf conductance
     4689!------------------------------------------------------------------------------!
     4690 SUBROUTINE rw_constant( rw_val, sai_present, gw )
    48754691
    48764692    IMPLICIT NONE
    4877 
    4878     !
    4879     !-- Input/output variables:
     4693!
     4694!-- Input/output variables:
    48804695    LOGICAL, INTENT(IN) ::  sai_present
    48814696
    4882     REAL(wp), INTENT(IN) ::  tsurf          !< surface temperature (C)
    4883     REAL(wp), INTENT(IN) ::  rh             !< relative humidity (%)
    4884 
    4885     REAL(wp), INTENT(OUT) ::  gw            !< external leaf conductance (m/s)
    4886 
    4887     !
    4888     !-- Local variables:
    4889     REAL(wp) ::  rw                         !< external leaf resistance (s/m)
    4890     REAL(wp) ::  sai_grass_haarweg          !< surface area index at experimental site Haarweg
    4891 
    4892     !
    4893     !-- Fix sai_grass at value valid for Haarweg data for which gamma_w parametrization is derived
    4894     sai_grass_haarweg = 3.5_wp
    4895 
    4896     !
    4897     !-- Calculation rw:
    4898     !                  100 - rh
    4899     !  rw = 2.0 * exp(----------)
    4900     !                    12
    4901 
    4902     IF ( sai_present )  THEN
    4903 
    4904        ! External resistance according to Sutton & Fowler, 1993
    4905        rw = 2.0_wp * exp( ( 100.0_wp - rh ) / 12.0_wp )
    4906        rw = sai_grass_haarweg * rw
    4907 
    4908        ! Frozen soil (from Depac v1):
    4909        IF ( tsurf < 0.0_wp ) rw = 200.0_wp
    4910 
    4911        ! Conductance:
    4912        gw = 1.0_wp / rw
    4913     ELSE
    4914        ! no vegetation:
    4915        gw = 0.0_wp
    4916     ENDIF
    4917 
    4918  END SUBROUTINE rw_nh3_sutton
    4919 
    4920 
    4921 
    4922  !-------------------------------------------------------------------
    4923  !> rw_constant: compute constant external leaf conductance
    4924  !-------------------------------------------------------------------
    4925  SUBROUTINE rw_constant( rw_val, sai_present, gw )
    4926 
    4927     IMPLICIT NONE
    4928 
    4929     !
    4930     !-- Input/output variables:
    4931     LOGICAL, INTENT(IN) ::  sai_present
    4932 
    49334697    REAL(wp), INTENT(IN) ::  rw_val       !< constant value of Rw   
    49344698
    49354699    REAL(wp), INTENT(OUT) ::  gw          !< wernal leaf conductance (m/s)
    4936 
    4937     !
    4938     !-- Compute conductance:
     4700!
     4701!-- Compute conductance:
    49394702    IF ( sai_present .AND. .NOT.missing(rw_val) )  THEN
    49404703       gw = 1.0_wp / rw_val
     
    49464709
    49474710
    4948 
    4949  !-------------------------------------------------------------------
    4950  !> rc_gstom: compute stomatal conductance
    4951  !-------------------------------------------------------------------
     4711!------------------------------------------------------------------------------!
     4712! Description:
     4713! ------------
     4714!> Subroutine to compute stomatal conductance
     4715!------------------------------------------------------------------------------!
    49524716 SUBROUTINE rc_gstom( icmp, compnam, lu, lai_present, lai, glrad, sinphi, t, rh, diffc, gstom, p )
    49534717
    49544718    IMPLICIT NONE
    4955 
    4956     !
    4957     !-- input/output variables:
    4958     CHARACTER(len=*), INTENT(IN) ::  compnam       !< component name
     4719!
     4720!-- input/output variables:
     4721    CHARACTER(LEN=*), INTENT(IN) ::  compnam       !< component name
    49594722
    49604723    INTEGER(iwp), INTENT(IN) ::  icmp              !< component index
     
    49734736
    49744737    REAL(wp), INTENT(OUT) ::  gstom                !< stomatal conductance (m/s)
    4975 
    4976     !
    4977     !-- Local variables
     4738!
     4739!-- Local variables
    49784740    REAL(wp) ::  vpd                               !< vapour pressure deficit (kPa)
    49794741
    49804742    REAL(wp), PARAMETER ::  dO3 = 0.13e-4          !< diffusion coefficient of ozon (m2/s)
    4981 
    49824743!
    49834744!-- Next line is to avoid compiler warning about unused variables
    49844745    IF ( icmp == 0 )  CONTINUE
    49854746
    4986     SELECT CASE( TRIM(compnam) )
     4747    SELECT CASE( TRIM( compnam ) )
    49874748
    49884749    CASE( 'NO', 'CO' )
    4989        !
    4990        !-- For no stomatal uptake is neglected:
     4750!
     4751!--    For no stomatal uptake is neglected:
    49914752       gstom = 0.0_wp
    49924753
    49934754    CASE( 'NO2', 'O3', 'SO2', 'NH3' )
    4994 
    4995        !
    4996        !-- if vegetation present:
     4755!
     4756!--    if vegetation present:
    49974757       IF ( lai_present )  THEN
    49984758
     
    50054765          ENDIF
    50064766       ELSE
    5007           !
    5008           !--no vegetation; zero conductance (infinite resistance):
     4767!
     4768!--       no vegetation; zero conductance (infinite resistance):
    50094769          gstom = 0.0_wp
    50104770       ENDIF
    50114771
    50124772    CASE default
    5013        message_string = 'Component '// TRIM(compnam) // ' not supported'
     4773       message_string = 'Component '// TRIM( compnam ) // ' not supported'
    50144774       CALL message( 'rc_gstom', 'CM0459', 1, 2, 0, 6, 0 )
    50154775    END SELECT
     
    50184778
    50194779
    5020 
    5021  !-------------------------------------------------------------------
    5022  !> rc_gstom_emb: stomatal conductance according to Emberson
    5023  !-------------------------------------------------------------------
     4780!------------------------------------------------------------------------------!
     4781! Description:
     4782! ------------
     4783!> Subroutine to compute stomatal conductance according to Emberson
     4784!------------------------------------------------------------------------------!
    50244785 SUBROUTINE rc_gstom_emb( lu, glrad, T, vpd, lai_present, lai, sinp, Gsto, p )
    5025     !
    5026     !-- History
    5027     !>   Original code from Lotos-Euros, TNO, M. Schaap
    5028     !>   2009-08, M.C. van Zanten, Rivm
    5029     !>     Updated and extended.
    5030     !>   2009-09, Arjo Segers, TNO
    5031     !>     Limitted temperature influence to range to avoid
    5032     !>     floating point exceptions.
    5033     !
    5034     !> Method
    5035     !
    5036     !>   Code based on Emberson et al, 2000, Env. Poll., 403-413
    5037     !>   Notation conform Unified EMEP Model Description Part 1, ch 8
    5038     !
    5039     !>   In the calculation of f_light the modification of L. Zhang 2001, AE to the PARshade and PARsun
    5040     !>   parametrizations of Norman 1982 are applied
    5041     !
    5042     !>   f_phen and f_SWP are set to 1
    5043     !
    5044     !>   Land use types DEPAC versus Emberson (Table 5.1, EMEP model description)
    5045     !>   DEPAC                     Emberson
    5046     !>     1 = grass                 GR = grassland
    5047     !>     2 = arable land           TC = temperate crops ( lai according to RC = rootcrops)
    5048     !<     3 = permanent crops       TC = temperate crops ( lai according to RC = rootcrops)
    5049     !<     4 = coniferous forest     CF = tempareate/boREAL(wp) coniferous forest
    5050     !>     5 = deciduous forest      DF = temperate/boREAL(wp) deciduous forest
    5051     !>     6 = water                 W  = water
    5052     !>     7 = urban                 U  = urban
    5053     !>     8 = other                 GR = grassland
    5054     !>     9 = desert                DE = desert
     4786!
     4787!>  History
     4788!>   Original code from Lotos-Euros, TNO, M. Schaap
     4789!>   2009-08, M.C. van Zanten, Rivm
     4790!>     Updated and extended.
     4791!>   2009-09, Arjo Segers, TNO
     4792!>     Limitted temperature influence to range to avoid
     4793!>     floating point exceptions.
     4794
     4795!> Method
     4796
     4797!>   Code based on Emberson et al, 2000, Env. Poll., 403-413
     4798!>   Notation conform Unified EMEP Model Description Part 1, ch 8
     4799!
     4800!>   In the calculation of f_light the modification of L. Zhang 2001, AE to the PARshade and PARsun
     4801!>   parametrizations of Norman 1982 are applied
     4802!>   f_phen and f_SWP are set to 1
     4803!
     4804!>   Land use types DEPAC versus Emberson (Table 5.1, EMEP model description)
     4805!>   DEPAC                     Emberson
     4806!>     1 = grass                 GR = grassland
     4807!>     2 = arable land           TC = temperate crops ( lai according to RC = rootcrops)
     4808!>     3 = permanent crops       TC = temperate crops ( lai according to RC = rootcrops)
     4809!>     4 = coniferous forest     CF = tempareate/boREAL(wp) coniferous forest
     4810!>     5 = deciduous forest      DF = temperate/boREAL(wp) deciduous forest
     4811!>     6 = water                 W  = water
     4812!>     7 = urban                 U  = urban
     4813!>     8 = other                 GR = grassland
     4814!>     9 = desert                DE = desert
    50554815
    50564816    IMPLICIT NONE
    5057 
    5058     !
    5059     !-- Emberson specific declarations
    5060 
    5061     !
    5062     !-- Input/output variables:
     4817!
     4818!-- Emberson specific declarations
     4819!
     4820!-- Input/output variables:
    50634821    INTEGER(iwp), INTENT(IN) ::  lu             !< land use type, lu = 1,...,nlu
    50644822
     
    50754833
    50764834    REAL(wp), INTENT(OUT) ::  gsto              !< stomatal conductance (m/s)
    5077 
    5078     !
    5079     !-- Local variables:
     4835!
     4836!-- Local variables:
    50804837    REAL(wp) ::  f_light
    50814838    REAL(wp) ::  f_phen
     
    50944851    REAL(wp) ::  pres
    50954852    REAL(wp), PARAMETER ::  p_sealevel = 1.01325e05    !< Pa
    5096 
    5097 
    5098     !
    5099     !-- Check whether vegetation is present:
     4853!
     4854!-- Check whether vegetation is present:
    51004855    IF ( lai_present )  THEN
    51014856
     
    51064861          sinphi = sinp
    51074862       END IF
    5108 
    5109        !
    5110        !-- ratio between actual and sea-level pressure is used
    5111        !-- to correct for height in the computation of par;
    5112        !-- should not exceed sea-level pressure therefore ...
     4863!
     4864!--    ratio between actual and sea-level pressure is used
     4865!--    to correct for height in the computation of par;
     4866!--    should not exceed sea-level pressure therefore ...
    51134867       IF (  present(p) )  THEN
    51144868          pres = min( p, p_sealevel )
     
    51164870          pres = p_sealevel
    51174871       ENDIF
    5118 
    5119        !
    5120        !-- direct and diffuse par, Photoactive (=visible) radiation:
     4872!
     4873!--    direct and diffuse par, Photoactive (=visible) radiation:
    51214874       CALL par_dir_diff( glrad, sinphi, pres, p_sealevel, pardir, pardiff )
    5122 
    5123        !
    5124        !-- par for shaded leaves (canopy averaged):
     4875!
     4876!--    par for shaded leaves (canopy averaged):
    51254877       parshade = pardiff * exp( -0.5 * lai**0.7 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi )     !< Norman,1982
    51264878       IF ( glrad > 200.0_wp .AND. lai > 2.5_wp )  THEN
    51274879          parshade = pardiff * exp( -0.5 * lai**0.8 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi )  !< Zhang et al., 2001
    51284880       END IF
    5129 
    5130        !
    5131        !-- par for sunlit leaves (canopy averaged):
    5132        !-- alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5
     4881!
     4882!--    par for sunlit leaves (canopy averaged):
     4883!--    alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5
    51334884       parsun = pardir * 0.5/sinphi + parshade             !< Norman, 1982
    51344885       IF ( glrad > 200.0_wp .AND. lai > 2.5_wp )  THEN
    51354886          parsun = pardir**0.8 * 0.5 / sinphi + parshade   !< Zhang et al., 2001
    51364887       END IF
    5137 
    5138        !
    5139        !-- leaf area index for sunlit and shaded leaves:
     4888!
     4889!--    leaf area index for sunlit and shaded leaves:
    51404890       IF ( sinphi > 0 )  THEN
    51414891          laisun = 2 * sinphi * ( 1 - exp( -0.5 * lai / sinphi ) )
     
    51504900
    51514901       f_light = MAX(f_light,f_min(lu))
    5152 
    5153        !
    5154        !-- temperature influence; only non-zero within range [tmin,tmax]:
     4902!
     4903!--    temperature influence; only non-zero within range [tmin,tmax]:
    51554904       IF ( ( tmin(lu) < t ) .AND. ( t < tmax(lu) ) )  THEN
    51564905          bt = ( tmax(lu) - topt(lu) ) / ( topt(lu) - tmin(lu) )
     
    51604909       END IF
    51614910       f_temp = MAX( f_temp, f_min(lu) )
    5162 
    5163        ! vapour pressure deficit influence
     4911!
     4912!--      vapour pressure deficit influence
    51644913       f_vpd = MIN( 1.0_wp, ( ( 1.0_wp - f_min(lu) ) * ( vpd_min(lu) - vpd ) / ( vpd_min(lu) - vpd_max(lu) ) + f_min(lu) ) )
    51654914       f_vpd = MAX( f_vpd, f_min(lu) )
    51664915
    51674916       f_swp = 1.0_wp
    5168 
    5169        ! influence of phenology on stom. conductance
    5170        ! ignored for now in DEPAC since influence of f_phen on lu classes in use is negligible.
    5171        ! When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore
     4917!
     4918!--      influence of phenology on stom. conductance
     4919!--      ignored for now in DEPAC since influence of f_phen on lu classes in use is negligible.
     4920!--      When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore
    51724921       f_phen = 1.0_wp
    5173 
    5174        ! evaluate total stomatal conductance
     4922!
     4923!--      evaluate total stomatal conductance
    51754924       f_env = f_temp * f_vpd * f_swp
    51764925       f_env = MAX( f_env,f_min(lu) )
    51774926       gsto = g_max(lu) * f_light * f_phen * f_env
    5178 
    5179        ! gstom expressed per m2 leafarea;
    5180        ! this is converted with lai to m2 surface.
     4927!
     4928!--      gstom expressed per m2 leafarea;
     4929!--      this is converted with lai to m2 surface.
    51814930       gsto = lai * gsto    ! in m/s
    51824931
     
    51884937
    51894938
    5190 
    51914939 !-------------------------------------------------------------------
    51924940 !> par_dir_diff
    5193  !
    51944941 !>     Weiss, A., Norman, J.M. (1985) Partitioning solar radiation into direct and
    51954942 !>     diffuse, visible and near-infrared components. Agric. Forest Meteorol.
    51964943 !>     34, 205-213.
    5197  !
    51984944 !>     From a SUBROUTINE obtained from Leiming Zhang,
    51994945 !>     Meteorological Service of Canada
    5200  !
    52014946 !>     Leiming uses solar irradiance. This should be equal to global radiation and
    52024947 !>     Willem Asman set it to global radiation
     
    52284973    REAL(wp) ::  ww                          !< water absorption in the near infrared for 10 mm of precipitable water
    52294974
    5230 
    5231 
    5232     !
    5233     !-- Calculate visible (PAR) direct beam radiation
    5234     !-- 600 W m-2 represents average amount of par (400-700 nm wavelength)
    5235     !-- at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2)
     4975!
     4976!-- Calculate visible (PAR) direct beam radiation
     4977!-- 600 W m-2 represents average amount of par (400-700 nm wavelength)
     4978!-- at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2)
    52364979    rdu = 600.0_wp* exp( -0.185_wp * ( pres / pres_0 ) / sinphi ) * sinphi
    5237 
    5238     !
    5239     !-- Calculate potential visible diffuse radiation
     4980!
     4981!-- Calculate potential visible diffuse radiation
    52404982    rdv = 0.4_wp * ( 600.0_wp - rdu ) * sinphi
    5241 
    5242     !
    5243     !-- Calculate the water absorption in the-near infrared
     4983!
     4984!-- Calculate the water absorption in the-near infrared
    52444985    ww = 1320 * 10**( -1.195_wp + 0.4459_wp * log10( 1.0_wp / sinphi ) - 0.0345_wp * ( log10( 1.0_wp / sinphi ) )**2 )
    5245 
    5246     !
    5247     !-- Calculate potential direct beam near-infrared radiation
     4986!
     4987!-- Calculate potential direct beam near-infrared radiation
    52484988    rdm = (720.0_wp * exp(-0.06_wp * (pres / pres_0) / sinphi ) - ww ) * sinphi     !< 720 = solar constant - 600
    5249 
    5250     !
    5251     !-- Calculate potential diffuse near-infrared radiation
     4989!
     4990!-- Calculate potential diffuse near-infrared radiation
    52524991    rdn = 0.6_wp * ( 720 - rdm - ww ) * sinphi
    5253 
    5254     !
    5255     !-- Compute visible and near-infrared radiation
     4992!
     4993!-- Compute visible and near-infrared radiation
    52564994    rv = MAX( 0.1_wp, rdu + rdv )
    52574995    rn = MAX( 0.01_wp, rdm + rdn )
    5258 
    5259     !
    5260     !-- Compute ratio between input global radiation and total radiation computed here
     4996!
     4997!-- Compute ratio between input global radiation and total radiation computed here
    52614998    ratio = MIN( 0.89_wp, glrad / ( rv + rn ) )
    5262 
    5263     !
    5264     !-- Calculate total visible radiation
     4999!
     5000!-- Calculate total visible radiation
    52655001    sv = ratio * rv
    5266 
    5267     !
    5268     !-- Calculate fraction of par in the direct beam
     5002!
     5003!-- Calculate fraction of par in the direct beam
    52695004    fv = MIN( 0.99_wp, ( 0.9_wp - ratio ) / 0.7_wp )              !< help variable
    52705005    fv = MAX( 0.01_wp, rdu / rv * ( 1.0_wp - fv**0.6667_wp ) )    !< fraction of par in the direct beam
    5271 
    5272     !
    5273     !-- Compute direct and diffuse parts of par
     5006!
     5007!-- Compute direct and diffuse parts of par
    52745008    par_dir = fv * sv
    52755009    par_diff = sv - par_dir
     
    52845018
    52855019    IMPLICIT NONE
    5286 
    5287     !
    5288     !-- Input/output variables:
     5020!
     5021!-- Input/output variables:
    52895022    REAL(wp), INTENT(IN) ::  temp    !< temperature (C)
    52905023    REAL(wp), INTENT(IN) ::  relh    !< relative humidity (%)
    52915024
    52925025    REAL(wp), INTENT(OUT) ::  vpd    !< vapour pressure deficit (kPa)
    5293 
    5294     !
    5295     !-- Local variables:
     5026!
     5027!-- Local variables:
    52965028    REAL(wp) ::  esat
    5297 
    5298     !
    5299     !-- fit parameters:
     5029!
     5030!-- fit parameters:
    53005031    REAL(wp), PARAMETER ::  a1 = 6.113718e-01
    53015032    REAL(wp), PARAMETER ::  a2 = 4.43839e-02
     
    53045035    REAL(wp), PARAMETER ::  a5 = 2.16e-07
    53055036    REAL(wp), PARAMETER ::  a6 = 3.0e-09
    5306 
    5307     !
    5308     !-- esat is saturation vapour pressure (kPa) at temp(C) following Monteith(1973)
     5037!
     5038!-- esat is saturation vapour pressure (kPa) at temp(C) following Monteith(1973)
    53095039    esat = a1 + a2 * temp + a3 * temp**2 + a4 * temp**3 + a5 * temp**4 + a6 * temp**5
    53105040    vpd  = esat * ( 1 - relh / 100 )
    53115041
    53125042 END SUBROUTINE rc_get_vpd
    5313 
    53145043
    53155044
     
    53205049
    53215050    IMPLICIT NONE
    5322 
    5323     !
    5324     !-- Input/output variables:
     5051!
     5052!-- Input/output variables:
    53255053    INTEGER(iwp), INTENT(IN) ::  icmp          !< component index
    53265054    INTEGER(iwp), INTENT(IN) ::  lu            !< land use type, lu = 1,..., nlu
     
    53295057                                               !< N.B. this routine cannot be called with nwet = 9,
    53305058                                               !< nwet = 9 should be handled outside this routine.
    5331 
    53325059    REAL(wp), INTENT(IN) ::  sai               !< surface area index
    53335060    REAL(wp), INTENT(IN) ::  ust               !< friction velocity (m/s)
    53345061    REAL(wp), INTENT(IN) ::  t                 !< temperature (C)
    5335 
    53365062    REAL(wp), INTENT(OUT) ::  gsoil_eff        !< effective soil conductance (m/s)
    5337 
    5338     !
    5339     !-- local variables:
     5063!
     5064!-- local variables:
    53405065    REAL(wp) ::  rinc                          !< in canopy resistance  (s/m)
    53415066    REAL(wp) ::  rsoil_eff                     !< effective soil resistance (s/m)
    5342 
    5343 
    5344     !
    5345     !-- Soil resistance (numbers matched with lu_classes and component numbers)
     5067!
     5068!-- Soil resistance (numbers matched with lu_classes and component numbers)
    53465069    !     grs    ara    crp    cnf    dec    wat    urb   oth    des    ice    sav    trf    wai    med    sem
    53475070    REAL(wp), PARAMETER ::  rsoil(nlu_dep,ncmp) = reshape( (/ &
     
    53575080         -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999. /),&  !< H2O2   
    53585081         (/nlu_dep,ncmp/) )
    5359 
    5360     !
    5361     !-- For                                          o3    so2   no2     no    nh3     co     no3    hno3   n2o5   h2o2
     5082!
     5083!-- For                                          o3    so2   no2     no    nh3     co     no3    hno3   n2o5   h2o2
    53625084    REAL(wp), PARAMETER ::  rsoil_wet(ncmp)    = (/2000., 10. , 2000., -999., 10.  , -999., -999., -999., -999., -999./)
    53635085    REAL(wp), PARAMETER ::  rsoil_frozen(ncmp) = (/2000., 500., 2000., -999., 1000., -999., -999., -999., -999., -999./)
    5364 
    5365 
    5366 
    5367     !
    5368     !-- Compute in canopy (in crop) resistance:
     5086!
     5087!-- Compute in canopy (in crop) resistance:
    53695088    CALL rc_rinc( lu, sai, ust, rinc )
    5370 
    5371     !
    5372     !-- Check for missing deposition path:
     5089!
     5090!-- Check for missing deposition path:
    53735091    IF ( missing(rinc) )  THEN
    53745092       rsoil_eff = -9999.0_wp
    53755093    ELSE
    5376        !
    5377        !-- Frozen soil (temperature below 0):
     5094!
     5095!--    Frozen soil (temperature below 0):
    53785096       IF ( t < 0.0_wp )  THEN
    53795097          IF ( missing( rsoil_frozen( icmp ) ) )  THEN
     
    53835101          ENDIF
    53845102       ELSE
    5385           !
    5386           !-- Non-frozen soil; dry:
     5103!
     5104!--      Non-frozen soil; dry:
    53875105          IF ( nwet == 0 )  THEN
    53885106             IF ( missing( rsoil( lu, icmp ) ) )  THEN
     
    53915109                rsoil_eff = rsoil( lu, icmp ) + rinc
    53925110             ENDIF
    5393 
    5394              !
    5395              !-- Non-frozen soil; wet:
     5111!
     5112!--       Non-frozen soil; wet:
    53965113          ELSEIF ( nwet == 1 )  THEN
    53975114             IF ( missing( rsoil_wet( icmp ) ) )  THEN
     
    54065123       ENDIF
    54075124    ENDIF
    5408 
    5409     !
    5410     !-- Compute conductance:
     5125!
     5126!-- Compute conductance:
    54115127    IF ( rsoil_eff > 0.0_wp )  THEN
    54125128       gsoil_eff = 1.0_wp / rsoil_eff
     
    54265142    IMPLICIT NONE
    54275143
    5428     !
    5429     !-- Input/output variables:
     5144!
     5145!-- Input/output variables:
    54305146    INTEGER(iwp), INTENT(IN) ::  lu          !< land use class, lu = 1, ..., nlu
    54315147
     
    54345150
    54355151    REAL(wp), INTENT(OUT) ::  rinc           !< in canopy resistance (s/m)
    5436 
    5437     !
    5438     !-- b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 if not applicable)
    5439     !-- h = vegetation height (m)                     gra  ara crop con dec wat   urb   oth   des   ice   sav   trf  wai  med semi
     5152!
     5153!-- b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 if not applicable)
     5154!-- h = vegetation height (m)                     gra  ara crop con dec wat   urb   oth   des   ice   sav   trf  wai  med semi
    54405155    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  b = (/ -999, 14, 14, 14, 14, -999, -999, -999, -999, -999, -999, 14, -999,  &
    54415156         14, 14 /)
    54425157    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  h = (/ -999, 1,  1,  20, 20, -999, -999, -999, -999, -999, -999, 20, -999,  &
    54435158         1 ,  1 /)
    5444 
    5445     !
    5446     !-- Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0:
     5159!
     5160!-- Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0:
    54475161    IF ( b(lu) > 0.0_wp )  THEN
    5448        !
    5449        !-- Check for u* > 0 (otherwise denominator = 0):
     5162!       !
     5163!--    Check for u* > 0 (otherwise denominator = 0):
    54505164       IF ( ust > 0.0_wp )  THEN
    54515165          rinc = b(lu) * h(lu) * sai/ust
     
    54645178
    54655179
    5466 
    54675180 !-------------------------------------------------------------------
    54685181 !> rc_rctot: compute total canopy (or surface) resistance Rc
     
    54725185    IMPLICIT NONE
    54735186
    5474     !
    5475     !-- Input/output variables:
     5187!
     5188!-- Input/output variables:
    54765189    REAL(wp), INTENT(IN) ::  gstom         !< stomatal conductance (s/m)
    54775190    REAL(wp), INTENT(IN) ::  gsoil_eff     !< effective soil conductance (s/m)
     
    54805193    REAL(wp), INTENT(OUT) ::  gc_tot       !< total canopy conductance (m/s)
    54815194    REAL(wp), INTENT(OUT) ::  rc_tot       !< total canopy resistance Rc (s/m)
    5482 
    5483     !
    5484     !-- Total conductance:
     5195!
     5196!-- Total conductance:
    54855197    gc_tot = gstom + gsoil_eff + gw
    5486 
    5487     !
    5488     !-- Total resistance (note: gw can be negative, but no total emission allowed here):
     5198!
     5199!-- Total resistance (note: gw can be negative, but no total emission allowed here):
    54895200    IF ( gc_tot <= 0.0_wp .OR. gw < 0.0_wp )  THEN
    54905201       rc_tot = -9999.0_wp
     
    54965207
    54975208
    5498 
    54995209 !-------------------------------------------------------------------
    55005210 !> rc_comp_point_rc_eff: calculate the effective resistance Rc
    55015211 !> based on one or more compensation points
    55025212 !-------------------------------------------------------------------
    5503  !
    55045213 !> NH3rc (see depac v3.6 is based on Avero workshop Marc Sutton. p. 173.
    55055214 !> Sutton 1998 AE 473-480)
     
    55585267 !>
    55595268 ! -------------------------------------------------------------------------------------------
    5560 
    55615269! SUBROUTINE rc_comp_point_rc_eff( ccomp_tot, catm, ra, rb, rc_tot, rc_eff )
    55625270!
    55635271!    IMPLICIT NONE
    55645272!
    5565 !    !
    5566 !    !-- Input/output variables:
     5273!!-- Input/output variables:
    55675274!    REAL(wp), INTENT(IN) ::  ccomp_tot    !< total compensation point (weighed average of separate compensation points) (ug/m3)
    55685275!    REAL(wp), INTENT(IN) ::  catm         !< atmospheric concentration (ug/m3)
     
    55745281!
    55755282!    !
    5576 !    !-- Compute effective resistance:
     5283!!-- Compute effective resistance:
    55775284!    IF (  ccomp_tot == 0.0_wp )  THEN
    55785285!       !
    5579 !       !-- trace with no compensiation point ( or compensation point equal to zero)
     5286!!--    trace with no compensiation point ( or compensation point equal to zero)
    55805287!       rc_eff = rc_tot
    55815288!
    55825289!    ELSE IF ( ccomp_tot > 0.0_wp .AND. ( abs( catm - ccomp_tot ) < 1.e-8 ) )  THEN
    55835290!       !
    5584 !       !--surface concentration (almost) equal to atmospheric concentration
    5585 !       !-- no exchange between surface and atmosphere, infinite RC --> vd=0
     5291!!--   surface concentration (almost) equal to atmospheric concentration
     5292!!--    no exchange between surface and atmosphere, infinite RC --> vd=0
    55865293!       rc_eff = 9999999999.0_wp
    55875294!
    55885295!    ELSE IF ( ccomp_tot > 0.0_wp )  THEN
    55895296!       !
    5590 !       !-- compensation point available, calculate effective resistance
     5297!!--    compensation point available, calculate effective resistance
    55915298!       rc_eff = ( ( ra + rb ) * ccomp_tot + rc_tot * catm ) / ( catm - ccomp_tot )
    55925299!
     
    56025309
    56035310
    5604 
    56055311 !-------------------------------------------------------------------
    56065312 !> missing: check for data that correspond with a missing deposition path
    56075313 !>          this data is represented by -999
    56085314 !-------------------------------------------------------------------
    5609 
    56105315 LOGICAL function missing( x )
    56115316
    56125317    REAL(wp), INTENT(IN) ::  x
    56135318
    5614     !
    5615     !-- bandwidth for checking (in)equalities of floats
     5319!
     5320!-- bandwidth for checking (in)equalities of floats
    56165321    REAL(wp), PARAMETER :: eps = 1.0e-5
    56175322
     
    56215326
    56225327
    5623 
    56245328 ELEMENTAL FUNCTION sedimentation_velocity( rhopart, partsize, slipcor, visc ) RESULT( vs )
    56255329
    5626 
    56275330    IMPLICIT NONE
    56285331
    5629     !
    5630     !-- in/out
     5332!
     5333!-- in/out
    56315334
    56325335    REAL(wp), INTENT(IN) ::  rhopart                 !< particle density (kg/m3)
     
    56365339
    56375340    REAL(wp) ::  vs
    5638 
    5639     !
    5640     !-- acceleration of gravity:
     5341!
     5342!-- acceleration of gravity:
    56415343    REAL(wp), PARAMETER         ::  grav = 9.80665   !< m/s2
    56425344
    5643 
    5644 
    5645     !
    5646     !-- sedimentation velocity
     5345!-- sedimentation velocity
    56475346    vs = rhopart * ( partsize**2.0_wp ) * grav * slipcor / ( 18.0_wp * visc )
    56485347
     
    56565355      luc, ftop_lu, ustar_lu )
    56575356
    5658 
    56595357    IMPLICIT NONE 
    56605358
    5661     ! -- in/out
     5359!
     5360!-- in/out
    56625361
    56635362    INTEGER(iwp), INTENT(IN) ::  nwet        !< 1=rain, 9=snowcover
     
    56755374    REAL(wp), INTENT(OUT) ::  vd             !< deposition velocity (m/s)
    56765375    REAL(wp), INTENT(OUT) ::  Rs             !< sedimentaion resistance (s/m)
    5677 
    5678 
    5679     !
    5680     !-- constants
     5376!
     5377!-- constants
    56815378
    56825379    REAL(wp), PARAMETER ::  grav     = 9.80665             !< acceleration of gravity (m/s2)
     
    56935390    REAL(wp), PARAMETER ::A_lu(nlu_dep) = &   
    56945391         (/3.0,  3.0,   2.0,  2.0,  7.0,  -99.,   10.0, 3.0, -99., -99.,  3.0, 7.0, -99., 2.0, -99./)
    5695     !
    5696     !--   grass  arabl crops conif decid  water  urba  othr  desr  ice   sav  trf   wai  med   sem     
    5697 
    5698     !
    5699     !-- local
    5700 
     5392!
     5393!--   grass  arabl crops conif decid  water  urba  othr  desr  ice   sav  trf   wai  med   sem     
     5394!
     5395!-- local
    57015396    REAL(wp) ::  kinvisc
    57025397    REAL(wp) ::  diff_part
     
    57075402    REAL(wp) ::  Einterc
    57085403    REAL(wp) ::  Reffic
    5709 
    5710 
    5711     !
    5712     !-- kinetic viscosity & diffusivity
     5404!
     5405!-- kinetic viscosity & diffusivity
    57135406    kinvisc = viscos1 / dens1    !< only needed at surface
    57145407
    57155408    diff_part = kb * tsurf * slipcor / ( 3 * pi * viscos1 * partsize )
    5716 
    5717     !
    5718     !-- Schmidt number
     5409!
     5410!-- Schmidt number
    57195411    schmidt = kinvisc / diff_part
    5720 
    5721     !
    5722     !-- calculate collection efficiencie E
     5412!
     5413!-- calculate collection efficiencie E
    57235414    Ebrown = Schmidt**( -gamma_lu(luc) )    !< Brownian diffusion
    5724     !
    5725     !-- determine Stokes number, interception efficiency
    5726     !-- and sticking efficiency R (1 = no rebound)
     5415!
     5416!-- determine Stokes number, interception efficiency
     5417!-- and sticking efficiency R (1 = no rebound)
    57275418    IF ( luc == ilu_ice .OR. nwet==9 .OR. luc == ilu_water_sea .OR. luc == ilu_water_inland )  THEN
    57285419       stokes = vs1 * ustar_lu**2 / ( grav * kinvisc )
     
    57385429       Reffic = exp( -Stokes**0.5_wp )
    57395430    END IF
    5740     !
    5741     !-- when surface is wet all particles do not rebound:
     5431!
     5432!-- when surface is wet all particles do not rebound:
    57425433    IF ( nwet==1 )  Reffic = 1.0_wp
    5743     !
    5744     !-- determine impaction efficiency:
     5434!
     5435!-- determine impaction efficiency:
    57455436    Eimpac = ( stokes / ( alfa_lu(luc) + stokes ) )**beta
    5746 
    5747     !
    5748     !-- sedimentation resistance:
     5437!
     5438!-- sedimentation resistance:
    57495439    Rs = 1.0_wp / ( epsilon0 * ustar_lu * ( Ebrown + Eimpac + Einterc ) * Reffic )
    57505440
    5751     !
    5752     !-- deposition velocity according to Seinfeld and Pandis (2006; eq 19.7):
    5753     !
    5754     !            1
    5755     !    vd = ------------------ + vs
    5756     !         Ra + Rs + Ra*Rs*vs
    5757     !
    5758     !-- where: Rs = Rb (in Seinfeld and Pandis, 2006)
    5759     !
    5760     !
     5441!-- deposition velocity according to Seinfeld and Pandis (2006; eq 19.7):
     5442!-- 
     5443!--              1
     5444!--      vd = ------------------ + vs
     5445!--           Ra + Rs + Ra*Rs*vs
     5446!-- 
     5447!-- where: Rs = Rb (in Seinfeld and Pandis, 2006)
    57615448
    57625449    vd = 1.0_wp / ( ftop_lu + Rs + ftop_lu * Rs * vs1) + vs1
     
    57725459 SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffc, Rb )   
    57735460
    5774    
     5461
    57755462    IMPLICIT NONE
    57765463
    5777     !-- in/out
     5464!
     5465!-- in/out
    57785466
    57795467    LOGICAL , INTENT(IN) ::  is_water
     
    57845472
    57855473    REAL(wp), INTENT(OUT) ::  Rb                  !< boundary layer resistance
    5786 
    5787     !
    5788     !-- const
     5474!
     5475!-- const
    57895476
    57905477    REAL(wp), PARAMETER ::  thk = 0.19e-4         !< thermal diffusivity of dry air 20 C
    57915478    REAL(wp), PARAMETER ::  kappa_stab = 0.35     !< von Karman constant
    5792 
    57935479!
    57945480!-- Next line is to avoid compiler warning about unused variable
    57955481    IF ( is_water  .OR.  ( z0h + kappa_stab ) > 0.0_wp )  CONTINUE
    5796 
    5797     !
    5798     !-- Use Simpson et al. (2003)
    5799     !-- @TODO: Check Rb over water calculation, until then leave commented lines
    5800     !IF ( is_water )  THEN
    5801     ! org: Rb = 1.0_wp / (kappa_stab*MAX(0.01_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.01_wp,ustar))
    5802     !      Rb = 1.0_wp / (kappa_stab*MAX(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.1_wp,ustar))
    5803     !ELSE
     5482!
     5483!-- Use Simpson et al. (2003)
     5484!-- @TODO: Check Rb over water calculation, until then leave commented lines
     5485!--  IF ( is_water )  THEN
     5486!--   org: Rb = 1.0_wp / (kappa_stab*MAX(0.01_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.01_wp,ustar))
     5487!--        Rb = 1.0_wp / (kappa_stab*MAX(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.1_wp,ustar))
     5488!--  ELSE
    58045489    Rb = 5.0_wp / MAX( 0.01_wp, ustar ) * ( thk / diffc )**0.67_wp
    5805     !END IF
     5490!--  END IF
    58065491
    58075492 END SUBROUTINE get_rb_cell
     
    58355520    IMPLICIT NONE
    58365521
    5837     !
    5838     !-- in/out
     5522!
     5523!-- in/out
    58395524
    58405525    REAL(wp), INTENT(IN) ::  q                      !< specific humidity [(kg water)/(kg air)]
     
    58425527
    58435528    REAL(wp) ::  p_w                                !< water vapor partial pressure [Pa]
    5844 
    5845     !
    5846     !-- const
     5529!
     5530!-- const
    58475531
    58485532    REAL(wp), PARAMETER  ::  eps = xm_h2o / xm_air  !< mole mass ratio ~ 0.622
    5849 
    5850     !
    5851     !-- partial pressure of water vapor:
     5533!
     5534!-- partial pressure of water vapor:
    58525535    p_w = p * q / eps
    58535536
    58545537 END function watervaporpartialpressure
    5855 
    58565538
    58575539
     
    58805562    IMPLICIT NONE
    58815563
    5882     !
    5883     !-- in/out
     5564!
     5565!-- in/out
    58845566
    58855567    REAL(wp), INTENT(IN) ::  t            !< temperature [K]
    58865568
    58875569    REAL(wp) ::  e_sat_w                  !< saturation vapor pressure  [Pa]
    5888 
    5889     !
    5890     !-- const
     5570!
     5571!-- const
    58915572    REAL(wp), PARAMETER ::  p0 = 611.2   !< base pressure [Pa]
    5892 
    5893     !
    5894     !-- saturation vapor pressure:
     5573!
     5574!-- saturation vapor pressure:
    58955575    e_sat_w = p0 * exp( 17.67_wp * ( t - 273.16_wp ) / ( t - 29.66_wp ) )    !< [Pa]
    58965576
    58975577 END FUNCTION saturationvaporpressure
    5898 
    58995578
    59005579
     
    59125591    IMPLICIT NONE
    59135592
    5914     !
    5915     !-- in/out
     5593!
     5594!-- in/out
    59165595
    59175596    REAL(wp), INTENT(IN) ::  q    !< specific humidity [(kg water)/(kg air)]
     
    59205599
    59215600    REAL(wp) ::  rh               !< relative humidity [%]
    5922 
    5923     !
    5924     !-- relative humidity:
     5601!
     5602!-- relative humidity:
    59255603    rh = watervaporpartialpressure( q, p ) / saturationvaporpressure( t ) * 100.0_wp
    59265604
    59275605 END FUNCTION relativehumidity_from_specifichumidity
    5928 
    59295606
    59305607     
Note: See TracChangeset for help on using the changeset viewer.