Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 3458 for palm/trunk


Ignore:
Timestamp:
Oct 30, 2018 2:51:23 PM (6 years ago)
Author:
kanani
Message:

Reintegrated fixes/changes from branch chemistry

Location:
palm/trunk
Files:
7 added
34 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/Makefile

    r3448 r3458  
    2525# -----------------
    2626# $Id$
     27# from chemistry branch r3443, banzhafs, Russo, forkel, basit:
     28# radiation_model_mod added to chemistry_model_mod
     29# Added some missing dependencies and replaced blanks with tabs
     30# Added chemistry emission module
     31# chemistry_model_mod added to flow_statistics
     32#
     33# 3448 2018-10-29 18:14:31Z kanani
    2734# Adjustment of biometeorology dependencies
    2835#
    … …  
    856863        modules.o \
    857864        netcdf_data_input_mod.o \
     865        radiation_model_mod.o \
    858866        surface_mod.o \
    859867        user_module.o

  • palm/trunk/SOURCE/chem_emissions_mod.f90

    r3373 r3458  
    2727! -----------------
    2828! $Id$
     29! from chemistry branch r3443, banzhafs, Russo:
     30! Additional correction for index of input file of pre-processed mode
     31! Removed atomic and molecular weights as now available in chem_modules and
     32! added par_emis_time_factor (formerly in netcdf_data_input_mod)
     33! Added correction for index of input file of pre-processed mode
     34! Added a condition for default mode necessary for information read from ncdf file
     35! in pre-processed and default mode
     36! Correction of multiplying factor necessary for scaling emission values in time
     37! Introduced correction for scaling units in the case of DEFAULT emission mode
     38!
     39! 3373 2018-10-18 15:25:56Z kanani
    2940! Fix wrong location of __netcdf directive
    3041!
    … …  
    4051! Authors:
    4152! --------
    42 ! Emmanuele Russo Fu-Berlin
    43 ! Sabine Banzhaf  FU-Berlin
    44 ! Martijn Schaap  FU-Berlin, TNO Utrecht
     53! @author Emmanuele Russo (Fu-Berlin)
     54! @author Sabine Banzhaf  (FU-Berlin)
     55! @author Martijn Schaap  (FU-Berlin, TNO Utrecht)
    4556!
    4657! Description:
    … …  
    4960!>
    5061!> @todo Check_parameters routine should be developed: for now it includes just one condition
    51 !> @todo Add option for capital or small letters in the matching routine
    5262!> @todo Use of Restart files not contempled at the moment
     63!> @todo revise indices of files read from the netcdf: preproc_emission_data and expert_emission_data
     64!> @todo for now emission data may be passed on a singular vertical level: need to be more flexible
    5365!> @note <Enter notes on the module>
    5466!> @bug  <Enter known bugs here>
    … …  
    6173
    6274    USE control_parameters,                                                    &
    63        ONLY:  initializing_actions,end_time, message_string,                   &
     75       ONLY:  initializing_actions, end_time, message_string,                  &
    6476              intermediate_timestep_count, dt_3d
    6577 
    … …  
    6779
    6880    USE kinds
     81
     82#if defined ( __netcdf )
    6983
    7084    USE netcdf_data_input_mod,                                                  &
    7185       ONLY: chem_emis_att_type, chem_emis_val_type
    7286
    73 #if defined ( __netcdf )
    7487    USE NETCDF
     88
    7589#endif
    7690
    7791    USE date_and_time_mod,                                                      &
    78        ONLY: time_default_indices, month_of_year, day_of_month, hour_of_day,    &
     92       ONLY: time_default_indices, time_preprocessed_indices,                  &
     93             month_of_year, day_of_month, hour_of_day,                          &
    7994             index_mm, index_dd, index_hh
    8095
    … …  
    92107   
    93108    CHARACTER (LEN=80)                               :: filename_emis                   !> Variable for the name of the netcdf input file
    94     CHARACTER (LEN=80), ALLOCATABLE, DIMENSION(:)    :: string_values                   !> Output of string variables read from netcdf input file
    95 
    96     INTEGER(iwp)                                     :: n_dims                          !> Number of dimensions of input variable netcdf file
     109
    97110    INTEGER(iwp)                                     :: i                               !> index 1st selected dimension (some dims are not spatial)
    98111    INTEGER(iwp)                                     :: j                               !> index 2nd selected dimension
    … …  
    103116    INTEGER(iwp)                                     :: z_start                         !> Index to start read variable from netcdf in additional dims
    104117    INTEGER(iwp)                                     :: z_end                           !> Index to end read variable from netcdf in additional     dims
    105     INTEGER(iwp), ALLOCATABLE, DIMENSION(:)          :: id_dims                         !> id of dimension of selected variable netcdf input
    106     INTEGER(iwp), ALLOCATABLE, DIMENSION(:)          :: dum_var                         !> value of variable read from netcdf input
    107     INTEGER(iwp)                                     :: id_ncfile                       !> id netcdf file
    108     INTEGER(iwp)                                     :: errno                           !> error number NF90_???? function
    109     INTEGER(iwp)                                     :: id_var                          !< variable id
    110118    INTEGER(iwp)                                     :: dt_emis                         !> Time Step Emissions
    111119    INTEGER(iwp)                                     :: len_index                       !> length of index (used for several indices)
    … …  
    124132    !> Dobson units:
    125133    REAL, PARAMETER        ::  Dobs = 2.68668e16    ! (mlc/cm2) / DU
    126    
    127     ! molar weights of components
    128  
    129     !> atomic weights: 
    130     REAL, PARAMETER        ::  xm_H     =    1.00790e-3               ! kg/mol
    131     REAL, PARAMETER        ::  xm_N     =   14.00670e-3               ! kg/mol
    132     REAL, PARAMETER        ::  xm_C     =   12.01115e-3               ! kg/mol
    133     REAL, PARAMETER        ::  xm_S     =   32.06400e-3               ! kg/mol
    134     REAL, PARAMETER        ::  xm_O     =   15.99940e-3               ! kg/mol
    135     REAL, PARAMETER        ::  xm_F     =   18.99840e-3               ! kg/mol
    136     REAL, PARAMETER        ::  xm_Na    =   22.98977e-3               ! kg/mol
    137     REAL, PARAMETER        ::  xm_Cl    =   35.45300e-3               ! kg/mol
    138     REAL, PARAMETER        ::  xm_Rn222 =  222.00000e-3               ! kg/mol
    139     REAL, PARAMETER        ::  xm_Pb210 =  210.00000e-3               ! kg/mol
    140     REAL, PARAMETER        ::  xm_Ca    =   40.07800e-3               ! kg/mol
    141     REAL, PARAMETER        ::  xm_K     =   39.09800e-3               ! kg/mol
    142     REAL, PARAMETER        ::  xm_Mg    =   24.30500e-3               ! kg/mol
    143     REAL, PARAMETER        ::  xm_Pb    =  207.20000e-3               ! kg/mol
    144     REAL, PARAMETER        ::  xm_Cd    =  112.41000e-3               ! kg/mol
    145     REAL, PARAMETER        ::  xm_Rh    =  102.90550e-3               ! kg/mol
    146  
    147     !> molecule weights:
    148     REAL, PARAMETER        ::  xm_h2o   = xm_H * 2 + xm_O             ! kg/mol
    149     REAL, PARAMETER        ::  xm_o3    = xm_O * 3                    ! kg/mol
    150     REAL, PARAMETER        ::  xm_N2O5  = xm_N * 2 + xm_O * 5         ! kg/mol
    151     REAL, PARAMETER        ::  xm_HNO3  = xm_H + xm_N + xm_O * 3      ! kg/mol
    152     REAL, PARAMETER        ::  xm_NH4   = xm_N + xm_H * 4             ! kg/mol
    153     REAL, PARAMETER        ::  xm_SO4   = xm_S + xm_O * 4             ! kg/mol
    154     REAL, PARAMETER        ::  xm_NO3   = xm_N + xm_O * 3             ! kg/mol
    155     REAL, PARAMETER        ::  xm_CO2   = xm_C + xm_O * 2             ! kg/mol
    156    
    157     !> mass of air
    158     REAL, PARAMETER        ::  xm_air   =  28.964e-3                  ! kg/mol
    159     REAL, PARAMETER        ::  xmair    =  28.94                      ! g/mol; old name!
    160134
    161135    !> sesalt composition:
    … …  
    263237!-- Public Variables
    264238
    265     PUBLIC id_ncfile,errno,id_dims,n_dims,dum_var, con_factor,   &
    266            len_index,len_index_voc,len_index_pm, string_values
     239    PUBLIC con_factor, len_index,len_index_voc,len_index_pm
     240
    267241 CONTAINS
    268242
    … …  
    275249
    276250
    277 !TBD: Where Should we put the call to chem_emissions_check_parameters? In chem_init or in check_parameters?
     251    !TBD: Where Should we put the call to chem_emissions_check_parameters? In chem_init or in check_parameters?
    278252
    279253    IMPLICIT NONE
    … …  
    548522                IF (len_index_voc>0) THEN
    549523                   ALLOCATE(match_spec_voc_model(len_index_voc))   !> contains indices of the VOC model species
    550                    ALLOCATE(match_spec_voc_input(len_index_voc))   !> In input there is only one value for VOCs in the DEFAULT mode. This array contains the indices of the different values of VOC compositions of the input variable VOC_composition
     524                   ALLOCATE(match_spec_voc_input(len_index_voc))   !> In input there is only one value for VOCs in the DEFAULT mode.
     525                                                                   !  This array contains the indices of the different values of VOC compositions of the input variable VOC_composition
    551526                ENDIF
    552527
    … …  
    661636                                 ' DO NOT MATCH'                                //          &
    662637                                 ' model chemical species'                      //          &
    663                                  ' Chemistry Emissions routine is not called'         !TBD: IMPORTANT:Add a condition in chem init that sets the emission_output_required to false when len_index==0 so that the chem_emissions modules are not used
     638                                 ' Chemistry Emissions routine is not called'         
    664639                CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 )
    665640
    … …  
    724699                DO  ispec = 1 , len_index
    725700
    726                    IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND.                         &
    727                         emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN
    728 
    729                       message_string = 'PARAMETERIZED emissions mode selected:'            //           &
     701                   IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND. emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN
     702
     703                      message_string = 'PARAMETERIZED emissions mode selected:'            //          &
    730704                                       ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES'        //          &
    731705                                       ' but values of scaling factors for street types'    //          &
    … …  
    780754 SUBROUTINE chem_emissions_init(emt_att,emt,nspec_out)
    781755
    782 
     756#if defined( __netcdf )
    783757
    784758    USE surface_mod,                                                           &
    … …  
    799773 
    800774    INTEGER(iwp)                                                      :: ispec     !> Index to go through the emission chemical species
    801 #if defined( __netcdf )
     775
    802776
    803777!-- Actions for initial runs : TBD: needs to be updated
    … …  
    854828     SELECT CASE(TRIM(mode_emis))   !TBD: Add the option for CApital or small letters
    855829
     830
     831        !> PRE-PROCESSED case
     832        CASE ("PRE-PROCESSED")
     833
     834           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out)) 
     835
     836           CALL location_message( 'emis_distribution array allocated in PRE-PROCESSED mode', .FALSE. )
     837 
     838           !> Calculate the values of the emissions at the first time step
     839           CALL chem_emissions_setup(emt_att,emt,nspec_out)
     840
    856841        !> Default case
    857842        CASE ("DEFAULT")
    … …  
    864849           CALL chem_emissions_setup(emt_att,emt,nspec_out)
    865850
    866 
    867         !> PRE-PROCESSED case
    868         CASE ("PRE-PROCESSED")
    869 
    870            IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out)) 
    871 
    872            CALL location_message( 'emis_distribution array allocated in PRE-PROCESSED mode', .FALSE. )
    873  
    874            !> Calculate the values of the emissions at the first time step
    875            CALL chem_emissions_setup(emt_att,emt,nspec_out)
    876 
     851        !> PARAMETERIZED case
    877852        CASE ("PARAMETERIZED")
    878853
    879854           CALL location_message( 'emis_distribution array allocated in PARAMETERIZED mode', .FALSE. )
    880855
    881            ! For now for PAR and DEF values only, first vertical level of emis_distribution is allocated, while for EXP all.
     856           ! For now for PAR and DEF values only, first vertical level of emis_distribution is allocated, while for PRE-PROCESSED all.
    882857           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(1,0:ny,0:nx,nspec_out))
    883858
    … …  
    919894 IMPLICIT NONE
    920895
    921 
     896#if defined( __netcdf )
    922897 
    923898    !--- IN/OUT
    … …  
    949924    REAL(wp),ALLOCATABLE, DIMENSION(:,:)                              ::  emis
    950925
     926    REAL(wp), DIMENSION(24)                                           :: par_emis_time_factor      !< time factors
     927                                                                                      !  for the parameterized mode: these are fixed for each hour
     928                                                                                      !  of a single day.
    951929    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr)                    ::  conv_to_ratio !> factor used for converting input
    952930                                                                                        !  to adimensional concentration ratio
    … …  
    964942    ! --- const -------------------------------
    965943    !-CONVERSION FACTORS: TIME
     944    ! number of sec per hour = 3600   
     945    REAL, PARAMETER   ::  s_per_hour = 3600.0  !  (s)/(hour)
     946    ! number of sec per day = 86400   
     947    REAL, PARAMETER   ::  s_per_day = 86400.0  !  (s)/(day)
    966948    ! number of hours in a year of 365 days:
    967     REAL, PARAMETER  ::  hour_per_year = 8760.0 !> TBD: What about leap years?
     949    REAL, PARAMETER   ::  hour_per_year = 8760.0 !> TBD: What about leap years?
     950    ! number of hours in a day:
     951    REAL, PARAMETER   ::  hour_per_day = 24.0
     952
    968953    ! conversion from hours to seconds (s/hour) = 1/3600.0 ~ 0.2777778   
    969     REAL, PARAMETER   ::  hour_to_s = 0.0002777778  !  (hour)/(s)
     954    REAL, PARAMETER   ::  hour_to_s = 1/s_per_hour  !  (hour)/(s)
    970955    ! conversion from day to seconds (s/day) = 1/86400.0 ~ 1.157407e-05   
    971     REAL, PARAMETER   ::  day_to_s = 1.157407e-05   !  (day)/(s)
     956    REAL, PARAMETER   ::  day_to_s = 1/s_per_day   !  (day)/(s)
    972957    ! conversion from year to sec (s/year) = 1/31536000.0 ~ 3.170979e-08   
    973     REAL, PARAMETER   ::  year_to_s = 3.170979e-08  !  (year)/(s)
     958    REAL, PARAMETER   ::  year_to_s = 1/(s_per_hour*hour_per_year)  !  (year)/(s)
    974959
    975960    !-CONVERSION FACTORS: WEIGHT
    … …  
    984969    REAL(wp), PARAMETER   ::  ratio2ppm  = 1.0e06_wp 
    985970    !------------------------------------------------------   
    986 #if defined( __netcdf )
     971
    987972    IF ( emission_output_required ) THEN
    988973
    989     !>  Set emis_dt  !TBD: this is the same as dt_chem. We should consider the fact that dt_emis should be the timestep of input emissions or better defined, the timestep at which the emission routines are calle: for now one hour. It should be made changeable.
     974    !>  Set emis_dt  !TBD: this is the same as dt_chem. We should consider the fact that dt_emis should be the timestep of input emissions or better defined, the timestep at which the emission routines are called: for now one hour. It should be made changeable.
    990975
    991976       IF ( call_chem_at_all_substeps )  THEN
    … …  
    10931078 
    10941079          message_string = 'No Units conversion required for units of chemistry emissions' // &
    1095                            ' of the PARAMETERIZED mode: units have to be provided in mole/m**2/s '
     1080                           ' of the PARAMETERIZED mode: units have to be provided in'     //  &
     1081                           ' micromole/m**2/day for GASES and'                            //  &
     1082                           ' kg/m**2/day for PMs'                     
    10961083          CALL message( 'chem_emissions_setup', 'CM0447', 0, 0, 0, 6, 0 )
    10971084
    … …  
    11101097             !  V/N=RT/P
    11111098
    1112              !>    m**3/Nmole             (J/mol)*K^-1           K                      Pa         
     1099             !>    m**3/Nmole              (J/mol)*K^-1           K                      Pa         
    11131100             conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) ) 
    11141101          ENDDO
    … …  
    11251112    !> PRE-PROCESSED MODE
    11261113       IF (TRIM(mode_emis)=="PRE-PROCESSED") THEN
     1114
     1115          !> Update time indices
     1116          CALL time_preprocessed_indices(index_hh)
    11271117
    11281118          CALL location_message( 'PRE-PROCESSED MODE: No time-factor specification required', .FALSE. )
    … …  
    11941184                                 ' 24 values for every day of the year ', .FALSE. )
    11951185       
     1186          !Assign Constant Values of time factors, diurnal time profile for traffic sector:
     1187          par_emis_time_factor( : ) = (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039, 0.056, 0.053, 0.051, 0.051, 0.052, 0.055, &
     1188                                                 0.059, 0.061, 0.064, 0.067, 0.069, 0.069, 0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /)
     1189         
    11961190          !> in this case allocate time factor each hour in a day
    11971191          IF (.NOT. ALLOCATED(time_factor)) ALLOCATE(time_factor(1))
    … …  
    12001194          index_hh=hour_of_day
    12011195
    1202           time_factor(1)=emt_att%par_emis_time_factor(index_hh)
     1196          time_factor(1) = par_emis_time_factor(index_hh)
    12031197
    12041198       ENDIF
    … …  
    12121206          DO ispec=1,nspec_out
    12131207
    1214              !> Values are still micromoles/(m**2*s). Units are not given by the user for this mode and are always micromoles (or micrograms for PMs)
    1215              emis_distribution(1,nys:nyn,nxl:nxr,ispec)=surface_csflux(match_spec_input(ispec))*time_factor(1)
     1208             !> Values are still micromoles/(m**2*s). Units in this case are always micromoles/m**2*day (or kilograms/m**2*day for PMs)
     1209             emis_distribution(1,nys:nyn,nxl:nxr,ispec)=surface_csflux(match_spec_input(ispec))*time_factor(1)*hour_to_s
    12161210
    12171211          ENDDO
    … …  
    12231217          DO ispec=1,nspec_out !> nspec_out represents the number of species in common between
    12241218                               !  the emission input data and the chemistry mechanism used
    1225 
    1226              emis_distribution(nzb:nzt+1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))%                                  &
    1227                                                                       preproc_emission_data(index_hh,nzb:nzt+1,nys:nyn,nxl:nxr)* &
    1228                                                                    con_factor         
    1229 
     1219   
     1220             emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))%                               &
     1221                                                                   preproc_emission_data(index_hh,1,nys+1:nyn+1,nxl+1:nxr+1)* &
     1222                                                                      con_factor
     1223         
    12301224          ENDDO
    12311225
    … …  
    12491243!TBD: The consideration of dt_emis of the input data is still missing. Basically the emissions could be provided every 10, 30 minutes and not always at one hour. This should be eventually solved in the date_and_time mode routine.
    12501244
     1245             !> the time factors are 24 for each day. When multiplied by a daily value, they allow to have an hourly value. Then to convert it to seconds, we still have to divide this value by 3600.
     1246             !  So given any units, we convert them to seconds and finally multiply them by 24 ((value/sec)*(24*3600)=value/day ---- (value/day)*time_factor=value/hour ---(value/hour)/(3600)=value/sec )
     1247             !                                                                                 ((value/sec)*(24*3600)*time_factor)/3600=24*(value/sec)*time_factor                         
     1248
    12511249             !> NOX Compositions
    12521250                IF (TRIM(spc_names(match_spec_model(ispec)))=="NO") THEN
    1253                 !>             Kg/m2                       kg/m2*s                                                   
    1254                    delta_emis(nys:nyn,nxl:nxr) =                                                             &
    1255                          emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor
     1251                !>             Kg/m2*s                   kg/m2*s                                                   
     1252                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor*hour_per_day
    12561253                   
    1257                    emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
    1258                          emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1254                   emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    12591255
    12601256                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="NO2") THEN
    12611257   
    1262                    delta_emis(nys:nyn,nxl:nxr) =                                                             &
    1263                          emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,2)*con_factor
    1264 
    1265                    emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
    1266                          emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1258                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,2)*con_factor*hour_per_day
     1259
     1260                   emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    12671261 
    12681262             !> SOX Compositions
    12691263                     
    12701264                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO2") THEN
    1271                      
    1272                    delta_emis(nys:nyn,nxl:nxr) =                                                             &
    1273                          emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor
    1274 
    1275                    emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
    1276                          emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1265                   !>             Kg/m2*s                   kg/m2*s                                                                       
     1266                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor*hour_per_day
     1267
     1268                   emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    12771269
    12781270                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO4") THEN
    1279    
    1280                    delta_emis(nys:nyn,nxl:nxr) =                                                             &
    1281                          emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor
    1282 
    1283                    emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
    1284                          emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1271                   !>             Kg/m2*s                   kg/m2*s                                                                         
     1272                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor*hour_per_day
     1273
     1274                   emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    12851275 
    12861276
    … …  
    12921282                   DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,1))
    12931283
    1294                       delta_emis(nys:nyn,nxl:nxr) =                                                          &
    1295                             emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,1)*con_factor
     1284                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,1)*con_factor*hour_per_day
    12961285                                                                                         
    12971286
    1298                       emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                           &
    1299                             emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1287                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    13001288 
    13011289                   ENDDO
    … …  
    13071295                   DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,2))
    13081296
    1309                       delta_emis(nys:nyn,nxl:nxr) =                                                          &
    1310                             emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,2)*con_factor
     1297                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,2)*con_factor*hour_per_day
    13111298                                                                                         
    13121299
    1313                       emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                           &
    1314                             emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1300                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    13151301 
    13161302                   ENDDO
    … …  
    13221308                   DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,3)) 
    13231309                       
    1324                       delta_emis(nys:nyn,nxl:nxr) =                                                          &
    1325                             emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,3)*con_factor
     1310                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,3)*con_factor*hour_per_day
    13261311                                                                                                 
    13271312
    1328                       emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                           &
    1329                             emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1313                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    13301314
    13311315                   ENDDO
    … …  
    13391323                      IF (TRIM(spc_names(match_spec_model(ispec)))==TRIM(emt_att%voc_name(ivoc))) THEN   
    13401324
    1341                          delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*               &
    1342                                                        emt_att%voc_comp(icat,match_spec_voc_input(ivoc))*con_factor
    1343 
    1344                          emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                         &
    1345                                emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1325                         delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*                                    &
     1326                                                       emt_att%voc_comp(icat,match_spec_voc_input(ivoc))*con_factor*hour_per_day
     1327
     1328                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    13461329
    13471330                      ENDIF                       
    … …  
    13521335                ELSE
    13531336
    1354                    delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*con_factor
    1355  
    1356                    emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                               &
    1357                          emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1337                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*con_factor*hour_per_day
     1338 
     1339                   emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    13581340
    13591341                ENDIF  ! IF (spc_names==)
    … …  
    13961378                      DO  ispec=1,nspec_out
    13971379
    1398                          !> PMs are already in mass units:micrograms:they have to be converted to kilograms
    1399                          IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"         &
    1400                              .OR.  TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
     1380                         !> PMs are already in mass units: kilograms
     1381                         IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
    14011382                             .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
    14021383
    … …  
    14121393                         !> Other Species: inputs are micromoles: they have to be converted in moles
    14131394                         !                 ppm/s *m *kg/m^3               
    1414                             surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))*  &
     1395                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))*   &
    14151396                         !                                                    micromoles/(m^2*s)
    1416                                                                           emis_distribution(1,j,i,ispec) *             &
    1417                          !                                                    m**3/Nmole
    1418                                                                           conv_to_ratio(k,j,i)*                        &
     1397                                                                          emis_distribution(1,j,i,ispec) *              &
     1398                         !                                                    m^3/Nmole
     1399                                                                          conv_to_ratio(k,j,i)*                         &       
    14191400                         !                                                    kg/m^3
    14201401                                                                          rho_air(k)   
    … …  
    14261407
    14271408
    1428                    ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND.                                            &
    1429                             street_type_f%var(j,i) < main_street_id )  THEN
     1409                   ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND. street_type_f%var(j,i) < main_street_id )  THEN
    14301410
    14311411                   !> Cycle over already matched species
    … …  
    14331413
    14341414                         !> PMs are already in mass units: micrograms
    1435                          IF (     TRIM(spc_names(match_spec_model(ispec)))=="PM1"   &
    1436                              .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
     1415                         IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
    14371416                             .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
    14381417
    14391418                            !              kg/(m^2*s) *kg/m^3                               
    1440                             surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) *  &
     1419                            surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) *   &
    14411420                            !                                                       kg/(m^2*s)
    14421421                                                                                emis_distribution(1,j,i,ispec)*        &
    … …  
    14481427                         !>Other Species: inputs are micromoles
    14491428                         !                 ppm/s *m *kg/m^3               
    1450                             surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * &
     1429                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) *   &
    14511430                         !                                                    micromoles/(m^2*s)
    1452                                                                           emis_distribution(1,j,i,ispec) *             &
    1453                          !                                                    m**3/Nmole
    1454                                                                           conv_to_ratio(k,j,i)*                        &
     1431                                                                          emis_distribution(1,j,i,ispec) *              &
     1432                         !                                                    m^3/Nmole
     1433                                                                          conv_to_ratio(k,j,i)*                         &       
    14551434                         !                                                    kg/m^3
    14561435                                                                          rho_air(k)   
    … …  
    14861465                   j = surf_def_h(0)%j(m)
    14871466
    1488                    !> Distinguish between PMs (no needing conversion in ppms),
    1489                    !  VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and
    1490                    ! other Species (still expressed in Kg/(m**2*s) at this point)
    1491 
    1492                    !> PMs
    1493                    IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"                                   &
    1494                          .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"                           &
    1495                          .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
     1467                   IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
     1468
     1469
     1470                      !> Distinguish between PMs (no needing conversion in ppms),
     1471                      !  VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and
     1472                      ! other Species (still expressed in Kg/(m**2*s) at this point)
     1473
     1474                      !> PMs
     1475                      IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
     1476                             .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
    14961477                   
    1497                       !            kg/(m^2*s) *kg/m^3                         kg/(m^2*s)                 
    1498                       surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)*   &
    1499                       !                                                  kg/m^3
    1500                                                                        rho_air(k)   
    1501  
    1502  
    1503                    ELSE
    1504 
    1505                       !> VOCs
    1506                       IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN
    1507                          !          ( ppm/s) * m * kg/m^3                         mole/(m**2/s)   
    1508                          surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *      &
    1509                                                                          !    m**3/mole          ppm
    1510                                                                           conv_to_ratio(k,j,i)*ratio2ppm *      &
    1511                       !                                                    kg/m^3
    1512                                                                           rho_air(k)   
    1513 
    1514 
    1515                       !> OTHER SPECIES
     1478                         !            kg/(m^2*s) *kg/m^3                         kg/(m^2*s)                 
     1479                         surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)*   &
     1480                         !                                                  kg/m^3
     1481                                                                          rho_air(nzb)   
     1482 
     1483 
    15161484                      ELSE
    15171485
    1518                          !               ( ppm/s )*m  * kg/m^3                      kg/(m**2/s)                     
    1519                          surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *      &
    1520                                                                             !  mole/Kg       
    1521                                                                           (1/emt_att%xm(ispec))*                &
    1522                          !                                                    m**3/mole          ppm
    1523                                                                           conv_to_ratio(k,j,i)*ratio2ppm*       &
    1524                          !                                                  kg/m^3
    1525                                                                           rho_air(k)   
    1526  
     1486                         !> VOCs
     1487                         IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN
     1488                            !          ( ppm/s) * m * kg/m^3                         mole/(m^2/s)   
     1489                            surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *      &
     1490                                                                            !    m^3/mole          ppm
     1491                                                                             conv_to_ratio(nzb,j,i)*ratio2ppm *      &
     1492                         !                                                    kg/m^3
     1493                                                                             rho_air(nzb)   
     1494
     1495
     1496                         !> OTHER SPECIES
     1497                         ELSE
     1498
     1499                            !               ( ppm/s )*m  * kg/m^3                      kg/(m^2/s)                     
     1500                            surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *      &
     1501                                                                               !  mole/Kg       
     1502                                                                             (1/emt_att%xm(ispec))*                &
     1503                            !                                                    m^3/mole          ppm
     1504                                                                             conv_to_ratio(nzb,j,i)*ratio2ppm*       &
     1505                            !                                                  kg/m^3
     1506                                                                             rho_air(nzb)   
     1507 
     1508
     1509                         ENDIF
    15271510
    15281511                      ENDIF
    … …  
    15431526                   k = surf_lsm_h%k(m)
    15441527
    1545                    !> Distinguish between PMs (no needing conversion in ppms),
    1546                    !  VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and
    1547                    ! other Species (still expressed in Kg/(m**2*s) at this point)
    1548 
    1549                    !> PMs
    1550                    IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"                                            &
    1551                          .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"                                    &
    1552                          .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
    1553 
    1554                       !         kg/(m^2*s) * kg/m^3                           kg/(m^2*s)           
    1555                       surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *              &
    1556                       !                                                  kg/m^3
    1557                                                                        rho_air(k)   
    1558  
    1559                    ELSE
    1560 
    1561                       !> VOCs
    1562                       IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN
    1563                          !          ( ppm/s) * m * kg/m^3                        mole/(m**2/s)   
    1564                          surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *      &
    1565                                                                        !    m**3/mole          ppm
    1566                                                                        conv_to_ratio(k,j,i)*ratio2ppm*    &
    1567                       !                                                 kg/m^3
    1568                                                                        rho_air(k)   
    1569 
    1570 
    1571                       !> OTHER SPECIES
     1528                   IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
     1529
     1530                      !> Distinguish between PMs (no needing conversion in ppms),
     1531                      !  VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and
     1532                      ! other Species (still expressed in Kg/(m**2*s) at this point)
     1533
     1534                      !> PMs
     1535                      IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
     1536                             .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
     1537   
     1538                         !         kg/(m^2*s) * kg/m^3                           kg/(m^2*s)           
     1539                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *              &
     1540                         !                                                  kg/m^3
     1541                                                                          rho_air(k)   
     1542 
    15721543                      ELSE
    15731544
    1574                          !         ( ppm/s) * m * kg/m^3                        kg/(m**2/s)                     
    1575                          surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *               &
    1576                                                                             !  mole/Kg   
    1577                                                                       (1/emt_att%xm(ispec))*                          &
    1578                          !                                                m**3/mole           ppm
    1579                                                                       conv_to_ratio(k,j,i)*ratio2ppm*                 &
    1580                          !                                            kg/m^3
    1581                                                                       rho_air(k)   
     1545                         !> VOCs
     1546                         IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN
     1547                            !          ( ppm/s) * m * kg/m^3                        mole/(m^2/s)   
     1548                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *      &
     1549                                                                          !    m^3/mole          ppm
     1550                                                                          conv_to_ratio(k,j,i)*ratio2ppm*    &
     1551                         !                                                 kg/m^3
     1552                                                                          rho_air(k)   
     1553
     1554                         !> OTHER SPECIES
     1555                         ELSE
     1556   
     1557                            !         ( ppm/s) * m * kg/m^3                        kg/(m^2/s)                     
     1558                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *               &
     1559                                                                               !  mole/Kg   
     1560                                                                         (1/emt_att%xm(ispec))*                          &
     1561                            !                                                m^3/mole           ppm
     1562                                                                         conv_to_ratio(k,j,i)*ratio2ppm*                 &
     1563                            !                                            kg/m^3
     1564                                                                         rho_air(k)   
    15821565                                                     
     1566                         ENDIF
     1567
    15831568                      ENDIF
    15841569
    15851570                   ENDIF
     1571
    15861572                ENDDO
    15871573
    … …  
    15981584                   k = surf_usm_h%k(m)
    15991585
    1600 
    1601                    !> PMs
    1602                    IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"                                     &
    1603                          .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"                             &
    1604                          .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
    1605                    
    1606                       !          kg/(m^2*s) *kg/m^3                             kg/(m^2*s)                     
    1607                       surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)*        &
    1608                       !                                              kg/m^3
    1609                                                                     rho_air(k)   
    1610 
    1611  
    1612                    ELSE
    1613 
    1614                       !> VOCs
    1615                       IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN
    1616                          !          ( ppm/s) * m * kg/m^3                        mole/(m**2/s)   
    1617                          surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *   &
    1618                                                                        !    m**3/mole          ppm
    1619                                                                        conv_to_ratio(k,j,i)*ratio2ppm*    &
    1620                       !                                                kg/m^3
     1586                   IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
     1587
     1588                      !> PMs
     1589                      IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
     1590                             .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
     1591                     
     1592                         !          kg/(m^2*s) *kg/m^3                             kg/(m^2*s)                     
     1593                         surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)*        &
     1594                         !                                              kg/m^3
    16211595                                                                       rho_air(k)   
    16221596
    1623                       !> OTHER SPECIES
     1597 
    16241598                      ELSE
    16251599
    1626 
    1627                       !            ( ppm/s) * m * kg/m^3                        kg/(m**2/s)                     
    1628                          surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *      &
    1629                                                                           !  mole/Kg   
    1630                                                                       (1/emt_att%xm(ispec))*                 &
    1631                          !                                                m**3/mole           ppm
    1632                                                                       conv_to_ratio(k,j,i)*ratio2ppm*        &
    1633                          !                                            kg/m^3
    1634                                                                       rho_air(k)   
    1635 
     1600                         !> VOCs
     1601                         IF ( len_index_voc .GT. 0 .AND. emt_att%species_name(match_spec_input(ispec))=="VOC" ) THEN
     1602                            !          ( ppm/s) * m * kg/m^3                        mole/(m^2/s)   
     1603                            surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *   &
     1604                                                                          !    m^3/mole          ppm
     1605                                                                          conv_to_ratio(k,j,i)*ratio2ppm*    &
     1606                         !                                                kg/m^3
     1607                                                                          rho_air(k)   
     1608
     1609                         !> OTHER SPECIES
     1610                         ELSE
     1611
     1612
     1613                         !            ( ppm/s) * m * kg/m^3                        kg/(m^2/s)                     
     1614                            surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *      &
     1615                                                                             !  mole/Kg   
     1616                                                                         (1/emt_att%xm(ispec))*                 &
     1617                            !                                                m^3/mole           ppm
     1618                                                                         conv_to_ratio(k,j,i)*ratio2ppm*        &
     1619                            !                                            kg/m^3
     1620                                                                         rho_air(k)   
     1621
     1622
     1623                         ENDIF
    16361624
    16371625                      ENDIF

  • palm/trunk/SOURCE/chem_modules.f90

    r3298 r3458  
    2727! -----------------
    2828! $Id$
     29! from chemistry branch r3443:
     30! ??
     31!
     32! 3298 2018-10-02 12:21:11Z kanani
    2933! - Minor formatting
    3034! - Introduced Variables for Chemistry emissions Module (Russo)
    … …  
    4347! @author Farah Kanani-Suehring
    4448! @author Basit Khan
     49! @author Sabine Banzhaf
     50! @author Emmanuele Russo
    4551!
    4652!------------------------------------------------------------------------------!
    … …  
    8894    LOGICAL       :: do_emis                               = .FALSE.               !< Flag for turning on chemistry emissions
    8995    LOGICAL       :: cs_pr_namelist_found                  = .FALSE.               !< Namelist parameter: Names of t
     96    LOGICAL       :: do_depo                               = .FALSE.               !< namelist parameter for activation of deposition calculation
    9097
    9198
    … …  
    161168                                                                                   ! matching between the model and the input files
    162169
     170   
     171    !-- Selected atomic/molecular weights:
     172   
     173    REAL, PARAMETER        ::  xm_H     =    1.00790e-3           ! kg/mol
     174    REAL, PARAMETER        ::  xm_N     =   14.00670e-3           ! kg/mol
     175    REAL, PARAMETER        ::  xm_C     =   12.01115e-3           ! kg/mol
     176    REAL, PARAMETER        ::  xm_S     =   32.06400e-3           ! kg/mol
     177    REAL, PARAMETER        ::  xm_O     =   15.99940e-3           ! kg/mol
     178    REAL, PARAMETER        ::  xm_F     =   18.99840e-3           ! kg/mol
     179    REAL, PARAMETER        ::  xm_Na    =   22.98977e-3           ! kg/mol
     180    REAL, PARAMETER        ::  xm_Cl    =   35.45300e-3           ! kg/mol
     181    REAL, PARAMETER        ::  xm_Rn222 =  222.00000e-3           ! kg/mol
     182    REAL, PARAMETER        ::  xm_Pb210 =  210.00000e-3           ! kg/mol
     183    REAL, PARAMETER        ::  xm_Ca    =   40.07800e-3           ! kg/mol
     184    REAL, PARAMETER        ::  xm_K     =   39.09800e-3           ! kg/mol
     185    REAL, PARAMETER        ::  xm_Mg    =   24.30500e-3           ! kg/mol
     186    REAL, PARAMETER        ::  xm_Pb    =  207.20000e-3           ! kg/mol
     187    REAL, PARAMETER        ::  xm_Cd    =  112.41000e-3           ! kg/mol
     188   
     189    REAL, PARAMETER        ::  xm_h2o   = xm_H * 2 + xm_O         ! kg/mol
     190    REAL, PARAMETER        ::  xm_o3    = xm_O * 3                ! kg/mol
     191    REAL, PARAMETER        ::  xm_N2O5  = xm_N * 2 + xm_O * 5     ! kg/mol
     192    REAL, PARAMETER        ::  xm_HNO3  = xm_H + xm_N + xm_O * 3  ! kg/mol
     193    REAL, PARAMETER        ::  xm_NH4   = xm_N + xm_H * 4         ! kg/mol
     194    REAL, PARAMETER        ::  xm_SO4   = xm_S + xm_O * 4         ! kg/mol
     195    REAL, PARAMETER        ::  xm_NO3   = xm_N + xm_O * 3         ! kg/mol
     196    REAL, PARAMETER        ::  xm_CO2   = xm_C + xm_O * 2         ! kg/mol
     197   
     198    ! mass of air
     199    REAL, PARAMETER        ::  xm_air   =  28.964e-3              ! kg/mol
     200       
     201    ! dummy weight, used for complex molecules:
     202    REAL, PARAMETER        ::  xm_dummy =  1000.0e-3              ! kg/mol
     203
     204   
    163205    SAVE
    164206 END MODULE chem_modules

  • palm/trunk/SOURCE/chemistry_model_mod.f90

    r3449 r3458  
    2727! -----------------
    2828! $Id$
     29! from chemistry branch r3443, banzhafs, basit:
     30! replace surf_lsm_h%qv1(m) by q(k,j,i) for mixing ratio in chem_depo
     31! bug fix in chem_depo: allow different surface fractions for one
     32! surface element and set lai to zero for non vegetated surfaces
     33! bug fixed in chem_data_output_2d
     34! bug fix in chem_depo subroutine
     35! added code on deposition of gases and particles
     36! removed cs_profile_name from chem_parin
     37! bug fixed in output profiles and code cleaned
     38!
     39! 3449 2018-10-29 19:36:56Z suehring
    2940! additional output - merged from branch resler
    3041!
    … …  
    151162! @author Klaus Ketelsen
    152163! @author Basit Khan
     164! @author Sabine Banzhaf
    153165!
    154166!
    … …  
    197209                                 initializing_actions, message_string,                             &
    198210                                 omega, tsc, intermediate_timestep_count,                          &
    199                                  intermediate_timestep_count_max,                                  &
     211                                 intermediate_timestep_count_max, max_pr_user,                     &
    200212                                 timestep_scheme, use_prescribed_profile_data, ws_scheme_sca 
    201213   USE arrays_3d,          ONLY: exner, hyp, pt, q, ql, rdf_sc, tend, zu
    … …  
    271283                                         ! 'fastj'  (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED
    272284
     285   
     286    !-----------------------------------------------------------------------
     287    ! Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
     288    !
     289    INTEGER(iwp), PARAMETER :: nlu_dep = 15                   ! Number of DEPAC landuse classes (lu's)
     290    INTEGER(iwp), PARAMETER :: ncmp = 10                      ! Number of DEPAC gas components
     291    !------------------------------------------------------------------------
     292    ! DEPAC landuse classes as defined in LOTOS-EUROS model v2.1                             
     293    !
     294    INTEGER(iwp)  ::  ilu_grass              = 1       
     295    INTEGER(iwp)  ::  ilu_arable             = 2       
     296    INTEGER(iwp)  ::  ilu_permanent_crops    = 3         
     297    INTEGER(iwp)  ::  ilu_coniferous_forest  = 4         
     298    INTEGER(iwp)  ::  ilu_deciduous_forest   = 5         
     299    INTEGER(iwp)  ::  ilu_water_sea          = 6       
     300    INTEGER(iwp)  ::  ilu_urban              = 7       
     301    INTEGER(iwp)  ::  ilu_other              = 8 
     302    INTEGER(iwp)  ::  ilu_desert             = 9 
     303    INTEGER(iwp)  ::  ilu_ice                = 10
     304    INTEGER(iwp)  ::  ilu_savanna            = 11
     305    INTEGER(iwp)  ::  ilu_tropical_forest    = 12
     306    INTEGER(iwp)  ::  ilu_water_inland       = 13
     307    INTEGER(iwp)  ::  ilu_mediterrean_scrub  = 14
     308    INTEGER(iwp)  ::  ilu_semi_natural_veg   = 15
     309    !
     310    !--------------------------------------------------------------------------
     311    ! NH3/SO2 ratio regimes:
     312    INTEGER, PARAMETER  ::  iratns_low  = 1       ! low ratio NH3/SO2
     313    INTEGER, PARAMETER  ::  iratns_high = 2       ! high ratio NH3/SO2
     314    INTEGER, PARAMETER  ::  iratns_very_low = 3   ! very low ratio NH3/SO2
     315    ! Default:
     316    INTEGER, PARAMETER  ::  iratns_default = iratns_low
     317    !
     318    ! Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) naar W m-2
     319    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: alpha   =(/0.009,0.009, 0.009,0.006,0.006, -999., -999.,0.009,-999.,-999.,0.009,0.006,-999.,0.009,0.008/)*4.57
     320    !
     321    ! Set temperatures per land use for F_temp
     322    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Tmin =   (/12.0, 12.0,  12.0,  0.0,  0.0, -999., -999., 12.0, -999., -999., 12.0,  0.0, -999., 12.0,  8.0/)
     323    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Topt =   (/26.0, 26.0,  26.0, 18.0, 20.0, -999., -999., 26.0, -999., -999., 26.0, 20.0, -999., 26.0, 24.0 /)
     324    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Tmax =   (/40.0, 40.0,  40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999., 40.0, 35.0, -999., 40.0, 39.0 /)
     325    !
     326    ! Set F_min:
     327    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,0.1,-999.,0.01, 0.04/)
     328   
     329    ! Set maximal conductance (m/s)
     330    ! (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from  mmol O3/m2/s to m/s
     331    ! Could be refined to a function of T and P. in Jones
     332    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: g_max   =(/270., 300.,  300., 140., 150.,  -999., -999.,270., -999.,-999.,270., 150.,-999.,300., 422./)/41000
     333    !
     334    ! Set max, min for vapour pressure deficit vpd;
     335    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,1.0,  -999.,0.9, 2.8/)
     336    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: vpd_min =(/3.0,  2.8,   2.8,  3.0,  3.25,  -999., -999.,3.0,  -999.,-999.,3.0,3.25, -999.,2.8, 4.5/)
     337    !
     338    ! 
     339    !------------------------------------------------------------------------
     340
     341   
    273342    PUBLIC nest_chemistry
    274343    PUBLIC nspec
    … …  
    357426    END INTERFACE chem_wrd_local
    358427
     428    INTERFACE chem_depo
     429       MODULE PROCEDURE chem_depo
     430    END INTERFACE chem_depo
     431
     432    INTERFACE drydepos_gas_depac
     433       MODULE PROCEDURE drydepos_gas_depac
     434    END INTERFACE drydepos_gas_depac
     435
     436    INTERFACE rc_special
     437       MODULE PROCEDURE rc_special
     438    END INTERFACE rc_special
     439
     440    INTERFACE  rc_gw
     441       MODULE PROCEDURE rc_gw
     442    END INTERFACE rc_gw
     443
     444    INTERFACE rw_so2
     445       MODULE PROCEDURE rw_so2 
     446    END INTERFACE rw_so2
     447
     448    INTERFACE rw_nh3_sutton
     449       MODULE PROCEDURE rw_nh3_sutton
     450    END INTERFACE rw_nh3_sutton
     451
     452    INTERFACE rw_constant
     453       MODULE PROCEDURE rw_constant
     454    END INTERFACE rw_constant
     455
     456    INTERFACE rc_gstom
     457       MODULE PROCEDURE rc_gstom
     458    END INTERFACE rc_gstom
     459
     460    INTERFACE rc_gstom_emb
     461       MODULE PROCEDURE rc_gstom_emb
     462    END INTERFACE rc_gstom_emb
     463
     464    INTERFACE par_dir_diff
     465       MODULE PROCEDURE par_dir_diff
     466    END INTERFACE par_dir_diff
     467
     468    INTERFACE rc_get_vpd
     469       MODULE PROCEDURE rc_get_vpd
     470    END INTERFACE rc_get_vpd
     471
     472    INTERFACE rc_gsoil_eff
     473       MODULE PROCEDURE rc_gsoil_eff
     474    END INTERFACE rc_gsoil_eff
     475
     476    INTERFACE rc_rinc
     477       MODULE PROCEDURE rc_rinc
     478    END INTERFACE rc_rinc
     479
     480    INTERFACE rc_rctot
     481       MODULE PROCEDURE rc_rctot 
     482    END INTERFACE rc_rctot
     483
     484    INTERFACE rc_comp_point_rc_eff
     485       MODULE PROCEDURE rc_comp_point_rc_eff
     486    END INTERFACE rc_comp_point_rc_eff
     487
     488    INTERFACE drydepo_aero_zhang_vd
     489       MODULE PROCEDURE drydepo_aero_zhang_vd
     490    END INTERFACE drydepo_aero_zhang_vd
     491
     492    INTERFACE get_rb_cell
     493       MODULE PROCEDURE  get_rb_cell
     494    END INTERFACE get_rb_cell
     495
     496   
    359497
    360498    PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, &
    … …  
    364502           chem_init_profiles, chem_integrate, chem_parin,                      &
    365503           chem_prognostic_equations, chem_rrd_local,                           &
    366            chem_statistics, chem_swap_timelevel, chem_wrd_local
     504           chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo
    367505
    368506
    … …  
    838976!            It is possible  to plant PM10 and PM25 into the gasphase chemistry code
    839977!            as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
    840 !            set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
     978!            set unit to kilograms per m**3 for PM10 and PM25 (PM2.5)
    841979             IF (spec_name(1:2) == 'PM')   THEN 
    842980               unit = 'kg m-3'
    … …  
    852990
    853991
    854       RETURN
     992       RETURN
    855993    END SUBROUTINE chem_check_data_output
    856994!
    … …  
    10721210                         DO  k = nzb_do, nzt_do
    10731211                              local_pf(i,j,k) = MERGE(                                                &
    1074                                                  chem_species(lsp)%conc(k,j,i),                       &
     1212                                                 chem_species(lsp)%conc_av(k,j,i),                    &
    10751213                                                 REAL( fill_value, KIND = wp ),                       &
    10761214                                                 BTEST( wall_flags_0(k,j,i), 0 ) )
    … …  
    11041242       USE surface_mod
    11051243
    1106        !USE chem_modules, ONLY: nvar
    11071244
    11081245       IMPLICIT NONE
    … …  
    12391376
    12401377       spec_name = TRIM( variable(4:) )
    1241        !av == 0
    12421378
    12431379       DO lsp=1,nspec
    … …  
    18702006                                        cs_name,                          &
    18712007                                        cs_profile,                       &
    1872                                         cs_profile_name,                  &
    18732008                                        cs_surface,                       &
    18742009                                        decycle_chem_lr,                  &
    18752010                                        decycle_chem_ns,                  &           
    1876                                         decycle_method,                   &
     2011                                        decycle_method,                   &
     2012                                        do_depo,                          &
    18772013                                        emiss_factor_main,                &
    18782014                                        emiss_factor_side,                &                     
    … …  
    23652501                DO lpr = 1, cs_pr_count
    23662502
    2367                    sums_l(k,pr_palm+lpr,tn) = sums_l(k,pr_palm+lpr,tn) +    &
     2503                   sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +    &
    23682504                                                 chem_species(cs_pr_index(lpr))%conc(k,j,i) *       &
    23692505                                                 rmask(j,i,sr)  *                                   &
    … …  
    24322568
    24332569    END SUBROUTINE chem_wrd_local
     2570
     2571
     2572
     2573!-------------------------------------------------------------------------------!
     2574!
     2575! Description:
     2576! ------------
     2577!> Subroutine to calculate the deposition of gases and PMs. For now deposition
     2578!> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
     2579!> considered. The deposition of particles is derived following Zhang et al.,
     2580!> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
     2581!>     
     2582!> @TODO: Consider deposition on vertical surfaces   
     2583!> @TODO: Consider overlaying horizontal surfaces
     2584!> @TODO: Check error messages
     2585!-------------------------------------------------------------------------------!
     2586     
     2587    SUBROUTINE chem_depo (i,j)
     2588
     2589
     2590      USE control_parameters,                                                 &   
     2591           ONLY:  dt_3d, intermediate_timestep_count, latitude
     2592
     2593      USE arrays_3d,                                                          &
     2594           ONLY:  dzw, rho_air_zw
     2595
     2596      USE date_and_time_mod,                                                  &
     2597           ONLY: day_of_year
     2598
     2599      USE surface_mod,                                                          &
     2600           ONLY: surf_lsm_h, surf_usm_h, surf_type, ind_veg_wall,               &
     2601           ind_pav_green, ind_wat_win
     2602
     2603      USE radiation_model_mod,                                                   &
     2604           ONLY: zenith
     2605
     2606
     2607
     2608      IMPLICIT NONE
     2609
     2610      INTEGER,INTENT(IN)       :: i,j
     2611
     2612      INTEGER(iwp) ::  k                                        ! matching k to surface m at i,j
     2613      INTEGER(iwp) ::  lsp                                      ! running index for chem spcs.
     2614      INTEGER(iwp) ::  lu                                       ! running index for landuse classes
     2615      INTEGER(iwp) ::  lu_palm                                  ! index of PALM LSM vegetation_type at current surface element
     2616      INTEGER(iwp) ::  lup_palm                                 ! index of PALM LSM pavement_type at current surface element
     2617      INTEGER(iwp) ::  luw_palm                                 ! index of PALM LSM water_type at current surface element
     2618      INTEGER(iwp) ::  luu_palm                                 ! index of PALM USM walls/roofs at current surface element
     2619      INTEGER(iwp) ::  lug_palm                                 ! index of PALM USM green walls/roofs at current surface element
     2620      INTEGER(iwp) ::  lud_palm                                 ! index of PALM USM windows at current surface element
     2621      INTEGER(iwp) ::  lu_dep                                   ! matching DEPAC LU to lu_palm
     2622      INTEGER(iwp) ::  lup_dep                                  ! matching DEPAC LU to lup_palm
     2623      INTEGER(iwp) ::  luw_dep                                  ! matching DEPAC LU to luw_palm
     2624      INTEGER(iwp) ::  luu_dep                                  ! matching DEPAC LU to luu_palm
     2625      INTEGER(iwp) ::  lug_dep                                  ! matching DEPAC LU to lug_palm
     2626      INTEGER(iwp) ::  lud_dep                                  ! matching DEPAC LU to lud_palm
     2627      INTEGER(iwp) ::  m                                        ! index for horizontal surfaces
     2628
     2629      INTEGER(iwp) ::  pspec                                    ! running index
     2630      INTEGER(iwp) ::  i_pspec                                  ! index for matching depac gas component
     2631
     2632
     2633      ! Vegetation                                              !Match to DEPAC
     2634      INTEGER(iwp)  ::  ind_lu_user = 0                         ! --> ERROR 
     2635      INTEGER(iwp)  ::  ind_lu_b_soil = 1                       ! --> ilu_desert
     2636      INTEGER(iwp)  ::  ind_lu_mixed_crops = 2                  ! --> ilu_arable
     2637      INTEGER(iwp)  ::  ind_lu_s_grass = 3                      ! --> ilu_grass
     2638      INTEGER(iwp)  ::  ind_lu_ev_needle_trees = 4              ! --> ilu_coniferous_forest
     2639      INTEGER(iwp)  ::  ind_lu_de_needle_trees = 5              ! --> ilu_coniferous_forest
     2640      INTEGER(iwp)  ::  ind_lu_ev_broad_trees = 6               ! --> ilu_tropical_forest
     2641      INTEGER(iwp)  ::  ind_lu_de_broad_trees = 7               ! --> ilu_deciduous_forest
     2642      INTEGER(iwp)  ::  ind_lu_t_grass = 8                      ! --> ilu_grass
     2643      INTEGER(iwp)  ::  ind_lu_desert = 9                       ! --> ilu_desert
     2644      INTEGER(iwp)  ::  ind_lu_tundra = 10                      ! --> ilu_other
     2645      INTEGER(iwp)  ::  ind_lu_irr_crops = 11                   ! --> ilu_arable
     2646      INTEGER(iwp)  ::  ind_lu_semidesert = 12                  ! --> ilu_other
     2647      INTEGER(iwp)  ::  ind_lu_ice = 13                         ! --> ilu_ice
     2648      INTEGER(iwp)  ::  ind_lu_marsh = 14                       ! --> ilu_other
     2649      INTEGER(iwp)  ::  ind_lu_ev_shrubs = 15                   ! --> ilu_mediterrean_scrub
     2650      INTEGER(iwp)  ::  ind_lu_de_shrubs = 16                   ! --> ilu_mediterrean_scrub
     2651      INTEGER(iwp)  ::  ind_lu_mixed_forest = 17                ! --> ilu_coniferous_forest (ave(decid+conif))
     2652      INTEGER(iwp)  ::  ind_lu_intrup_forest = 18               ! --> ilu_other (ave(other+decid))
     2653
     2654      ! Water
     2655      INTEGER(iwp)  ::  ind_luw_user = 0                        ! --> ERROR
     2656      INTEGER(iwp)  ::  ind_luw_lake = 1                        ! --> ilu_water_inland
     2657      INTEGER(iwp)  ::  ind_luw_river = 2                       ! --> ilu_water_inland
     2658      INTEGER(iwp)  ::  ind_luw_ocean = 3                       ! --> ilu_water_sea
     2659      INTEGER(iwp)  ::  ind_luw_pond = 4                        ! --> ilu_water_inland
     2660      INTEGER(iwp)  ::  ind_luw_fountain = 5                    ! --> ilu_water_inland
     2661
     2662      ! Pavement
     2663      INTEGER(iwp)  ::  ind_lup_user = 0                        ! --> ERROR
     2664      INTEGER(iwp)  ::  ind_lup_asph_conc = 1                   ! --> ilu_desert
     2665      INTEGER(iwp)  ::  ind_lup_asph = 2                        ! --> ilu_desert
     2666      INTEGER(iwp)  ::  ind_lup_conc = 3                        ! --> ilu_desert
     2667      INTEGER(iwp)  ::  ind_lup_sett = 4                        ! --> ilu_desert
     2668      INTEGER(iwp)  ::  ind_lup_pav_stones = 5                  ! --> ilu_desert
     2669      INTEGER(iwp)  ::  ind_lup_cobblest = 6                    ! --> ilu_desert
     2670      INTEGER(iwp)  ::  ind_lup_metal = 7                       ! --> ilu_desert
     2671      INTEGER(iwp)  ::  ind_lup_wood = 8                        ! --> ilu_desert
     2672      INTEGER(iwp)  ::  ind_lup_gravel = 9                      ! --> ilu_desert
     2673      INTEGER(iwp)  ::  ind_lup_f_gravel = 10                   ! --> ilu_desert
     2674      INTEGER(iwp)  ::  ind_lup_pebblest = 11                   ! --> ilu_desert
     2675      INTEGER(iwp)  ::  ind_lup_woodchips = 12                  ! --> ilu_desert
     2676      INTEGER(iwp)  ::  ind_lup_tartan = 13                     ! --> ilu_desert
     2677      INTEGER(iwp)  ::  ind_lup_art_turf = 14                   ! --> ilu_desert
     2678      INTEGER(iwp)  ::  ind_lup_clay = 15                       ! --> ilu_desert
     2679
     2680
     2681
     2682      !-- Particle parameters according to the respective aerosol classes (PM25, PM10)
     2683      INTEGER(iwp) ::  ind_p_size = 1     !< index for partsize in particle_pars
     2684      INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
     2685      INTEGER(iwp) ::  ind_p_slip = 3     !< index for slipcor in particle_pars
     2686
     2687      INTEGER(iwp) ::  part_type
     2688
     2689      INTEGER(iwp) :: nwet                ! wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     2690
     2691      REAL(kind=wp)  :: dt_chem
     2692      REAL(kind=wp)  :: dh               
     2693      REAL(kind=wp)  :: inv_dh
     2694      REAL(kind=wp)  :: dt_dh
     2695
     2696      REAL(kind=wp)  :: dens              ! Density at lowest layer at i,j 
     2697      REAL(kind=wp)  :: r_aero_lu         ! aerodynamic resistance (s/m) at current surface element
     2698      REAL(kind=wp)  :: ustar_lu          ! ustar at current surface element
     2699      REAL(kind=wp)  :: z0h_lu            ! roughness length for heat at current surface element
     2700      REAL(kind=wp)  :: glrad             ! rad_sw_in at current surface element
     2701      REAL(kind=wp)  :: ppm_to_ugm3       ! conversion factor
     2702      REAL(kind=wp)  :: relh              ! relative humidity at current surface element
     2703      REAL(kind=wp)  :: lai               ! leaf area index at current surface element
     2704      REAL(kind=wp)  :: sai               ! surface area index at current surface element assumed to be lai + 1
     2705
     2706      REAL(kind=wp)  :: slinnfac       
     2707      REAL(kind=wp)  :: visc              ! Viscosity
     2708      REAL(kind=wp)  :: vs                ! Sedimentation velocity
     2709      REAL(kind=wp)  :: vd_lu             ! deposition velocity (m/s)
     2710      REAL(kind=wp)  :: Rs                ! Sedimentaion resistance (s/m)
     2711      REAL(kind=wp)  :: Rb                ! quasi-laminar boundary layer resistance (s/m)
     2712      REAL(kind=wp)  :: Rc_tot            ! total canopy resistance Rc (s/m)
     2713
     2714      REAL(kind=wp)  :: catm              ! gasconc. at i,j,k=1 in ug/m3
     2715      REAL(kind=wp)  :: diffc             ! Diffusivity
     2716
     2717
     2718      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_lu       ! budget for LSM vegetation type at current surface element
     2719      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_lup      ! budget for LSM pavement type at current surface element
     2720      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_luw      ! budget for LSM water type at current surface element
     2721      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_luu       ! budget for USM walls/roofs at current surface element
     2722      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_lug      ! budget for USM green surfaces at current surface element
     2723      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_lud      ! budget for USM windows at current surface element
     2724      REAL(kind=wp),DIMENSION(nspec)            :: bud_now          ! overall budget at current surface element
     2725      REAL(kind=wp),DIMENSION(nspec)            :: cc               ! concentration i,j,k
     2726      REAL(kind=wp),DIMENSION(nspec)            :: ccomp_tot        ! total compensation point (ug/m3) (= 0 for species that don't have a compensation)
     2727      ! For now kept to zero for all species!
     2728
     2729      REAL(kind=wp)                             :: ttemp          ! temperatur at i,j,k
     2730      REAL(kind=wp)                             :: ts             ! surface temperatur in degrees celsius
     2731      REAL(kind=wp)                             :: qv_tmp         ! surface mixing ratio at current surface element
     2732
     2733
     2734      !-- Particle parameters (PM10 (1), PM25 (2))
     2735      !-- partsize (diameter in m) ,    rhopart (density in kg/m3),     slipcor (slip correction factor DIMENSIONless, Seinfeld and Pandis 2006, Table 9.3)
     2736      REAL(wp), DIMENSION(1:3,1:2), PARAMETER   :: particle_pars = RESHAPE( (/ &
     2737           8.0e-6_wp, 1.14e3_wp, 1.016_wp, & !  1
     2738           0.7e-6_wp, 1.14e3_wp, 1.082_wp &  !  2
     2739           /), (/ 3, 2 /) )
     2740
     2741
     2742      LOGICAL                                   ::  match_lsm     ! flag indicating natural-type surface
     2743      LOGICAL                                   ::  match_usm     ! flag indicating urban-type surface
     2744
     2745
     2746
     2747      !-- List of names of possible tracers
     2748      CHARACTER(len=*), PARAMETER  ::  pspecnames(69) = (/ &
     2749           'NO2           ', &    ! NO2
     2750           'NO            ', &    ! NO
     2751           'O3            ', &    ! O3
     2752           'CO            ', &    ! CO
     2753           'form          ', &    ! FORM
     2754           'ald           ', &    ! ALD
     2755           'pan           ', &    ! PAN
     2756           'mgly          ', &    ! MGLY
     2757           'par           ', &    ! PAR
     2758           'ole           ', &    ! OLE
     2759           'eth           ', &    ! ETH
     2760           'tol           ', &    ! TOL
     2761           'cres          ', &    ! CRES
     2762           'xyl           ', &    ! XYL
     2763           'SO4a_f        ', &    ! SO4a_f
     2764           'SO2           ', &    ! SO2
     2765           'HNO2          ', &    ! HNO2
     2766           'CH4           ', &    ! CH4
     2767           'NH3           ', &    ! NH3
     2768           'NO3           ', &    ! NO3
     2769           'OH            ', &    ! OH
     2770           'HO2           ', &    ! HO2
     2771           'N2O5          ', &    ! N2O5
     2772           'SO4a_c        ', &    ! SO4a_c
     2773           'NH4a_f        ', &    ! NH4a_f
     2774           'NO3a_f        ', &    ! NO3a_f
     2775           'NO3a_c        ', &    ! NO3a_c
     2776           'C2O3          ', &    ! C2O3
     2777           'XO2           ', &    ! XO2
     2778           'XO2N          ', &    ! XO2N
     2779           'cro           ', &    ! CRO
     2780           'HNO3          ', &    ! HNO3
     2781           'H2O2          ', &    ! H2O2
     2782           'iso           ', &    ! ISO
     2783           'ispd          ', &    ! ISPD
     2784           'to2           ', &    ! TO2
     2785           'open          ', &    ! OPEN
     2786           'terp          ', &    ! TERP
     2787           'ec_f          ', &    ! EC_f
     2788           'ec_c          ', &    ! EC_c
     2789           'pom_f         ', &    ! POM_f
     2790           'pom_c         ', &    ! POM_c
     2791           'ppm_f         ', &    ! PPM_f
     2792           'ppm_c         ', &    ! PPM_c
     2793           'na_ff         ', &    ! Na_ff
     2794           'na_f          ', &    ! Na_f
     2795           'na_c          ', &    ! Na_c
     2796           'na_cc         ', &    ! Na_cc
     2797           'na_ccc        ', &    ! Na_ccc
     2798           'dust_ff       ', &    ! dust_ff
     2799           'dust_f        ', &    ! dust_f
     2800           'dust_c        ', &    ! dust_c
     2801           'dust_cc       ', &    ! dust_cc
     2802           'dust_ccc      ', &    ! dust_ccc
     2803           'tpm10         ', &    ! tpm10
     2804           'tpm25         ', &    ! tpm25
     2805           'tss           ', &    ! tss
     2806           'tdust         ', &    ! tdust
     2807           'tc            ', &    ! tc
     2808           'tcg           ', &    ! tcg
     2809           'tsoa          ', &    ! tsoa
     2810           'tnmvoc        ', &    ! tnmvoc
     2811           'SOxa          ', &    ! SOxa
     2812           'NOya          ', &    ! NOya
     2813           'NHxa          ', &    ! NHxa
     2814           'NO2_obs       ', &    ! NO2_obs
     2815           'tpm10_biascorr', &    ! tpm10_biascorr
     2816           'tpm25_biascorr', &    ! tpm25_biascorr
     2817           'O3_biascorr   ' /)    ! o3_biascorr
     2818
     2819
     2820
     2821      ! tracer mole mass:
     2822      REAL, PARAMETER  ::  specmolm(69) = (/ &
     2823           xm_O * 2 + xm_N, &                         ! NO2
     2824           xm_O + xm_N, &                             ! NO
     2825           xm_O * 3, &                                ! O3
     2826           xm_C + xm_O, &                             ! CO
     2827           xm_H * 2 + xm_C + xm_O, &                  ! FORM
     2828           xm_H * 3 + xm_C * 2 + xm_O, &              ! ALD
     2829           xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, &   ! PAN
     2830           xm_H * 4 + xm_C * 3 + xm_O * 2, &          ! MGLY
     2831           xm_H * 3 + xm_C, &                         ! PAR
     2832           xm_H * 3 + xm_C * 2, &                     ! OLE
     2833           xm_H * 4 + xm_C * 2, &                     ! ETH
     2834           xm_H * 8 + xm_C * 7, &                     ! TOL
     2835           xm_H * 8 + xm_C * 7 + xm_O, &              ! CRES
     2836           xm_H * 10 + xm_C * 8, &                    ! XYL
     2837           xm_S + xm_O * 4, &                         ! SO4a_f
     2838           xm_S + xm_O * 2, &                         ! SO2
     2839           xm_H + xm_O * 2 + xm_N, &                  ! HNO2
     2840           xm_H * 4 + xm_C, &                         ! CH4
     2841           xm_H * 3 + xm_N, &                         ! NH3
     2842           xm_O * 3 + xm_N, &                         ! NO3
     2843           xm_H + xm_O, &                             ! OH
     2844           xm_H + xm_O * 2, &                         ! HO2
     2845           xm_O * 5 + xm_N * 2, &                     ! N2O5
     2846           xm_S + xm_O * 4, &                         ! SO4a_c
     2847           xm_H * 4 + xm_N, &                         ! NH4a_f
     2848           xm_O * 3 + xm_N, &                         ! NO3a_f
     2849           xm_O * 3 + xm_N, &                         ! NO3a_c
     2850           xm_C * 2 + xm_O * 3, &                     ! C2O3
     2851           xm_dummy, &                                ! XO2
     2852           xm_dummy, &                                ! XO2N
     2853           xm_dummy, &                                ! CRO
     2854           xm_H + xm_O * 3 + xm_N, &                  ! HNO3
     2855           xm_H * 2 + xm_O * 2, &                     ! H2O2
     2856           xm_H * 8 + xm_C * 5, &                     ! ISO
     2857           xm_dummy, &                                ! ISPD
     2858           xm_dummy, &                                ! TO2
     2859           xm_dummy, &                                ! OPEN
     2860           xm_H * 16 + xm_C * 10, &                   ! TERP
     2861           xm_dummy, &                                ! EC_f
     2862           xm_dummy, &                                ! EC_c
     2863           xm_dummy, &                                ! POM_f
     2864           xm_dummy, &                                ! POM_c
     2865           xm_dummy, &                                ! PPM_f
     2866           xm_dummy, &                                ! PPM_c
     2867           xm_Na, &                                   ! Na_ff
     2868           xm_Na, &                                   ! Na_f
     2869           xm_Na, &                                   ! Na_c
     2870           xm_Na, &                                   ! Na_cc
     2871           xm_Na, &                                   ! Na_ccc
     2872           xm_dummy, &                                ! dust_ff
     2873           xm_dummy, &                                ! dust_f
     2874           xm_dummy, &                                ! dust_c
     2875           xm_dummy, &                                ! dust_cc
     2876           xm_dummy, &                                ! dust_ccc
     2877           xm_dummy, &                                ! tpm10
     2878           xm_dummy, &                                ! tpm25
     2879           xm_dummy, &                                ! tss
     2880           xm_dummy, &                                ! tdust
     2881           xm_dummy, &                                ! tc
     2882           xm_dummy, &                                ! tcg
     2883           xm_dummy, &                                ! tsoa
     2884           xm_dummy, &                                ! tnmvoc
     2885           xm_dummy, &                                ! SOxa
     2886           xm_dummy, &                                ! NOya
     2887           xm_dummy, &                                ! NHxa
     2888           xm_O * 2 + xm_N, &                         ! NO2_obs
     2889           xm_dummy, &                                ! tpm10_biascorr
     2890           xm_dummy, &                                ! tpm25_biascorr
     2891           xm_O * 3 /)                                ! o3_biascorr
     2892
     2893
     2894
     2895      ! Initialize surface element m
     2896      m = 0
     2897      k = 0
     2898      ! LSM or USM surface present at i,j:
     2899      ! Default surfaces are NOT considered for deposition
     2900      match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
     2901      match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
     2902
     2903
     2904      !!!!!!!
     2905      ! LSM !
     2906      !!!!!!!
     2907
     2908      IF ( match_lsm )  THEN
     2909
     2910         ! Get surface element information at i,j:
     2911         m = surf_lsm_h%start_index(j,i)
     2912
     2913         ! Get corresponding k
     2914         k = surf_lsm_h%k(m)
     2915
     2916         ! Get needed variables for surface element m
     2917         ustar_lu  = surf_lsm_h%us(m)
     2918         z0h_lu    = surf_lsm_h%z0h(m)
     2919         r_aero_lu = surf_lsm_h%r_a(m)
     2920         glrad     = surf_lsm_h%rad_sw_in(m)
     2921         lai = surf_lsm_h%lai(m)
     2922         sai = lai + 1
     2923
     2924         ! For small grid spacing neglect R_a
     2925         IF (dzw(k) <= 1.0) THEN
     2926            r_aero_lu = 0.0_wp
     2927         ENDIF
     2928
     2929         !Initialize lu's
     2930         lu_palm = 0
     2931         lu_dep = 0
     2932         lup_palm = 0
     2933         lup_dep = 0
     2934         luw_palm = 0
     2935         luw_dep = 0
     2936
     2937         !Initialize budgets
     2938         bud_now_lu  = 0.0_wp
     2939         bud_now_lup = 0.0_wp
     2940         bud_now_luw = 0.0_wp
     2941
     2942
     2943         ! Get land use for i,j and assign to DEPAC lu
     2944         IF (surf_lsm_h%frac(ind_veg_wall,m) > 0) THEN
     2945            lu_palm = surf_lsm_h%vegetation_type(m)
     2946            IF (lu_palm == ind_lu_user) THEN
     2947               message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
     2948               CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 )
     2949            ELSEIF (lu_palm == ind_lu_b_soil) THEN
     2950               lu_dep = 9
     2951            ELSEIF (lu_palm == ind_lu_mixed_crops) THEN
     2952               lu_dep = 2
     2953            ELSEIF (lu_palm == ind_lu_s_grass) THEN
     2954               lu_dep = 1
     2955            ELSEIF (lu_palm == ind_lu_ev_needle_trees) THEN
     2956               lu_dep = 4
     2957            ELSEIF (lu_palm == ind_lu_de_needle_trees) THEN
     2958               lu_dep = 4
     2959            ELSEIF (lu_palm == ind_lu_ev_broad_trees) THEN
     2960               lu_dep = 12
     2961            ELSEIF (lu_palm == ind_lu_de_broad_trees) THEN
     2962               lu_dep = 5
     2963            ELSEIF (lu_palm == ind_lu_t_grass) THEN
     2964               lu_dep = 1
     2965            ELSEIF (lu_palm == ind_lu_desert) THEN
     2966               lu_dep = 9
     2967            ELSEIF (lu_palm == ind_lu_tundra ) THEN
     2968               lu_dep = 8
     2969            ELSEIF (lu_palm == ind_lu_irr_crops) THEN
     2970               lu_dep = 2
     2971            ELSEIF (lu_palm == ind_lu_semidesert) THEN
     2972               lu_dep = 8
     2973            ELSEIF (lu_palm == ind_lu_ice) THEN
     2974               lu_dep = 10
     2975            ELSEIF (lu_palm == ind_lu_marsh) THEN
     2976               lu_dep = 8
     2977            ELSEIF (lu_palm == ind_lu_ev_shrubs) THEN
     2978               lu_dep = 14
     2979            ELSEIF (lu_palm == ind_lu_de_shrubs ) THEN
     2980               lu_dep = 14
     2981            ELSEIF (lu_palm == ind_lu_mixed_forest) THEN
     2982               lu_dep = 4
     2983            ELSEIF (lu_palm == ind_lu_intrup_forest) THEN
     2984               lu_dep = 8     
     2985            ENDIF
     2986         ENDIF
     2987
     2988         IF (surf_lsm_h%frac(ind_pav_green,m) > 0) THEN
     2989            lup_palm = surf_lsm_h%pavement_type(m)
     2990            IF (lup_palm == ind_lup_user) THEN
     2991               message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
     2992               CALL message( 'chem_depo', 'CM0452', 1, 2, 0, 6, 0 )
     2993            ELSEIF (lup_palm == ind_lup_asph_conc) THEN
     2994               lup_dep = 9
     2995            ELSEIF (lup_palm == ind_lup_asph) THEN
     2996               lup_dep = 9
     2997            ELSEIF (lup_palm ==  ind_lup_conc) THEN
     2998               lup_dep = 9
     2999            ELSEIF (lup_palm ==  ind_lup_sett) THEN
     3000               lup_dep = 9
     3001            ELSEIF (lup_palm == ind_lup_pav_stones) THEN
     3002               lup_dep = 9
     3003            ELSEIF (lup_palm == ind_lup_cobblest) THEN
     3004               lup_dep = 9       
     3005            ELSEIF (lup_palm == ind_lup_metal) THEN
     3006               lup_dep = 9
     3007            ELSEIF (lup_palm == ind_lup_wood) THEN
     3008               lup_dep = 9   
     3009            ELSEIF (lup_palm == ind_lup_gravel) THEN
     3010               lup_dep = 9
     3011            ELSEIF (lup_palm == ind_lup_f_gravel) THEN
     3012               lup_dep = 9
     3013            ELSEIF (lup_palm == ind_lup_pebblest) THEN
     3014               lup_dep = 9
     3015            ELSEIF (lup_palm == ind_lup_woodchips) THEN
     3016               lup_dep = 9
     3017            ELSEIF (lup_palm == ind_lup_tartan) THEN
     3018               lup_dep = 9
     3019            ELSEIF (lup_palm == ind_lup_art_turf) THEN
     3020               lup_dep = 9
     3021            ELSEIF (lup_palm == ind_lup_clay) THEN
     3022               lup_dep = 9
     3023            ENDIF
     3024         ENDIF
     3025
     3026         IF (surf_lsm_h%frac(ind_wat_win,m) > 0) THEN
     3027            luw_palm = surf_lsm_h%water_type(m)     
     3028            IF (luw_palm == ind_luw_user) THEN
     3029               message_string = 'No water type defined. Please define water type to enable deposition calculation'
     3030               CALL message( 'chem_depo', 'CM0453', 1, 2, 0, 6, 0 )
     3031            ELSEIF (luw_palm ==  ind_luw_lake) THEN
     3032               luw_dep = 13
     3033            ELSEIF (luw_palm == ind_luw_river) THEN
     3034               luw_dep = 13
     3035            ELSEIF (luw_palm == ind_luw_ocean) THEN
     3036               luw_dep = 6
     3037            ELSEIF (luw_palm == ind_luw_pond) THEN
     3038               luw_dep = 13
     3039            ELSEIF (luw_palm == ind_luw_fountain) THEN
     3040               luw_dep = 13
     3041            ENDIF
     3042         ENDIF
     3043
     3044
     3045
     3046         ! Set wetness indicator to dry or wet for lsm vegetation or pavement
     3047         IF (surf_lsm_h%c_liq(m) > 0) THEN
     3048            nwet = 1
     3049         ELSE
     3050            nwet = 0
     3051         ENDIF
     3052
     3053         ! Compute length of time step
     3054         IF ( call_chem_at_all_substeps )  THEN
     3055            dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
     3056         ELSE
     3057            dt_chem = dt_3d
     3058         ENDIF
     3059
     3060
     3061         dh = dzw(k)
     3062         inv_dh = 1.0_wp / dh
     3063         dt_dh = dt_chem/dh
     3064
     3065         ! Concentration at i,j,k
     3066         DO lsp = 1, nspec
     3067            cc(lsp) = chem_species(lsp)%conc(k,j,i)
     3068         ENDDO
     3069
     3070
     3071         ! Temperature column at i,j
     3072         ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
     3073
     3074         ! Surface temperature in degrees Celcius:
     3075         ts       = ttemp - 273.15 ! in degrees celcius
     3076
     3077         ! Viscosity of air in lowest layer
     3078         visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
     3079
     3080         ! Air density in lowest layer
     3081         dens = rho_air_zw(k)
     3082
     3083         ! Calculate relative humidity from specific humidity for DEPAC
     3084         qv_tmp = max(q(k,j,i),0.0_wp)
     3085         relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(k) )
     3086
     3087
     3088
     3089         ! Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
     3090         ! for each surface fraction. Then derive overall budget taking into account the surface fractions.
     3091
     3092         ! Vegetation
     3093         IF (surf_lsm_h%frac(ind_veg_wall,m) > 0) THEN
     3094
     3095
     3096            slinnfac = 1.0_wp
     3097
     3098            ! Get vd
     3099            DO lsp = 1, nvar
     3100               !Initialize
     3101               vs = 0.0_wp
     3102               vd_lu = 0.0_wp
     3103               Rs = 0.0_wp
     3104               Rb = 0.0_wp
     3105               Rc_tot = 0.0_wp
     3106               IF (spc_names(lsp) == 'PM10') THEN
     3107                  part_type = 1
     3108                  ! sedimentation velocity in lowest layer
     3109                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3110                       particle_pars(ind_p_size, part_type), &
     3111                       particle_pars(ind_p_slip, part_type), &
     3112                       visc)
     3113
     3114                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3115                       vs, &
     3116                       particle_pars(ind_p_size, part_type), &
     3117                       particle_pars(ind_p_slip, part_type), &
     3118                       nwet, ttemp, dens, visc, &
     3119                       lu_dep,  &
     3120                       r_aero_lu, ustar_lu)
     3121
     3122                  bud_now_lu(lsp) = - cc(lsp) * &
     3123                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     3124
     3125
     3126               ELSEIF (spc_names(lsp) == 'PM25') THEN
     3127                  part_type = 2
     3128                  ! sedimentation velocity in lowest layer:
     3129                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3130                       particle_pars(ind_p_size, part_type), &
     3131                       particle_pars(ind_p_slip, part_type), &
     3132                       visc)
     3133
     3134
     3135                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3136                       vs, &
     3137                       particle_pars(ind_p_size, part_type), &
     3138                       particle_pars(ind_p_slip, part_type), &
     3139                       nwet, ttemp, dens, visc, &
     3140                       lu_dep , &
     3141                       r_aero_lu, ustar_lu)
     3142
     3143                  bud_now_lu(lsp) = - cc(lsp) * &
     3144                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     3145
     3146
     3147               ELSE
     3148
     3149                  ! GASES
     3150
     3151                  ! Read spc_name of current species for gas parameter
     3152
     3153                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
     3154                     i_pspec = 0
     3155                     DO pspec = 1, 69
     3156                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
     3157                           i_pspec = pspec
     3158                        END IF
     3159                     ENDDO
     3160
     3161                  ELSE
     3162                     ! Species for now not deposited
     3163                     CYCLE
     3164                  ENDIF
     3165
     3166                  ! Factor used for conversion from ppb to ug/m3 :
     3167                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3168                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3169                  !    c           1e-9              xm_tracer         1e9       /       xm_air            dens
     3170                  ! thus:
     3171                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3172                  ! Use density at lowest layer:
     3173
     3174                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
     3175
     3176                  ! atmospheric concentration in DEPAC is requested in ug/m3:
     3177                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     3178                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3179
     3180                  ! Diffusivity for DEPAC relevant gases
     3181                  ! Use default value
     3182                  diffc            = 0.11e-4
     3183                  ! overwrite with known coefficients of diffusivity from Massman (1998)
     3184                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
     3185                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     3186                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
     3187                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
     3188                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
     3189                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
     3190                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
     3191
     3192
     3193                  ! Get quasi-laminar boundary layer resistance Rb:
     3194                  CALL get_rb_cell( (lu_dep == ilu_water_sea) .or. (lu_dep == ilu_water_inland), &
     3195                       z0h_lu, ustar_lu, diffc, &
     3196                       Rb)
     3197
     3198                  ! Get Rc_tot
     3199                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     3200                       relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
     3201                       r_aero_lu , Rb)
     3202
     3203
     3204                  ! Calculate budget
     3205                  IF (Rc_tot .le. 0.0) THEN
     3206
     3207                     bud_now_lu(lsp) = 0.0_wp
     3208
     3209                  ELSE
     3210
     3211                     ! Compute exchange velocity for current lu:
     3212                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
     3213
     3214                     bud_now_lu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3215                          (1.0 - exp(-vd_lu*dt_dh ))*dh
     3216                  ENDIF
     3217
     3218               ENDIF
     3219            ENDDO
     3220         ENDIF
     3221
     3222
     3223         ! Pavement
     3224         IF (surf_lsm_h%frac(ind_pav_green,m) > 0) THEN
     3225
     3226
     3227            ! No vegetation on pavements:
     3228            lai = 0.0_wp
     3229            sai = 0.0_wp
     3230           
     3231            slinnfac = 1.0_wp
     3232
     3233            ! Get vd
     3234            DO lsp = 1, nvar
     3235               !Initialize
     3236               vs = 0.0_wp
     3237               vd_lu = 0.0_wp
     3238               Rs = 0.0_wp
     3239               Rb = 0.0_wp
     3240               Rc_tot = 0.0_wp
     3241               IF (spc_names(lsp) == 'PM10') THEN
     3242                  part_type = 1
     3243                  ! sedimentation velocity in lowest layer:
     3244                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3245                       particle_pars(ind_p_size, part_type), &
     3246                       particle_pars(ind_p_slip, part_type), &
     3247                       visc)
     3248
     3249                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3250                       vs, &
     3251                       particle_pars(ind_p_size, part_type), &
     3252                       particle_pars(ind_p_slip, part_type), &
     3253                       nwet, ttemp, dens, visc, &
     3254                       lup_dep,  &
     3255                       r_aero_lu, ustar_lu)
     3256
     3257                  bud_now_lup(lsp) = - cc(lsp) * &
     3258                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     3259
     3260
     3261               ELSEIF (spc_names(lsp) == 'PM25') THEN
     3262                  part_type = 2
     3263                  ! sedimentation velocity in lowest layer:
     3264                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3265                       particle_pars(ind_p_size, part_type), &
     3266                       particle_pars(ind_p_slip, part_type), &
     3267                       visc)
     3268
     3269
     3270                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3271                       vs, &
     3272                       particle_pars(ind_p_size, part_type), &
     3273                       particle_pars(ind_p_slip, part_type), &
     3274                       nwet, ttemp, dens, visc, &
     3275                       lup_dep, &
     3276                       r_aero_lu, ustar_lu)
     3277
     3278                  bud_now_lup(lsp) = - cc(lsp) * &
     3279                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     3280
     3281
     3282               ELSE
     3283
     3284                  ! GASES
     3285
     3286                  ! Read spc_name of current species for gas parameter
     3287
     3288                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
     3289                     i_pspec = 0
     3290                     DO pspec = 1, 69
     3291                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
     3292                           i_pspec = pspec
     3293                        END IF
     3294                     ENDDO
     3295
     3296                  ELSE
     3297                     ! Species for now not deposited
     3298                     CYCLE
     3299                  ENDIF
     3300
     3301                  ! Factor used for conversion from ppb to ug/m3 :
     3302                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3303                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3304                  !    c           1e-9               xm_tracer         1e9       /       xm_air            dens
     3305                  ! thus:
     3306                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3307                  ! Use density at lowest layer:
     3308
     3309                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
     3310
     3311                  ! atmospheric concentration in DEPAC is requested in ug/m3:
     3312                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     3313                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3314
     3315                  ! Diffusivity for DEPAC relevant gases
     3316                  ! Use default value
     3317                  diffc            = 0.11e-4
     3318                  ! overwrite with known coefficients of diffusivity from Massman (1998)
     3319                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
     3320                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     3321                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
     3322                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
     3323                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
     3324                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
     3325                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
     3326
     3327
     3328                  ! Get quasi-laminar boundary layer resistance Rb:
     3329                  CALL get_rb_cell( (lup_dep == ilu_water_sea) .or. (lup_dep == ilu_water_inland), &
     3330                       z0h_lu, ustar_lu, diffc, &
     3331                       Rb)
     3332
     3333
     3334                  ! Get Rc_tot
     3335                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     3336                       relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
     3337                       r_aero_lu , Rb)
     3338
     3339
     3340                  ! Calculate budget
     3341                  IF (Rc_tot .le. 0.0) THEN
     3342
     3343                     bud_now_lup(lsp) = 0.0_wp
     3344
     3345                  ELSE
     3346
     3347                     ! Compute exchange velocity for current lu:
     3348                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
     3349
     3350                     bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3351                          (1.0 - exp(-vd_lu*dt_dh ))*dh
     3352                  ENDIF
     3353
     3354
     3355               ENDIF
     3356            ENDDO
     3357         ENDIF
     3358
     3359
     3360         ! Water
     3361         IF (surf_lsm_h%frac(ind_wat_win,m) > 0) THEN
     3362
     3363
     3364            ! No vegetation on water:
     3365            lai = 0.0_wp
     3366            sai = 0.0_wp
     3367           
     3368            slinnfac = 1.0_wp
     3369
     3370            ! Get vd
     3371            DO lsp = 1, nvar
     3372               !Initialize
     3373               vs = 0.0_wp
     3374               vd_lu = 0.0_wp
     3375               Rs = 0.0_wp
     3376               Rb = 0.0_wp
     3377               Rc_tot = 0.0_wp
     3378               IF (spc_names(lsp) == 'PM10') THEN
     3379                  part_type = 1
     3380                  ! sedimentation velocity in lowest layer:
     3381                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3382                       particle_pars(ind_p_size, part_type), &
     3383                       particle_pars(ind_p_slip, part_type), &
     3384                       visc)
     3385
     3386                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3387                       vs, &
     3388                       particle_pars(ind_p_size, part_type), &
     3389                       particle_pars(ind_p_slip, part_type), &
     3390                       nwet, ttemp, dens, visc, &
     3391                       luw_dep, &
     3392                       r_aero_lu, ustar_lu)
     3393
     3394                  bud_now_luw(lsp) = - cc(lsp) * &
     3395                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     3396
     3397
     3398               ELSEIF (spc_names(lsp) == 'PM25') THEN
     3399                  part_type = 2
     3400                  ! sedimentation velocity in lowest layer:
     3401                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3402                       particle_pars(ind_p_size, part_type), &
     3403                       particle_pars(ind_p_slip, part_type), &
     3404                       visc)
     3405
     3406
     3407                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3408                       vs, &
     3409                       particle_pars(ind_p_size, part_type), &
     3410                       particle_pars(ind_p_slip, part_type), &
     3411                       nwet, ttemp, dens, visc, &
     3412                       luw_dep, &
     3413                       r_aero_lu, ustar_lu)
     3414
     3415                  bud_now_luw(lsp) = - cc(lsp) * &
     3416                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     3417
     3418
     3419               ELSE
     3420
     3421                  ! GASES
     3422
     3423                  ! Read spc_name of current species for gas PARAMETER
     3424
     3425                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
     3426                     i_pspec = 0
     3427                     DO pspec = 1, 69
     3428                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
     3429                           i_pspec = pspec
     3430                        END IF
     3431                     ENDDO
     3432
     3433                  ELSE
     3434                     ! Species for now not deposited
     3435                     CYCLE
     3436                  ENDIF
     3437
     3438                  ! Factor used for conversion from ppb to ug/m3 :
     3439                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3440                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3441                  !    c           1e-9               xm_tracer         1e9       /       xm_air            dens
     3442                  ! thus:
     3443                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3444                  ! Use density at lowest layer:
     3445
     3446                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
     3447
     3448                  ! atmospheric concentration in DEPAC is requested in ug/m3:
     3449                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     3450                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3451
     3452                  ! Diffusivity for DEPAC relevant gases
     3453                  ! Use default value
     3454                  diffc            = 0.11e-4
     3455                  ! overwrite with known coefficients of diffusivity from Massman (1998)
     3456                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
     3457                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     3458                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
     3459                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
     3460                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
     3461                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
     3462                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
     3463
     3464
     3465                  ! Get quasi-laminar boundary layer resistance Rb:
     3466                  CALL get_rb_cell( (luw_dep == ilu_water_sea) .or. (luw_dep == ilu_water_inland), &
     3467                       z0h_lu, ustar_lu, diffc, &
     3468                       Rb)
     3469
     3470                  ! Get Rc_tot
     3471                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     3472                       relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
     3473                       r_aero_lu , Rb)
     3474
     3475
     3476                  ! Calculate budget
     3477                  IF (Rc_tot .le. 0.0) THEN
     3478
     3479                     bud_now_luw(lsp) = 0.0_wp
     3480
     3481                  ELSE
     3482
     3483                     ! Compute exchange velocity for current lu:
     3484                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
     3485
     3486                     bud_now_luw(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3487                          (1.0 - exp(-vd_lu*dt_dh ))*dh
     3488                  ENDIF
     3489
     3490               ENDIF
     3491            ENDDO
     3492         ENDIF
     3493
     3494
     3495         bud_now = 0.0_wp
     3496         ! Calculate budget for surface m and adapt concentration
     3497         DO lsp = 1, nspec
     3498
     3499
     3500            bud_now(lsp) = surf_lsm_h%frac(ind_veg_wall,m)*bud_now_lu(lsp) + &
     3501                 surf_lsm_h%frac(ind_pav_green,m)*bud_now_lup(lsp) + &
     3502                 surf_lsm_h%frac(ind_wat_win,m)*bud_now_luw(lsp)
     3503
     3504            ! Compute new concentration, add contribution of all exchange processes:
     3505            cc(lsp) = cc(lsp) + bud_now(lsp)*inv_dh
     3506           
     3507            chem_species(lsp)%conc(k,j,i) = max(0.0_wp, cc(lsp))
     3508 
     3509         ENDDO
     3510
     3511      ENDIF
     3512
     3513
     3514
     3515      !!!!!!!
     3516      ! USM !
     3517      !!!!!!!   
     3518
     3519      IF ( match_usm )  THEN
     3520
     3521         ! Get surface element information at i,j:
     3522         m = surf_usm_h%start_index(j,i)
     3523
     3524         k = surf_usm_h%k(m)
     3525
     3526         ! Get needed variables for surface element m
     3527         ustar_lu  = surf_usm_h%us(m)
     3528         z0h_lu    = surf_usm_h%z0h(m)
     3529         r_aero_lu = surf_usm_h%r_a(m)
     3530         glrad     = surf_usm_h%rad_sw_in(m)
     3531         lai = surf_usm_h%lai(m)
     3532         sai = lai + 1
     3533
     3534         ! For small grid spacing neglect R_a
     3535         IF (dzw(k) <= 1.0) THEN
     3536            r_aero_lu = 0.0_wp
     3537         ENDIF
     3538
     3539         !Initialize lu's
     3540         luu_palm = 0
     3541         luu_dep = 0
     3542         lug_palm = 0
     3543         lug_dep = 0
     3544         lud_palm = 0
     3545         lud_dep = 0
     3546
     3547         !Initialize budgets
     3548         bud_now_luu  = 0.0_wp
     3549         bud_now_lug = 0.0_wp
     3550         bud_now_lud = 0.0_wp
     3551
     3552
     3553         ! Get land use for i,j and assign to DEPAC lu
     3554         IF (surf_usm_h%frac(ind_pav_green,m) > 0) THEN
     3555            ! For green urban surfaces (e.g. green roofs
     3556            ! assume LU short grass
     3557            lug_palm = ind_lu_s_grass
     3558            IF (lug_palm == ind_lu_user) THEN
     3559               message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
     3560               CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 )
     3561            ELSEIF (lug_palm == ind_lu_b_soil) THEN
     3562               lug_dep = 9
     3563            ELSEIF (lug_palm == ind_lu_mixed_crops) THEN
     3564               lug_dep = 2
     3565            ELSEIF (lug_palm == ind_lu_s_grass) THEN
     3566               lug_dep = 1
     3567            ELSEIF (lug_palm == ind_lu_ev_needle_trees) THEN
     3568               lug_dep = 4
     3569            ELSEIF (lug_palm == ind_lu_de_needle_trees) THEN
     3570               lug_dep = 4
     3571            ELSEIF (lug_palm == ind_lu_ev_broad_trees) THEN
     3572               lug_dep = 12
     3573            ELSEIF (lug_palm == ind_lu_de_broad_trees) THEN
     3574               lug_dep = 5
     3575            ELSEIF (lug_palm == ind_lu_t_grass) THEN
     3576               lug_dep = 1
     3577            ELSEIF (lug_palm == ind_lu_desert) THEN
     3578               lug_dep = 9
     3579            ELSEIF (lug_palm == ind_lu_tundra ) THEN
     3580               lug_dep = 8
     3581            ELSEIF (lug_palm == ind_lu_irr_crops) THEN
     3582               lug_dep = 2
     3583            ELSEIF (lug_palm == ind_lu_semidesert) THEN
     3584               lug_dep = 8
     3585            ELSEIF (lug_palm == ind_lu_ice) THEN
     3586               lug_dep = 10
     3587            ELSEIF (lug_palm == ind_lu_marsh) THEN
     3588               lug_dep = 8
     3589            ELSEIF (lug_palm == ind_lu_ev_shrubs) THEN
     3590               lug_dep = 14
     3591            ELSEIF (lug_palm == ind_lu_de_shrubs ) THEN
     3592               lug_dep = 14
     3593            ELSEIF (lug_palm == ind_lu_mixed_forest) THEN
     3594               lug_dep = 4
     3595            ELSEIF (lug_palm == ind_lu_intrup_forest) THEN
     3596               lug_dep = 8     
     3597            ENDIF
     3598         ENDIF
     3599
     3600         IF (surf_usm_h%frac(ind_veg_wall,m) > 0) THEN
     3601            ! For walls in USM assume concrete walls/roofs,
     3602            ! assumed LU class desert as also assumed for
     3603            ! pavements in LSM
     3604            luu_palm = ind_lup_conc
     3605            IF (luu_palm == ind_lup_user) THEN
     3606               message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
     3607               CALL message( 'chem_depo', 'CM0455', 1, 2, 0, 6, 0 )
     3608            ELSEIF (luu_palm == ind_lup_asph_conc) THEN
     3609               luu_dep = 9
     3610            ELSEIF (luu_palm == ind_lup_asph) THEN
     3611               luu_dep = 9
     3612            ELSEIF (luu_palm ==  ind_lup_conc) THEN
     3613               luu_dep = 9
     3614            ELSEIF (luu_palm ==  ind_lup_sett) THEN
     3615               luu_dep = 9
     3616            ELSEIF (luu_palm == ind_lup_pav_stones) THEN
     3617               luu_dep = 9
     3618            ELSEIF (luu_palm == ind_lup_cobblest) THEN
     3619               luu_dep = 9       
     3620            ELSEIF (luu_palm == ind_lup_metal) THEN
     3621               luu_dep = 9
     3622            ELSEIF (luu_palm == ind_lup_wood) THEN
     3623               luu_dep = 9   
     3624            ELSEIF (luu_palm == ind_lup_gravel) THEN
     3625               luu_dep = 9
     3626            ELSEIF (luu_palm == ind_lup_f_gravel) THEN
     3627               luu_dep = 9
     3628            ELSEIF (luu_palm == ind_lup_pebblest) THEN
     3629               luu_dep = 9
     3630            ELSEIF (luu_palm == ind_lup_woodchips) THEN
     3631               luu_dep = 9
     3632            ELSEIF (luu_palm == ind_lup_tartan) THEN
     3633               luu_dep = 9
     3634            ELSEIF (luu_palm == ind_lup_art_turf) THEN
     3635               luu_dep = 9
     3636            ELSEIF (luu_palm == ind_lup_clay) THEN
     3637               luu_dep = 9
     3638            ENDIF
     3639         ENDIF
     3640
     3641         IF (surf_usm_h%frac(ind_wat_win,m) > 0) THEN
     3642            ! For windows in USM assume metal as this is
     3643            ! as close as we get, assumed LU class desert
     3644            ! as also assumed for pavements in LSM
     3645            lud_palm = ind_lup_metal     
     3646            IF (lud_palm == ind_lup_user) THEN
     3647               message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
     3648               CALL message( 'chem_depo', 'CM0456', 1, 2, 0, 6, 0 )
     3649            ELSEIF (lud_palm == ind_lup_asph_conc) THEN
     3650               lud_dep = 9
     3651            ELSEIF (lud_palm == ind_lup_asph) THEN
     3652               lud_dep = 9
     3653            ELSEIF (lud_palm ==  ind_lup_conc) THEN
     3654               lud_dep = 9
     3655            ELSEIF (lud_palm ==  ind_lup_sett) THEN
     3656               lud_dep = 9
     3657            ELSEIF (lud_palm == ind_lup_pav_stones) THEN
     3658               lud_dep = 9
     3659            ELSEIF (lud_palm == ind_lup_cobblest) THEN
     3660               lud_dep = 9       
     3661            ELSEIF (lud_palm == ind_lup_metal) THEN
     3662               lud_dep = 9
     3663            ELSEIF (lud_palm == ind_lup_wood) THEN
     3664               lud_dep = 9   
     3665            ELSEIF (lud_palm == ind_lup_gravel) THEN
     3666               lud_dep = 9
     3667            ELSEIF (lud_palm == ind_lup_f_gravel) THEN
     3668               lud_dep = 9
     3669            ELSEIF (lud_palm == ind_lup_pebblest) THEN
     3670               lud_dep = 9
     3671            ELSEIF (lud_palm == ind_lup_woodchips) THEN
     3672               lud_dep = 9
     3673            ELSEIF (lud_palm == ind_lup_tartan) THEN
     3674               lud_dep = 9
     3675            ELSEIF (lud_palm == ind_lup_art_turf) THEN
     3676               lud_dep = 9
     3677            ELSEIF (lud_palm == ind_lup_clay) THEN
     3678               lud_dep = 9
     3679            ENDIF
     3680         ENDIF
     3681
     3682
     3683         ! TODO: Activate these lines as soon as new ebsolver branch is merged:
     3684         ! Set wetness indicator to dry or wet for usm vegetation or pavement
     3685         !IF (surf_usm_h%c_liq(m) > 0) THEN
     3686         !   nwet = 1
     3687         !ELSE
     3688         nwet = 0
     3689         !ENDIF
     3690
     3691         ! Compute length of time step
     3692         IF ( call_chem_at_all_substeps )  THEN
     3693            dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
     3694         ELSE
     3695            dt_chem = dt_3d
     3696         ENDIF
     3697
     3698
     3699         dh = dzw(k)
     3700         inv_dh = 1.0_wp / dh
     3701         dt_dh = dt_chem/dh
     3702
     3703         ! Concentration at i,j,k
     3704         DO lsp = 1, nspec
     3705            cc(lsp) = chem_species(lsp)%conc(k,j,i)
     3706         ENDDO
     3707
     3708         ! Temperature at i,j,k
     3709         ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
     3710
     3711         ! Surface temperature in degrees Celcius:
     3712         ts       = ttemp - 273.15 ! in degrees celcius
     3713
     3714         ! Viscosity of air in lowest layer
     3715         visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
     3716
     3717         ! Air density in lowest layer
     3718         dens = rho_air_zw(k)
     3719
     3720         ! Calculate relative humidity from specific humidity for DEPAC
     3721         qv_tmp = max(q(k,j,i),0.0_wp)
     3722         relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(nzb) )
     3723
     3724
     3725
     3726         ! Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
     3727         ! for each surface fraction. Then derive overall budget taking into account the surface fractions.
     3728
     3729         ! Walls/roofs
     3730         IF (surf_usm_h%frac(ind_veg_wall,m) > 0) THEN
     3731
     3732
     3733            ! No vegetation on non-green walls:
     3734            lai = 0.0_wp
     3735            sai = 0.0_wp
     3736           
     3737            slinnfac = 1.0_wp
     3738
     3739            ! Get vd
     3740            DO lsp = 1, nvar
     3741               !Initialize
     3742               vs = 0.0_wp
     3743               vd_lu = 0.0_wp
     3744               Rs = 0.0_wp
     3745               Rb = 0.0_wp
     3746               Rc_tot = 0.0_wp
     3747               IF (spc_names(lsp) == 'PM10') THEN
     3748                  part_type = 1
     3749                  ! sedimentation velocity in lowest layer
     3750                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3751                       particle_pars(ind_p_size, part_type), &
     3752                       particle_pars(ind_p_slip, part_type), &
     3753                       visc)
     3754
     3755                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3756                       vs, &
     3757                       particle_pars(ind_p_size, part_type), &
     3758                       particle_pars(ind_p_slip, part_type), &
     3759                       nwet, ttemp, dens, visc, &
     3760                       luu_dep,  &
     3761                       r_aero_lu, ustar_lu)
     3762
     3763                  bud_now_luu(lsp) = - cc(lsp) * &
     3764                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     3765
     3766
     3767               ELSEIF (spc_names(lsp) == 'PM25') THEN
     3768                  part_type = 2
     3769                  ! sedimentation velocity in lowest layer:
     3770                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3771                       particle_pars(ind_p_size, part_type), &
     3772                       particle_pars(ind_p_slip, part_type), &
     3773                       visc)
     3774
     3775
     3776                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3777                       vs, &
     3778                       particle_pars(ind_p_size, part_type), &
     3779                       particle_pars(ind_p_slip, part_type), &
     3780                       nwet, ttemp, dens, visc, &
     3781                       luu_dep , &
     3782                       r_aero_lu, ustar_lu)
     3783
     3784                  bud_now_luu(lsp) = - cc(lsp) * &
     3785                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     3786
     3787
     3788               ELSE
     3789
     3790                  ! GASES
     3791
     3792                  ! Read spc_name of current species for gas parameter
     3793
     3794                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
     3795                     i_pspec = 0
     3796                     DO pspec = 1, 69
     3797                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
     3798                           i_pspec = pspec
     3799                        END IF
     3800                     ENDDO
     3801
     3802                  ELSE
     3803                     ! Species for now not deposited
     3804                     CYCLE
     3805                  ENDIF
     3806
     3807                  ! Factor used for conversion from ppb to ug/m3 :
     3808                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3809                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3810                  !    c           1e-9              xm_tracer         1e9       /       xm_air            dens
     3811                  ! thus:
     3812                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3813                  ! Use density at lowest layer:
     3814
     3815                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
     3816
     3817                  ! atmospheric concentration in DEPAC is requested in ug/m3:
     3818                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     3819                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3820
     3821                  ! Diffusivity for DEPAC relevant gases
     3822                  ! Use default value
     3823                  diffc            = 0.11e-4
     3824                  ! overwrite with known coefficients of diffusivity from Massman (1998)
     3825                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
     3826                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     3827                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
     3828                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
     3829                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
     3830                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
     3831                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
     3832
     3833
     3834                  ! Get quasi-laminar boundary layer resistance Rb:
     3835                  CALL get_rb_cell( (luu_dep == ilu_water_sea) .or. (luu_dep == ilu_water_inland), &
     3836                       z0h_lu, ustar_lu, diffc, &
     3837                       Rb)
     3838
     3839                  ! Get Rc_tot
     3840                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     3841                       relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
     3842                       r_aero_lu , Rb)
     3843
     3844
     3845                  ! Calculate budget
     3846                  IF (Rc_tot .le. 0.0) THEN
     3847
     3848                     bud_now_luu(lsp) = 0.0_wp
     3849
     3850                  ELSE
     3851
     3852                     ! Compute exchange velocity for current lu:
     3853                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
     3854
     3855                     bud_now_luu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3856                          (1.0 - exp(-vd_lu*dt_dh ))*dh
     3857                  ENDIF
     3858
     3859               ENDIF
     3860            ENDDO
     3861         ENDIF
     3862
     3863
     3864         ! Green usm surfaces
     3865         IF (surf_usm_h%frac(ind_pav_green,m) > 0) THEN
     3866
     3867
     3868            slinnfac = 1.0_wp
     3869
     3870            ! Get vd
     3871            DO lsp = 1, nvar
     3872               !Initialize
     3873               vs = 0.0_wp
     3874               vd_lu = 0.0_wp
     3875               Rs = 0.0_wp
     3876               Rb = 0.0_wp
     3877               Rc_tot = 0.0_wp
     3878               IF (spc_names(lsp) == 'PM10') THEN
     3879                  part_type = 1
     3880                  ! sedimentation velocity in lowest layer:
     3881                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3882                       particle_pars(ind_p_size, part_type), &
     3883                       particle_pars(ind_p_slip, part_type), &
     3884                       visc)
     3885
     3886                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3887                       vs, &
     3888                       particle_pars(ind_p_size, part_type), &
     3889                       particle_pars(ind_p_slip, part_type), &
     3890                       nwet, ttemp, dens, visc, &
     3891                       lug_dep,  &
     3892                       r_aero_lu, ustar_lu)
     3893
     3894                  bud_now_lug(lsp) = - cc(lsp) * &
     3895                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     3896
     3897
     3898               ELSEIF (spc_names(lsp) == 'PM25') THEN
     3899                  part_type = 2
     3900                  ! sedimentation velocity in lowest layer:
     3901                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3902                       particle_pars(ind_p_size, part_type), &
     3903                       particle_pars(ind_p_slip, part_type), &
     3904                       visc)
     3905
     3906
     3907                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3908                       vs, &
     3909                       particle_pars(ind_p_size, part_type), &
     3910                       particle_pars(ind_p_slip, part_type), &
     3911                       nwet, ttemp, dens, visc, &
     3912                       lug_dep, &
     3913                       r_aero_lu, ustar_lu)
     3914
     3915                  bud_now_lug(lsp) = - cc(lsp) * &
     3916                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     3917
     3918
     3919               ELSE
     3920
     3921                  ! GASES
     3922
     3923                  ! Read spc_name of current species for gas parameter
     3924
     3925                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
     3926                     i_pspec = 0
     3927                     DO pspec = 1, 69
     3928                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
     3929                           i_pspec = pspec
     3930                        END IF
     3931                     ENDDO
     3932
     3933                  ELSE
     3934                     ! Species for now not deposited
     3935                     CYCLE
     3936                  ENDIF
     3937
     3938                  ! Factor used for conversion from ppb to ug/m3 :
     3939                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3940                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3941                  !    c           1e-9               xm_tracer         1e9       /       xm_air            dens
     3942                  ! thus:
     3943                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3944                  ! Use density at lowest layer:
     3945
     3946                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
     3947
     3948                  ! atmospheric concentration in DEPAC is requested in ug/m3:
     3949                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     3950                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3951
     3952                  ! Diffusivity for DEPAC relevant gases
     3953                  ! Use default value
     3954                  diffc            = 0.11e-4
     3955                  ! overwrite with known coefficients of diffusivity from Massman (1998)
     3956                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
     3957                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     3958                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
     3959                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
     3960                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
     3961                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
     3962                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
     3963
     3964
     3965                  ! Get quasi-laminar boundary layer resistance Rb:
     3966                  CALL get_rb_cell( (lug_dep == ilu_water_sea) .or. (lug_dep == ilu_water_inland), &
     3967                       z0h_lu, ustar_lu, diffc, &
     3968                       Rb)
     3969
     3970
     3971                  ! Get Rc_tot
     3972                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     3973                       relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
     3974                       r_aero_lu , Rb)
     3975
     3976
     3977                  ! Calculate budget
     3978                  IF (Rc_tot .le. 0.0) THEN
     3979
     3980                     bud_now_lug(lsp) = 0.0_wp
     3981
     3982                  ELSE
     3983
     3984                     ! Compute exchange velocity for current lu:
     3985                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
     3986
     3987                     bud_now_lug(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3988                          (1.0 - exp(-vd_lu*dt_dh ))*dh
     3989                  ENDIF
     3990
     3991
     3992               ENDIF
     3993            ENDDO
     3994         ENDIF
     3995
     3996
     3997         ! Windows
     3998         IF (surf_usm_h%frac(ind_wat_win,m) > 0) THEN
     3999
     4000
     4001            ! No vegetation on windows:
     4002            lai = 0.0_wp
     4003            sai = 0.0_wp
     4004           
     4005            slinnfac = 1.0_wp
     4006
     4007            ! Get vd
     4008            DO lsp = 1, nvar
     4009               !Initialize
     4010               vs = 0.0_wp
     4011               vd_lu = 0.0_wp
     4012               Rs = 0.0_wp
     4013               Rb = 0.0_wp
     4014               Rc_tot = 0.0_wp
     4015               IF (spc_names(lsp) == 'PM10') THEN
     4016                  part_type = 1
     4017                  ! sedimentation velocity in lowest layer:
     4018                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     4019                       particle_pars(ind_p_size, part_type), &
     4020                       particle_pars(ind_p_slip, part_type), &
     4021                       visc)
     4022
     4023                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     4024                       vs, &
     4025                       particle_pars(ind_p_size, part_type), &
     4026                       particle_pars(ind_p_slip, part_type), &
     4027                       nwet, ttemp, dens, visc, &
     4028                       lud_dep, &
     4029                       r_aero_lu, ustar_lu)
     4030
     4031                  bud_now_lud(lsp) = - cc(lsp) * &
     4032                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     4033
     4034
     4035               ELSEIF (spc_names(lsp) == 'PM25') THEN
     4036                  part_type = 2
     4037                  ! sedimentation velocity in lowest layer:
     4038                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     4039                       particle_pars(ind_p_size, part_type), &
     4040                       particle_pars(ind_p_slip, part_type), &
     4041                       visc)
     4042
     4043
     4044                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     4045                       vs, &
     4046                       particle_pars(ind_p_size, part_type), &
     4047                       particle_pars(ind_p_slip, part_type), &
     4048                       nwet, ttemp, dens, visc, &
     4049                       lud_dep, &
     4050                       r_aero_lu, ustar_lu)
     4051
     4052                  bud_now_lud(lsp) = - cc(lsp) * &
     4053                       (1.0 - exp(-vd_lu*dt_dh ))*dh
     4054
     4055
     4056               ELSE
     4057
     4058                  ! GASES
     4059
     4060                  ! Read spc_name of current species for gas PARAMETER
     4061
     4062                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
     4063                     i_pspec = 0
     4064                     DO pspec = 1, 69
     4065                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
     4066                           i_pspec = pspec
     4067                        END IF
     4068                     ENDDO
     4069
     4070                  ELSE
     4071                     ! Species for now not deposited
     4072                     CYCLE
     4073                  ENDIF
     4074
     4075                  ! Factor used for conversion from ppb to ug/m3 :
     4076                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     4077                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     4078                  !    c           1e-9               xm_tracer         1e9       /       xm_air            dens
     4079                  ! thus:
     4080                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     4081                  ! Use density at lowest layer:
     4082
     4083                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
     4084
     4085                  ! atmospheric concentration in DEPAC is requested in ug/m3:
     4086                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     4087                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     4088
     4089                  ! Diffusivity for DEPAC relevant gases
     4090                  ! Use default value
     4091                  diffc            = 0.11e-4
     4092                  ! overwrite with known coefficients of diffusivity from Massman (1998)
     4093                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
     4094                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     4095                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
     4096                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
     4097                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
     4098                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
     4099                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
     4100
     4101
     4102                  ! Get quasi-laminar boundary layer resistance Rb:
     4103                  CALL get_rb_cell( (lud_dep == ilu_water_sea) .or. (lud_dep == ilu_water_inland), &
     4104                       z0h_lu, ustar_lu, diffc, &
     4105                       Rb)
     4106
     4107                  ! Get Rc_tot
     4108                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     4109                       relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
     4110                       r_aero_lu , Rb)
     4111
     4112
     4113                  ! Calculate budget
     4114                  IF (Rc_tot .le. 0.0) THEN
     4115
     4116                     bud_now_lud(lsp) = 0.0_wp
     4117
     4118                  ELSE
     4119
     4120                     ! Compute exchange velocity for current lu:
     4121                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
     4122
     4123                     bud_now_lud(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     4124                          (1.0 - exp(-vd_lu*dt_dh ))*dh
     4125                  ENDIF
     4126
     4127               ENDIF
     4128            ENDDO
     4129         ENDIF
     4130
     4131
     4132         bud_now = 0.0_wp
     4133         ! Calculate budget for surface m and adapt concentration
     4134         DO lsp = 1, nspec
     4135
     4136
     4137            bud_now(lsp) = surf_usm_h%frac(ind_veg_wall,m)*bud_now_luu(lsp) + &
     4138                 surf_usm_h%frac(ind_pav_green,m)*bud_now_lug(lsp) + &
     4139                 surf_usm_h%frac(ind_wat_win,m)*bud_now_lud(lsp)
     4140
     4141            ! Compute new concentration, add contribution of all exchange processes:
     4142            cc(lsp) = cc(lsp) + bud_now(lsp)*inv_dh
     4143
     4144            chem_species(lsp)%conc(k,j,i) = max(0.0_wp, cc(lsp))
     4145
     4146         ENDDO
     4147
     4148      ENDIF
     4149
     4150
     4151
     4152    END SUBROUTINE chem_depo
     4153
     4154 
     4155
     4156!----------------------------------------------------------------------------------
     4157!
     4158! DEPAC:
     4159! Code of the DEPAC routine and corresponding subroutines below from the DEPAC
     4160! module of the LOTOS-EUROS model (Manders et al., 2017)
     4161!   
     4162! Original DEPAC routines by RIVM and TNO (2015), for Documentation see
     4163! van Zanten et al., 2010.
     4164!---------------------------------------------------------------------------------
     4165!   
     4166!----------------------------------------------------------------------------------   
     4167!
     4168! depac: compute total canopy (or surface) resistance Rc for gases
     4169!----------------------------------------------------------------------------------
     4170
     4171    SUBROUTINE drydepos_gas_depac(compnam, day_of_year, lat, t, ust, glrad, sinphi, &
     4172         rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, catm, diffc, &
     4173         ra, rb) 
     4174
     4175
     4176      ! The last four rows of depac arguments are OPTIONAL:
     4177      !
     4178      ! A. compute Rc_tot without compensation points (ccomp_tot will be zero):
     4179      !     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi])
     4180      !
     4181      ! B. compute Rc_tot with compensation points (used for LOTOS-EUROS):
     4182      !     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
     4183      !                  c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
     4184      !
     4185      ! C. compute effective Rc based on compensation points (used for OPS):
     4186      !     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
     4187      !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
     4188      !                 ra, rb, rc_eff)
     4189      ! X1. Extra (OPTIONAL) output variables:
     4190      !     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
     4191      !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
     4192      !                 ra, rb, rc_eff, &
     4193      !                 gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out)
     4194      ! X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves):
     4195      !     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
     4196      !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
     4197      !                 ra, rb, rc_eff, &
     4198      !                 gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, &
     4199      !                 calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
     4200
     4201
     4202      IMPLICIT NONE
     4203
     4204      CHARACTER(len=*)          , INTENT(in)       :: compnam          ! component name
     4205      ! 'HNO3','NO','NO2','O3','SO2','NH3'
     4206      INTEGER(iwp)              , INTENT(in)       :: day_of_year      ! day of year, 1 ... 365 (366)
     4207      REAL(kind=wp)             , INTENT(in)       :: lat                   ! latitude Northern hemisphere (degrees) (DEPAC cannot be used for S. hemisphere)
     4208      REAL(kind=wp)            , INTENT(in)        :: t                ! temperature (C)
     4209      ! NB discussion issue is temp T_2m or T_surf or T_leaf?
     4210      REAL(kind=wp)             , INTENT(in)       :: ust              ! friction velocity (m/s)
     4211      REAL(kind=wp)             , INTENT(in)       :: glrad            ! global radiation (W/m2)
     4212      REAL(kind=wp)             , INTENT(in)       :: sinphi           ! sin of solar elevation angle
     4213      REAL(kind=wp)             , INTENT(in)       :: rh               ! relative humidity (%)
     4214      REAL(kind=wp)             , INTENT(in)       :: lai              ! one-sidedleaf area index (-)
     4215      REAL(kind=wp)             , INTENT(in)       :: sai              ! surface area index (-) (lai + branches and stems)
     4216      REAL(kind=wp)             , INTENT(in)       :: diffc            ! Diffusivity
     4217      INTEGER(iwp)              , INTENT(in)           :: nwet             ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     4218      INTEGER(iwp)              , INTENT(in)       :: lu               ! land use type, lu = 1,...,nlu
     4219      INTEGER(iwp)              , INTENT(in)       :: iratns           ! index for NH3/SO2 ratio used for SO2;
     4220      ! iratns = 1: low NH3/SO2
     4221      ! iratns = 2: high NH3/SO2
     4222      ! iratns = 3: very low NH3/SO2
     4223      REAL(kind=wp)             , INTENT(out)      :: rc_tot           ! total canopy resistance Rc (s/m)
     4224      REAL(kind=wp)             , INTENT(out)      :: ccomp_tot        ! total compensation point (ug/m3) [= 0 for species that don't have a compensation
     4225      REAL(kind=wp)             ,INTENT(in)        :: p                ! pressure (Pa)
     4226      REAL(kind=wp), INTENT(in)  :: catm            ! actual atmospheric concentration (ug/m3)
     4227      REAL(kind=wp), INTENT(in)  :: ra              ! aerodynamic resistance (s/m)
     4228      REAL(kind=wp), INTENT(in)  :: rb              ! boundary layer resistance (s/m)
     4229
     4230      ! Local variables:
     4231      REAL(kind=wp)                           :: laimax          ! maximum leaf area index (-)
     4232      LOGICAL                        :: ready           ! Rc has been set
     4233      ! = 1 -> constant Rc
     4234      ! = 2 -> temperature dependent Rc
     4235      REAL(kind=wp)                           :: gw              ! external leaf conductance (m/s)
     4236      REAL(kind=wp)                           :: gstom           ! stomatal conductance (m/s)
     4237      REAL(kind=wp)                           :: gsoil_eff       ! effective soil conductance (m/s)
     4238      REAL(kind=wp)                           :: gc_tot          ! total canopy conductance (m/s)
     4239      REAL(kind=wp)                           :: cw              ! external leaf surface compensation point (ug/m3)
     4240      REAL(kind=wp)                           :: cstom           ! stomatal compensation point (ug/m3)
     4241      REAL(kind=wp)                           :: csoil           ! soil compensation point (ug/m3)
     4242
     4243      ! Vegetation indicators:
     4244      LOGICAL            :: LAI_present ! leaves are present for current land use type
     4245      LOGICAL            :: SAI_present ! vegetation is present for current land use type
     4246
     4247      ! Component number taken from component name, paramteres matched with include files
     4248      INTEGER(iwp)           ::  icmp
     4249
     4250      ! component numbers:
     4251      INTEGER(iwp), PARAMETER  ::  icmp_o3   = 1
     4252      INTEGER(iwp), PARAMETER  ::  icmp_so2  = 2
     4253      INTEGER(iwp), PARAMETER  ::  icmp_no2  = 3
     4254      INTEGER(iwp), PARAMETER  ::  icmp_no   = 4
     4255      INTEGER(iwp), PARAMETER  ::  icmp_nh3  = 5
     4256      INTEGER(iwp), PARAMETER  ::  icmp_co   = 6
     4257      INTEGER(iwp), PARAMETER  ::  icmp_no3  = 7
     4258      INTEGER(iwp), PARAMETER  ::  icmp_hno3 = 8
     4259      INTEGER(iwp), PARAMETER  ::  icmp_n2o5 = 9
     4260      INTEGER(iwp), PARAMETER  ::  icmp_h2o2 = 10
     4261
     4262
     4263
     4264      ! Define component number
     4265      SELECT CASE ( trim(compnam) )
     4266
     4267      CASE ( 'O3', 'o3' )
     4268         icmp = icmp_o3
     4269
     4270      CASE ( 'SO2', 'so2' )
     4271         icmp = icmp_so2
     4272
     4273      CASE ( 'NO2', 'no2' )
     4274         icmp = icmp_no2
     4275
     4276      CASE ( 'NO', 'no' )
     4277         icmp = icmp_no
     4278
     4279      CASE ( 'NH3', 'nh3' )
     4280         icmp = icmp_nh3
     4281
     4282      CASE ( 'CO', 'co' )
     4283         icmp = icmp_co
     4284
     4285      CASE ( 'NO3', 'no3' )
     4286         icmp = icmp_no3
     4287
     4288      CASE ( 'HNO3', 'hno3' )
     4289         icmp = icmp_hno3
     4290
     4291      CASE ( 'N2O5', 'n2o5' )
     4292         icmp = icmp_n2o5
     4293
     4294      CASE ( 'H2O2', 'h2o2' )
     4295         icmp = icmp_h2o2
     4296
     4297      CASE default
     4298         !Component not part of DEPAC --> not deposited
     4299         RETURN
     4300
     4301      END SELECT
     4302
     4303      ! inititalize
     4304      gw        = 0.
     4305      gstom     = 0.
     4306      gsoil_eff = 0.
     4307      gc_tot    = 0.
     4308      cw        = 0.
     4309      cstom     = 0.
     4310      csoil     = 0.
     4311
     4312
     4313      ! Check whether vegetation is present (in that CASE the leaf or surface area index > 0):
     4314      LAI_present = (lai .gt. 0.0)
     4315      SAI_present = (sai .gt. 0.0)
     4316
     4317      ! Set Rc (i.e. rc_tot) in special cases:
     4318      CALL rc_special(icmp,compnam,lu,t, nwet,rc_tot,ready,ccomp_tot)
     4319
     4320      ! set effective resistance equal to total resistance, this will be changed for compensation points
     4321      !IF ( present(rc_eff) ) then
     4322      !   rc_eff = rc_tot
     4323      !END IF
     4324
     4325      ! IF Rc is not set:
     4326      IF (.not. ready) then
     4327
     4328         ! External conductance:
     4329         CALL rc_gw(compnam,iratns,t,rh,nwet,SAI_present,sai,gw)         
     4330
     4331         ! Stomatal conductance:
     4332         ! CALL rc_gstom(icmp,compnam,lu,LAI_present,lai,glrad,sinphi,t,rh,diffc,gstom,p,&
     4333         !     smi=smi,calc_stom_o3flux=calc_stom_o3flux)
     4334         CALL rc_gstom(icmp,compnam,lu,LAI_present,lai,glrad,sinphi,t,rh,diffc,gstom,p)
     4335         ! Effective soil conductance:
     4336         CALL rc_gsoil_eff(icmp,lu,sai,ust,nwet,t,gsoil_eff)
     4337
     4338         ! Total canopy conductance (gc_tot) and resistance Rc (rc_tot):
     4339         CALL rc_rctot(gstom,gsoil_eff,gw,gc_tot,rc_tot)
     4340
     4341         ! Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3):
     4342         ! CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,&
     4343         !     LAI_present, SAI_present, &
     4344         !     ccomp_tot, &
     4345         !     catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, &
     4346         !     c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, &
     4347         !     tsea=tsea,cw=cw,cstom=cstom,csoil=csoil )
     4348
     4349         ! Effective Rc based on compensation points:
     4350         ! IF ( present(rc_eff) ) then
     4351         ! check on required arguments:
     4352         !IF ( (.not. present(catm)) .or. (.not. present(ra)) .or. (.not. present(rb)) ) then
     4353         !   stop 'output argument rc_eff requires input arguments catm, ra and rb'
     4354         !END IF
     4355         ! compute rc_eff :
     4356         !CALL rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff)
     4357         !ENDIF
     4358      ENDIF
     4359
     4360    END SUBROUTINE drydepos_gas_depac
     4361
     4362
     4363
     4364!-------------------------------------------------------------------
     4365! rc_special: compute total canopy resistance in special CASEs
     4366!-------------------------------------------------------------------
     4367    SUBROUTINE rc_special(icmp,compnam,lu,t,nwet,rc_tot,ready,ccomp_tot)
     4368
     4369      INTEGER(iwp)         , INTENT(in)  :: icmp            ! component index
     4370      CHARACTER(len=*)     , INTENT(in)  :: compnam         ! component name
     4371      INTEGER(iwp)         , INTENT(in)  :: lu              ! land use type, lu = 1,...,nlu
     4372      REAL(kind=wp)        , INTENT(in)  :: t               ! temperature (C)
     4373      ! = 1 -> constant Rc
     4374      ! = 2 -> temperature dependent Rc
     4375      INTEGER(iwp)         , INTENT(in)  :: nwet            ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     4376      REAL(kind=wp)        , INTENT(out) :: rc_tot          ! total canopy resistance Rc (s/m)
     4377      LOGICAL              , INTENT(out) :: ready           ! Rc has been set
     4378      REAL(kind=wp)        , INTENT(out) :: ccomp_tot       ! total compensation point (ug/m3)
     4379
     4380      ! rc_tot is not yet set:
     4381      ready = .false.
     4382
     4383      ! Default compensation point in special CASEs = 0:
     4384      ccomp_tot = 0.0
     4385
     4386      SELECT CASE(trim(compnam))
     4387      CASE('HNO3','N2O5','NO3','H2O2')
     4388         ! No separate resistances for HNO3; just one total canopy resistance:
     4389         IF (t .lt. -5.0 .and. nwet .eq. 9) then
     4390            ! T < 5 C and snow:
     4391            rc_tot = 50.
     4392         ELSE
     4393            ! all other circumstances:
     4394            rc_tot = 10.0
     4395         ENDIF
     4396         ready = .true.
     4397
     4398      CASE('NO','CO')
     4399         IF (lu .eq. ilu_water_sea .or. lu .eq. ilu_water_inland) then       ! water
     4400            rc_tot = 2000.
     4401            ready = .true.
     4402         ELSEIF (nwet .eq. 1) then ! wet
     4403            rc_tot = 2000.
     4404            ready = .true.
     4405         ENDIF
     4406      CASE('NO2','O3','SO2','NH3')
     4407         ! snow surface:
     4408         IF (nwet.eq.9) then
     4409            !CALL rc_snow(ipar_snow(icmp),t,rc_tot)
     4410            ready = .true.
     4411         ENDIF
     4412      CASE default
     4413         message_string = 'Component '// trim(compnam) // ' not supported'
     4414         CALL message( 'rc_special', 'CM0457', 1, 2, 0, 6, 0 )
     4415      END SELECT
     4416
     4417    END SUBROUTINE rc_special
     4418
     4419
     4420
     4421!-------------------------------------------------------------------
     4422! rc_gw: compute external conductance
     4423!-------------------------------------------------------------------
     4424    SUBROUTINE rc_gw( compnam, iratns,t,rh,nwet, SAI_present, sai, gw )
     4425
     4426      ! Input/output variables:
     4427      CHARACTER(len=*)     , INTENT(in)  :: compnam ! component name
     4428      INTEGER(iwp)         , INTENT(in)  :: iratns  ! index for NH3/SO2 ratio;
     4429      ! iratns = 1: low NH3/SO2
     4430      ! iratns = 2: high NH3/SO2
     4431      ! iratns = 3: very low NH3/SO2
     4432      REAL(kind=wp)         , INTENT(in)  :: t       ! temperature (C)
     4433      REAL(kind=wp)         , INTENT(in)  :: rh      ! relative humidity (%)
     4434      INTEGER(iwp)          , INTENT(in)  :: nwet    ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     4435      LOGICAL               , INTENT(in)  :: SAI_present
     4436      REAL(kind=wp)         , INTENT(in)  :: sai     ! one-sided leaf area index (-)
     4437      REAL(kind=wp)         , INTENT(out) :: gw      ! external leaf conductance (m/s)
     4438
     4439      SELECT CASE(trim(compnam))
     4440
     4441         !CASE('HNO3','N2O5','NO3','H2O2') this routine is not CALLed for HNO3
     4442
     4443      CASE('NO2')
     4444         CALL rw_constant(2000.0_wp,SAI_present,gw)
     4445
     4446      CASE('NO','CO')
     4447         CALL rw_constant(-9999.0_wp,SAI_present,gw)  ! see Erisman et al, 1994 section 3.2.3
     4448
     4449      CASE('O3')
     4450         CALL rw_constant(2500.0_wp,SAI_present,gw)
     4451
     4452      CASE('SO2')
     4453         CALL rw_so2( t, nwet, rh, iratns, SAI_present, gw )
     4454
     4455      CASE('NH3')
     4456         CALL rw_nh3_sutton(t,rh,SAI_present,gw)
     4457
     4458         ! conversion from leaf resistance to canopy resistance by multiplying with SAI:
     4459         Gw = sai*gw
     4460
     4461      CASE default
     4462         message_string = 'Component '// trim(compnam) // ' not supported'
     4463         CALL message( 'rc_gw', 'CM0458', 1, 2, 0, 6, 0 )
     4464      END SELECT
     4465
     4466    END SUBROUTINE rc_gw
     4467
     4468
     4469
     4470!-------------------------------------------------------------------
     4471! rw_so2: compute external leaf conductance for SO2
     4472!-------------------------------------------------------------------
     4473    SUBROUTINE rw_so2( t, nwet, rh, iratns, SAI_present, gw )
     4474
     4475      ! Input/output variables:
     4476      REAL(kind=wp)   , INTENT(in)  :: t      ! temperature (C)
     4477      INTEGER(iwp)    , INTENT(in)  :: nwet   ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     4478      REAL(kind=wp)   , INTENT(in)  :: rh     ! relative humidity (%)
     4479      INTEGER(iwp)    , INTENT(in)  :: iratns ! index for NH3/SO2 ratio;
     4480      ! iratns = 1: low NH3/SO2
     4481      ! iratns = 2: high NH3/SO2
     4482      ! iratns = 3: very low NH3/SO2
     4483      LOGICAL, INTENT(in)  :: SAI_present
     4484      REAL(kind=wp)   , INTENT(out) :: gw     ! external leaf conductance (m/s)
     4485
     4486      ! Variables from module:
     4487      ! SAI_present: vegetation is present
     4488
     4489      ! Local variables:
     4490      REAL(kind=wp)                 :: rw     ! external leaf resistance (s/m)
     4491
     4492      ! Check IF vegetation present:
     4493      IF (SAI_present) then
     4494
     4495         IF (nwet .eq. 0) then
     4496            !--------------------------
     4497            ! dry surface
     4498            !--------------------------
     4499            ! T > -1 C
     4500            IF (t .gt. -1.0) then
     4501               IF (rh .lt. 81.3) then
     4502                  rw = 25000*exp(-0.0693*rh)
     4503               ELSE
     4504                  rw = 0.58e12*exp(-0.278*rh) + 10.
     4505               ENDIF
     4506            ELSE
     4507               ! -5 C < T <= -1 C
     4508               IF (t .gt. -5.0) then
     4509                  rw=200
     4510               ELSE
     4511                  ! T <= -5 C
     4512                  rw=500
     4513               ENDIF
     4514            ENDIF
     4515         ELSE
     4516            !--------------------------
     4517            ! wet surface
     4518            !--------------------------
     4519            rw = 10. !see Table 5, Erisman et al, 1994 Atm. Environment, 0 is impl. as 10
     4520         ENDIF
     4521
     4522         ! very low NH3/SO2 ratio:
     4523         IF (iratns == iratns_very_low ) rw = rw+50.
     4524
     4525         ! Conductance:
     4526         gw = 1./rw
     4527      ELSE
     4528         ! no vegetation:
     4529         gw = 0.0
     4530      ENDIF
     4531
     4532    END SUBROUTINE rw_so2
     4533
     4534
     4535
     4536!-------------------------------------------------------------------
     4537! rw_nh3_sutton: compute external leaf conductance for NH3,
     4538!                  following Sutton & Fowler, 1993
     4539!-------------------------------------------------------------------
     4540    SUBROUTINE rw_nh3_sutton(tsurf,rh,SAI_present,gw)
     4541
     4542      ! Input/output variables:
     4543      REAL(kind=wp)   , INTENT(in)  :: tsurf     ! surface temperature (C)
     4544      REAL(kind=wp)   , INTENT(in)  :: rh        ! relative humidity (%)
     4545      LOGICAL, INTENT(in)  :: SAI_present
     4546      REAL(kind=wp)   , INTENT(out) :: gw        ! external leaf conductance (m/s)
     4547
     4548      ! Variables from module:
     4549      ! SAI_present: vegetation is present
     4550
     4551      ! Local variables:
     4552      REAL(kind=wp)                 :: rw                ! external leaf resistance (s/m)
     4553      REAL(kind=wp)                 :: sai_grass_haarweg ! surface area index at experimental site Haarweg
     4554
     4555      ! Fix SAI_grass at value valid for Haarweg data for which gamma_w parametrization is derived
     4556      sai_grass_haarweg = 3.5
     4557
     4558      !                  100 - rh
     4559      !  rw = 2.0 * exp(----------)
     4560      !                    12
     4561
     4562      IF (SAI_present) then
     4563
     4564         ! External resistance according to Sutton & Fowler, 1993
     4565         rw = 2.0 * exp((100.0 - rh)/12.0)
     4566         rw = sai_grass_haarweg * rw
     4567
     4568         ! Frozen soil (from Depac v1):
     4569         IF (tsurf .lt. 0) rw = 200
     4570
     4571         ! Conductance:
     4572         gw = 1./rw
     4573      ELSE
     4574         ! no vegetation:
     4575         gw = 0.0
     4576      ENDIF
     4577
     4578    END SUBROUTINE rw_nh3_sutton
     4579
     4580
     4581
     4582!-------------------------------------------------------------------
     4583! rw_constant: compute constant external leaf conductance
     4584!-------------------------------------------------------------------
     4585    SUBROUTINE rw_constant(rw_val,SAI_present,gw)
     4586
     4587      ! Input/output variables:
     4588      REAL(kind=wp)   , INTENT(in)  :: rw_val      ! constant value of Rw
     4589      LOGICAL         , INTENT(in)  :: SAI_present
     4590      REAL(kind=wp)   , INTENT(out) :: gw          ! wernal leaf conductance (m/s)
     4591
     4592      ! Variables from module:
     4593      ! SAI_present: vegetation is present
     4594
     4595      ! Compute conductance:
     4596      IF (SAI_present .and. .not.missing(rw_val)) then
     4597         gw = 1./rw_val
     4598      ELSE
     4599         gw = 0.
     4600      ENDIF
     4601
     4602    END SUBROUTINE rw_constant
     4603
     4604
     4605
     4606!-------------------------------------------------------------------
     4607! rc_gstom: compute stomatal conductance
     4608!-------------------------------------------------------------------
     4609    SUBROUTINE rc_gstom( icmp, compnam, lu, LAI_present, lai, glrad, sinphi, t, rh, diffc, &
     4610         gstom, &
     4611         p)
     4612
     4613      ! input/output variables:
     4614      INTEGER(iwp),           INTENT(in) :: icmp           ! component index
     4615      CHARACTER(len=*),       INTENT(in) :: compnam        ! component name
     4616      INTEGER(iwp),           INTENT(in) :: lu             ! land use type , lu = 1,...,nlu
     4617      LOGICAL,                INTENT(in) :: LAI_present
     4618      REAL(kind=wp),          INTENT(in) :: lai            ! one-sided leaf area index
     4619      REAL(kind=wp),          INTENT(in) :: glrad          ! global radiation (W/m2)
     4620      REAL(kind=wp),          INTENT(in) :: sinphi         ! sin of solar elevation angle
     4621      REAL(kind=wp),          INTENT(in) :: t              ! temperature (C)
     4622      REAL(kind=wp),          INTENT(in) :: rh             ! relative humidity (%)
     4623      REAL(kind=wp),          INTENT(in) :: diffc              ! diffusion coefficient of the gas involved
     4624      REAL(kind=wp),          INTENT(out):: gstom              ! stomatal conductance (m/s)
     4625      REAL(kind=wp), OPTIONAL,INTENT(in) :: p              ! pressure (Pa)
     4626
     4627
     4628      ! Local variables
     4629      REAL(kind=wp)                      :: vpd            ! vapour pressure deficit (kPa)
     4630
     4631      REAL(kind=wp), PARAMETER           :: dO3 = 0.13e-4  ! diffusion coefficient of ozon (m2/s)
     4632
     4633
     4634      SELECT CASE(trim(compnam))
     4635
     4636         !CASE('HNO3','N2O5','NO3','H2O2') this routine is not CALLed for HNO3
     4637
     4638      CASE('NO','CO')
     4639         ! for NO stomatal uptake is neglected:
     4640         gstom = 0.0
     4641
     4642      CASE('NO2','O3','SO2','NH3')
     4643
     4644         ! IF vegetation present:
     4645         IF (LAI_present) then
     4646
     4647            IF (glrad .gt. 0.0) then
     4648               CALL rc_get_vpd(t,rh,vpd)
     4649               CALL rc_gstom_emb( lu, glrad, t, vpd, LAI_present, lai, sinphi, gstom, p )
     4650               gstom = gstom*diffc/dO3       ! Gstom of Emberson is derived for ozone
     4651            ELSE
     4652               gstom = 0.0
     4653            ENDIF
     4654         ELSE
     4655            ! no vegetation; zero conductance (infinite resistance):
     4656            gstom = 0.0
     4657         ENDIF
     4658
     4659      CASE default
     4660         message_string = 'Component '// trim(compnam) // ' not supported'
     4661         CALL message( 'rc_gstom', 'CM0459', 1, 2, 0, 6, 0 )
     4662      END SELECT
     4663
     4664    END SUBROUTINE rc_gstom
     4665
     4666
     4667
     4668!-------------------------------------------------------------------
     4669! rc_gstom_emb: stomatal conductance according to Emberson
     4670!-------------------------------------------------------------------
     4671    SUBROUTINE rc_gstom_emb( lu, glrad, T, vpd, LAI_present, lai, sinp, &
     4672         Gsto, p )
     4673      ! History
     4674      !   Original code from Lotos-Euros, TNO, M. Schaap
     4675      !   2009-08, M.C. van Zanten, Rivm
     4676      !     Updated and extended.
     4677      !   2009-09, Arjo Segers, TNO
     4678      !     Limitted temperature influence to range to avoid
     4679      !     floating point exceptions.
     4680      !
     4681      ! Method
     4682      !
     4683      !   Code based on Emberson et al, 2000, Env. Poll., 403-413
     4684      !   Notation conform Unified EMEP Model Description Part 1, ch 8
     4685      !
     4686      !   In the calculation of F_light the modification of L. Zhang 2001, AE to the PARshade and PARsun
     4687      !   parametrizations of Norman 1982 are applied
     4688      !
     4689      !   F_phen and F_SWP are set to 1
     4690      !
     4691      !   Land use types DEPAC versus Emberson (Table 5.1, EMEP model description)
     4692      !   DEPAC                     Emberson
     4693      !     1 = grass                 GR = grassland
     4694      !     2 = arable land           TC = temperate crops ( LAI according to RC = rootcrops)
     4695      !     3 = permanent crops       TC = temperate crops ( LAI according to RC = rootcrops)
     4696      !     4 = coniferous forest     CF = tempareate/boREAL(kind=wp) coniferous forest
     4697      !     5 = deciduous forest      DF = temperate/boREAL(kind=wp) deciduous forest
     4698      !     6 = water                 W  = water
     4699      !     7 = urban                 U  = urban
     4700      !     8 = other                 GR = grassland
     4701      !     9 = desert                DE = desert
     4702      !
     4703
     4704      ! Emberson specific declarations
     4705
     4706      ! Input/output variables:
     4707      INTEGER(iwp),               INTENT(in) :: lu                ! land use type, lu = 1,...,nlu
     4708      REAL(kind=wp),              INTENT(in) :: glrad             ! global radiation (W/m2)
     4709      REAL(kind=wp),              INTENT(in) :: T                 ! temperature (C)
     4710      REAL(kind=wp),              INTENT(in) :: vpd               ! vapour pressure deficit (kPa)
     4711      LOGICAL,                    INTENT(in) :: LAI_present
     4712      REAL(kind=wp),              INTENT(in) :: lai               ! one-sided leaf area index
     4713      REAL(kind=wp),              INTENT(in) :: sinp              ! sin of solar elevation angle
     4714      REAL(kind=wp),              INTENT(out):: Gsto              ! stomatal conductance (m/s)
     4715      REAL(kind=wp), OPTIONAL,    INTENT(in) :: p                 ! pressure (Pa)
     4716
     4717      ! Local variables:
     4718      REAL(kind=wp)                          ::  F_light, F_phen, F_temp, F_vpd, F_swp, bt, F_env
     4719      REAL(kind=wp)                          ::  pardir, pardiff, parshade, parsun, LAIsun, LAIshade, sinphi
     4720      REAL(kind=wp)                          ::  pres
     4721      REAL(kind=wp), PARAMETER               ::  p_sealevel = 1.01325e5    ! Pa
     4722
     4723
     4724      ! Check whether vegetation is present:
     4725      IF (LAI_present) then
     4726
     4727         ! calculation of correction factors for stomatal conductance
     4728         IF (sinp .le. 0.0) then 
     4729            sinphi = 0.0001
     4730         ELSE
     4731            sinphi = sinp
     4732         END IF
     4733
     4734         ! ratio between actual and sea-level pressure is used
     4735         ! to correct for height in the computation of par ;
     4736         ! should not exceed sea-level pressure therefore ...
     4737         IF ( present(p) ) then
     4738            pres = min( p, p_sealevel )
     4739         ELSE
     4740            pres = p_sealevel
     4741         ENDIF
     4742
     4743         ! Direct and diffuse par, Photoactive (=visible) radiation:
     4744         CALL par_dir_diff( glrad, sinphi, pres, p_sealevel, pardir, pardiff )
     4745
     4746         ! par for shaded leaves (canopy averaged):
     4747         parshade = pardiff*exp(-0.5*LAI**0.7)+0.07*pardir*(1.1-0.1*LAI)*exp(-sinphi) ! Norman, 1982
     4748         IF (glrad .gt. 200 .and. LAI .gt. 2.5) then
     4749            parshade = pardiff*exp(-0.5*LAI**0.8)+0.07*pardir*(1.1-0.1*LAI)*exp(-sinphi) ! Zhang et al., 2001
     4750         END IF
     4751
     4752         ! par for sunlit leaves (canopy averaged):
     4753         ! alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5
     4754         parsun = pardir*0.5/sinphi + parshade ! Norman, 1982
     4755         IF (glrad .gt. 200 .and. LAI .gt. 2.5) then
     4756            parsun = pardir**0.8*0.5/sinphi + parshade ! Zhang et al., 2001
     4757         END IF
     4758
     4759         ! leaf area index for sunlit and shaded leaves:
     4760         IF (sinphi .gt. 0) then
     4761            LAIsun = 2*sinphi*(1-exp(-0.5*LAI/sinphi ))
     4762            LAIshade = LAI -LAIsun
     4763         ELSE
     4764            LAIsun = 0
     4765            LAIshade = LAI
     4766         END IF
     4767
     4768         F_light = (LAIsun*(1 - exp(-1.*alpha(lu)*parsun)) + LAIshade*(1 - exp(-1.*alpha(lu)*parshade)))/LAI
     4769
     4770         F_light = max(F_light,F_min(lu))
     4771
     4772         !  temperature influence; only non-zero within range [Tmin,Tmax]:
     4773         IF ( (Tmin(lu) < T) .and. (T < Tmax(lu)) ) then
     4774            BT = (Tmax(lu)-Topt(lu))/(Topt(lu)-Tmin(lu))
     4775            F_temp = ((T-Tmin(lu))/(Topt(lu)-Tmin(lu))) * ((Tmax(lu)-T)/(Tmax(lu)-Topt(lu)))**BT
     4776         ELSE
     4777            F_temp = 0.0
     4778         END IF
     4779         F_temp = max(F_temp,F_min(lu))
     4780
     4781         ! vapour pressure deficit influence
     4782         F_vpd = min(1.,((1.-F_min(lu))*(vpd_min(lu)-vpd)/(vpd_min(lu)-vpd_max(lu)) + F_min(lu) ))
     4783         F_vpd = max(F_vpd,F_min(lu))
     4784
     4785         F_swp = 1.
     4786
     4787         ! influence of phenology on stom. conductance
     4788         ! ignored for now in DEPAC since influence of F_phen on lu classes in use is negligible.
     4789         ! When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore
     4790         F_phen = 1.
     4791
     4792         ! evaluate total stomatal conductance
     4793         F_env = F_temp*F_vpd*F_swp
     4794         F_env = max(F_env,F_min(lu))
     4795         gsto = G_max(lu) * F_light * F_phen * F_env
     4796
     4797         ! gstom expressed per m2 leafarea;
     4798         ! this is converted with LAI to m2 surface.
     4799         Gsto = lai*gsto    ! in m/s
     4800
     4801      ELSE
     4802         GSto = 0.0
     4803      ENDIF
     4804
     4805    END SUBROUTINE rc_gstom_emb
     4806
     4807
     4808
     4809!-------------------------------------------------------------------
     4810! par_dir_diff
     4811!-------------------------------------------------------------------
     4812    SUBROUTINE par_dir_diff(glrad,sinphi,pres,pres_0,par_dir,par_diff)
     4813      !
     4814      !     Weiss, A., Norman, J.M. (1985) Partitioning solar radiation into direct and
     4815      !     diffuse, visible and near-infrared components. Agric. Forest Meteorol.
     4816      !     34, 205-213.
     4817      !
     4818      !     From a SUBROUTINE obtained from Leiming Zhang (via Willem Asman),
     4819      !     MeteoroLOGICAL Service of Canada
     4820      !     e-mail: leiming.zhang@ec.gc.ca
     4821      !
     4822      !     Leiming uses solar irradiance. This should be equal to global radiation and
     4823      !     Willem Asman set it to global radiation
     4824
     4825      REAL(kind=wp),    INTENT(in)  :: glrad             ! global radiation (W m-2)
     4826      REAL(kind=wp),    INTENT(in)  :: sinphi            ! sine of the solar elevation
     4827      REAL(kind=wp),    INTENT(in)  :: pres              ! actual pressure (to correct for height) (Pa)
     4828      REAL(kind=wp),    INTENT(in)  :: pres_0            ! pressure at sea level (Pa)
     4829      REAL(kind=wp),    INTENT(out) :: par_dir           ! PAR direct : visible (photoactive) direct beam radiation (W m-2)
     4830      REAL(kind=wp),    INTENT(out) :: par_diff          ! PAR diffuse: visible (photoactive) diffuse radiation (W m-2)
     4831
     4832      !     fn              = near-infrared direct beam fraction (DIMENSIONless)
     4833      !     fv              = PAR direct beam fraction (DIMENSIONless)
     4834      !     ratio           = ratio measured to potential solar radiation (DIMENSIONless)
     4835      !     rdm             = potential direct beam near-infrared radiation (W m-2); "potential" means clear-sky
     4836      !     rdn             = potential diffuse near-infrared radiation (W m-2)
     4837      !     rdu             = visible (PAR) direct beam radiation (W m-2)
     4838      !     rdv             = potential visible (PAR) diffuse radiation (W m-2)
     4839      !     rn              = near-infrared radiation (W m-2)
     4840      !     rv              = visible radiation (W m-2)
     4841      !     ww              = water absorption in the near infrared for 10 mm of precipitable water
     4842
     4843      REAL(kind=wp) :: rdu,rdv,ww,rdm,rdn,rv,rn,ratio,sv,fv
     4844
     4845      !     Calculate visible (PAR) direct beam radiation
     4846      !     600 W m-2 represents average amount of PAR (400-700 nm wavelength)
     4847      !     at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2)
     4848      rdu=600.*exp(-0.185*(pres/pres_0)/sinphi)*sinphi
     4849
     4850      !     Calculate potential visible diffuse radiation
     4851      rdv=0.4*(600.- rdu)*sinphi
     4852
     4853      !     Calculate the water absorption in the-near infrared
     4854      ww=1320*10**( -1.195+0.4459*log10(1./sinphi)-0.0345*(log10(1./sinphi))**2 )
     4855
     4856      !     Calculate potential direct beam near-infrared radiation
     4857      rdm=(720.*exp(-0.06*(pres/pres_0)/sinphi)-ww)*sinphi     !720 = solar constant - 600
     4858
     4859      !     Calculate potential diffuse near-infrared radiation
     4860      rdn=0.6*(720-rdm-ww)*sinphi
     4861
     4862      ! Compute visible and near-infrared radiation
     4863      rv=max(0.1,rdu+rdv)
     4864      rn=max(0.01,rdm+rdn)
     4865
     4866      ! Compute ratio between input global radiation and total radiation computed here
     4867      ratio=min(0.9,glrad/(rv+rn))
     4868
     4869      !     Calculate total visible radiation
     4870      sv=ratio*rv
     4871
     4872      !     Calculate fraction of PAR in the direct beam
     4873      fv=min(0.99, (0.9-ratio)/0.7)            ! help variable
     4874      fv=max(0.01,rdu/rv*(1.0-fv**0.6667))     ! fraction of PAR in the direct beam
     4875
     4876      ! Compute direct and diffuse parts of PAR
     4877      par_dir=fv*sv
     4878      par_diff=sv-par_dir
     4879
     4880    END SUBROUTINE par_dir_diff
     4881
     4882!-------------------------------------------------------------------
     4883! rc_get_vpd: get vapour pressure deficit (kPa)
     4884!-------------------------------------------------------------------
     4885    SUBROUTINE rc_get_vpd(temp, relh,vpd)
     4886
     4887      IMPLICIT NONE
     4888
     4889      ! Input/output variables:
     4890      REAL(kind=wp)    , INTENT(in)  :: temp   ! temperature (C)
     4891      REAL(kind=wp)    , INTENT(in)  :: relh   ! relative humidity (%)
     4892      REAL(kind=wp)    , INTENT(out) :: vpd    ! vapour pressure deficit (kPa)
     4893
     4894      ! Local variables:
     4895      REAL(kind=wp) :: esat
     4896
     4897      ! fit PARAMETERs:
     4898      REAL(kind=wp), PARAMETER :: a1 = 6.113718e-1
     4899      REAL(kind=wp), PARAMETER :: a2 = 4.43839e-2
     4900      REAL(kind=wp), PARAMETER :: a3 = 1.39817e-3
     4901      REAL(kind=wp), PARAMETER :: a4 = 2.9295e-5
     4902      REAL(kind=wp), PARAMETER :: a5 = 2.16e-7
     4903      REAL(kind=wp), PARAMETER :: a6 = 3.0e-9
     4904
     4905      ! esat is saturation vapour pressure (kPa) at temp(C) following monteith(1973)
     4906      esat = a1 + a2*temp + a3*temp**2 + a4*temp**3 + a5*temp**4 + a6*temp**5
     4907      vpd  = esat * (1-relh/100)
     4908
     4909    END SUBROUTINE rc_get_vpd
     4910
     4911
     4912
     4913!-------------------------------------------------------------------
     4914! rc_gsoil_eff: compute effective soil conductance
     4915!-------------------------------------------------------------------
     4916    SUBROUTINE rc_gsoil_eff(icmp,lu,sai,ust,nwet,t,gsoil_eff)
     4917
     4918      ! Input/output variables:
     4919      INTEGER(iwp)    , INTENT(in)  :: icmp           ! component index
     4920      INTEGER(iwp)    , INTENT(in)  :: lu             ! land use type, lu = 1,...,nlu
     4921      REAL(kind=wp)   , INTENT(in)  :: sai            ! surface area index
     4922      REAL(kind=wp)   , INTENT(in)  :: ust            ! friction velocity (m/s)
     4923      INTEGER(iwp)    , INTENT(in)  :: nwet           ! index for wetness
     4924      ! nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     4925      ! N.B. this routine cannot be CALLed with nwet = 9,
     4926      ! nwet = 9 should be handled outside this routine.
     4927      REAL(kind=wp)   , INTENT(in)  :: t              ! temperature (C)
     4928      REAL(kind=wp)   , INTENT(out) :: gsoil_eff      ! effective soil conductance (m/s)
     4929
     4930      ! local variables:
     4931      REAL(kind=wp)                 :: rinc           ! in canopy resistance  (s/m)
     4932      REAL(kind=wp)                 :: rsoil_eff      ! effective soil resistance (s/m)
     4933
     4934
     4935
     4936      ! Soil resistance (numbers matched with lu_classes and component numbers)
     4937      REAL(kind=wp), PARAMETER     :: rsoil(nlu_dep,ncmp) = reshape( (/ &
     4938                                ! grs    ara    crp    cnf    dec    wat    urb    oth    des    ice    sav    trf    wai    med    sem
     4939           1000.,  200.,  200.,  200.,  200., 2000.,  400., 1000., 2000., 2000., 1000.,  200., 2000.,  200.,  400., & ! O3
     4940           1000., 1000., 1000., 1000., 1000.,   10., 1000., 1000., 1000.,  500., 1000., 1000.,   10., 1000., 1000., & ! SO2
     4941           1000., 1000., 1000., 1000., 1000., 2000., 1000., 1000., 1000., 2000., 1000., 1000., 2000., 1000., 1000., & ! NO2
     4942           -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., & ! NO
     4943           100.,  100.,  100.,  100.,  100.,   10.,  100.,  100.,  100., 1000.,  100.,  100.,   10.,  100.,  100.,  & ! NH3
     4944           -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., & ! CO
     4945           -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! NO3
     4946           -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! HNO3
     4947           -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! N2O5
     4948           -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999. /),& ! H2O2   
     4949           (/nlu_dep,ncmp/) )
     4950
     4951      !  o3       so2       no2        no       nh3        co    no3   hno3   n2o5   h2o2
     4952
     4953
     4954      REAL(kind=wp), PARAMETER     :: rsoil_wet(ncmp)    = (/2000.,      10.,    2000.,    -999.,      10.,    -999., -999., -999., -999., -999./)
     4955      REAL(kind=wp), PARAMETER     :: rsoil_frozen(ncmp) = (/2000.,     500.,    2000.,    -999.,    1000.,    -999., -999., -999., -999., -999./)
     4956
     4957
     4958
     4959      ! Compute in canopy (in crop) resistance:
     4960      CALL rc_rinc(lu,sai,ust,rinc)
     4961
     4962      ! Check for missing deposition path:
     4963      IF (missing(rinc)) then
     4964         rsoil_eff = -9999.
     4965      ELSE
     4966         ! Frozen soil (temperature below 0):
     4967         IF (t .lt. 0.0) then
     4968            IF (missing(rsoil_frozen(icmp))) then
     4969               rsoil_eff = -9999.
     4970            ELSE
     4971               rsoil_eff = rsoil_frozen(icmp) + rinc
     4972            ENDIF
     4973         ELSE
     4974            ! Non-frozen soil; dry:
     4975            IF (nwet .eq. 0) then
     4976               IF (missing(rsoil(lu,icmp))) then
     4977                  rsoil_eff = -9999.
     4978               ELSE
     4979                  rsoil_eff = rsoil(lu,icmp) + rinc
     4980               ENDIF
     4981
     4982               ! Non-frozen soil; wet:
     4983            ELSEIF (nwet .eq. 1) then
     4984               IF (missing(rsoil_wet(icmp))) then
     4985                  rsoil_eff = -9999.
     4986               ELSE
     4987                  rsoil_eff = rsoil_wet(icmp) + rinc
     4988               ENDIF
     4989            ELSE
     4990               message_string = 'nwet can only be 0 or 1'
     4991               CALL message( 'rc_gsoil_eff', 'CM0460', 1, 2, 0, 6, 0 )
     4992            ENDIF
     4993         ENDIF
     4994      ENDIF
     4995
     4996      ! Compute conductance:
     4997      IF (rsoil_eff .gt. 0.0) then
     4998         gsoil_eff = 1./rsoil_eff
     4999      ELSE
     5000         gsoil_eff = 0.0
     5001      ENDIF
     5002
     5003    END SUBROUTINE rc_gsoil_eff
     5004
     5005
     5006
     5007!-------------------------------------------------------------------
     5008! rc_rinc: compute in canopy (or in crop) resistance
     5009! van Pul and Jacobs, 1993, BLM
     5010!-------------------------------------------------------------------
     5011    SUBROUTINE rc_rinc(lu,sai,ust,rinc)
     5012
     5013
     5014      ! Input/output variables:
     5015      INTEGER(iwp)    , INTENT(in)  :: lu          ! land use class, lu = 1, ..., nlu
     5016      REAL(kind=wp)   , INTENT(in)  :: sai         ! surface area index
     5017      REAL(kind=wp)   , INTENT(in)  :: ust         ! friction velocity (m/s)
     5018      REAL(kind=wp)   , INTENT(out) :: rinc        ! in canopy resistance (s/m)
     5019
     5020
     5021      ! b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 IF not applicable)
     5022      ! h = vegetation height (m)                             grass arabl  crops conIF decid   water  urba   othr  desr   ice  sav  trf   wai   med semi
     5023      REAL(kind=wp), DIMENSION(nlu_dep), PARAMETER :: b = (/ -999,  14,    14,   14,   14,    -999,  -999, -999, -999, -999, -999, 14, -999, 14, 14 /)
     5024      REAL(kind=wp), DIMENSION(nlu_dep), PARAMETER :: h = (/ -999,  1,     1,    20,   20,    -999,  -999, -999, -999, -999, -999, 20,  -999, 1 ,  1 /)
     5025
     5026      ! Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0:
     5027      IF (b(lu) .gt. 0.0) then
     5028
     5029         ! Check for u* > 0 (otherwise denominator = 0):
     5030         IF (ust .gt. 0.0) then
     5031            rinc = b(lu)*h(lu)*sai/ust
     5032         ELSE
     5033            rinc = 1000.0
     5034         ENDIF
     5035      ELSE
     5036         IF (lu .eq. ilu_grass .or. lu .eq. ilu_other ) then
     5037            rinc = -999. ! no deposition path for grass, other, and semi-natural
     5038         ELSE
     5039            rinc = 0.0   ! no in-canopy resistance
     5040         ENDIF
     5041      ENDIF
     5042
     5043    END SUBROUTINE rc_rinc
     5044
     5045
     5046
     5047!-------------------------------------------------------------------
     5048! rc_rctot: compute total canopy (or surface) resistance Rc
     5049!-------------------------------------------------------------------
     5050    SUBROUTINE rc_rctot(gstom,gsoil_eff,gw,gc_tot,rc_tot)
     5051
     5052      ! Input/output variables:
     5053      REAL(kind=wp), INTENT(in)  :: gstom        ! stomatal conductance (s/m)
     5054      REAL(kind=wp), INTENT(in)  :: gsoil_eff    ! effective soil conductance (s/m)
     5055      REAL(kind=wp), INTENT(in)  :: gw           ! external leaf conductance (s/m)
     5056      REAL(kind=wp), INTENT(out) :: gc_tot       ! total canopy conductance (m/s)
     5057      REAL(kind=wp), INTENT(out) :: rc_tot       ! total canopy resistance Rc (s/m)
     5058
     5059      ! Total conductance:
     5060      gc_tot = gstom + gsoil_eff + gw
     5061
     5062      ! Total resistance (note: gw can be negative, but no total emission allowed here):
     5063      IF (gc_tot .le. 0.0 .or. gw .lt. 0.0) then
     5064         rc_tot = -9999.
     5065      ELSE
     5066         rc_tot = 1./gc_tot
     5067      ENDIF
     5068
     5069    END SUBROUTINE rc_rctot
     5070
     5071
     5072
     5073    !-------------------------------------------------------------------
     5074    ! rc_comp_point_rc_eff: calculate the effective resistance Rc
     5075    ! based on one or more compensation points
     5076    !-------------------------------------------------------------------
     5077    !
     5078    ! NH3rc (see depac v3.6 is based on Avero workshop Marc Sutton. p. 173.
     5079    ! Sutton 1998 AE 473-480)
     5080    !
     5081    ! Documentation by Ferd Sauter, 2008; see also documentation block in header of depac subroutine.
     5082    ! FS 2009-01-29: variable names made consistent with DEPAC
     5083    ! FS 2009-03-04: use total compensation point
     5084    !
     5085    ! C: with total compensation point   ! D: approximation of C
     5086    !                                    !    with classical approach
     5087    !  zr --------- Catm                 !  zr --------- Catm   
     5088    !         |                          !         |       
     5089    !         Ra                         !         Ra     
     5090    !         |                          !         |       
     5091    !         Rb                         !         Rb     
     5092    !         |                          !         |       
     5093    !  z0 --------- Cc                   !  z0 --------- Cc
     5094    !         |                          !         |             
     5095    !        Rc                          !        Rc_eff         
     5096    !         |                          !         |             
     5097    !     --------- Ccomp_tot            !     --------- C=0
     5098    !
     5099    !
     5100    ! The effective Rc is defined such that instead of using
     5101    !
     5102    !   F = -vd*[Catm - Ccomp_tot]                                    (1)
     5103    !
     5104    ! we can use the 'normal' flux formula
     5105    !
     5106    !   F = -vd'*Catm,                                                (2)
     5107    !
     5108    ! with vd' = 1/(Ra + Rb + Rc')                                    (3)
     5109    !
     5110    ! and Rc' the effective Rc (rc_eff).
     5111    !                                                (Catm - Ccomp_tot)
     5112    ! vd'*Catm = vd*(Catm - Ccomp_tot) <=> vd' = vd* ------------------
     5113    !                                                      Catm
     5114    !
     5115    !                                        (Catm - Ccomp_tot)
     5116    ! 1/(Ra + Rb + Rc') = (1/Ra + Rb + Rc) * ------------------
     5117    !                                              Catm
     5118    !
     5119    !                                          Catm
     5120    ! (Ra + Rb + Rc') = (Ra + Rb + Rc) * ------------------
     5121    !                                     (Catm - Ccomp_tot)
     5122    !
     5123    !                              Catm
     5124    ! Rc' = (Ra + Rb + Rc) * ------------------ - Ra - Rb
     5125    !                        (Catm - Ccomp_tot)
     5126    !
     5127    !                        Catm                           Catm
     5128    ! Rc' = (Ra + Rb) [------------------ - 1 ] + Rc * ------------------
     5129    !                  (Catm - Ccomp_tot)              (Catm - Ccomp_tot)
     5130    !
     5131    ! Rc' = [(Ra + Rb)*Ccomp_tot + Rc*Catm ] / (Catm - Ccomp_tot)
     5132    !
     5133    ! This is not the most efficient way to DO this;
     5134    ! in the current LE version, (1) is used directly in the CALLing routine
     5135    ! and this routine is not used anymore.
     5136    !
     5137    ! -------------------------------------------------------------------------------------------
     5138    ! In fact it is the question IF this correct; note the difference in differential equations
     5139    !
     5140    !   dCatm                              !   dCatm
     5141    ! H ----- = F = -vd'*Catm,             ! H ----- = F = -vd*[Catm - Ccomp_tot],
     5142    !    dt                                !    dt
     5143    !
     5144    ! where in the left colum vd' is a function of Catm and not a constant anymore.
     5145    !
     5146    ! Another problem is that IF Catm -> 0, we end up with an infinite value of the exchange
     5147    ! velocity vd.
     5148    ! -------------------------------------------------------------------------------------------
     5149
     5150    SUBROUTINE rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff)
     5151
     5152      IMPLICIT NONE
     5153
     5154      ! Input/output variables:
     5155      REAL(kind=wp), INTENT(in)  :: ccomp_tot  ! total compensation point (weighed average of separate compensation points) (ug/m3)
     5156      REAL(kind=wp), INTENT(in)  :: catm       ! atmospheric concentration (ug/m3)
     5157      REAL(kind=wp), INTENT(in)  :: ra         ! aerodynamic resistance (s/m)
     5158      REAL(kind=wp), INTENT(in)  :: rb         ! boundary layer resistance (s/m)
     5159      REAL(kind=wp), INTENT(in)  :: rc_tot     ! total canopy resistance (s/m)
     5160      REAL(kind=wp), INTENT(out) :: rc_eff     ! effective total canopy resistance (s/m)
     5161
     5162      ! Compute effective resistance:
     5163      IF ( ccomp_tot == 0.0 ) then
     5164         ! trace with no compensiation point ( or compensation point equal to zero)
     5165         rc_eff = rc_tot
     5166
     5167      ELSE IF ( ccomp_tot > 0 .and. ( abs( catm - ccomp_tot ) < 1.e-8 ) ) then
     5168         ! surface concentration (almoast) equal to atmospheric concentration
     5169         ! no exchange between surface and atmosphere, infinite RC --> vd=0
     5170         rc_eff = 9999999999.
     5171
     5172      ELSE IF ( ccomp_tot > 0 ) then
     5173         ! compensation point available, calculate effective restisance
     5174         rc_eff = ((ra + rb)*ccomp_tot + rc_tot*catm)/(catm-ccomp_tot)
     5175
     5176      ELSE
     5177         rc_eff = -999.
     5178         message_string = 'This should not be possible, check ccomp_tot'
     5179         CALL message( 'rc_comp_point_rc_eff', 'CM0461', 1, 2, 0, 6, 0 )
     5180      ENDIF
     5181
     5182      return
     5183    END SUBROUTINE rc_comp_point_rc_eff
     5184
     5185
     5186
     5187!-------------------------------------------------------------------
     5188! missing: check for data that correspond with a missing deposition path
     5189!          this data is represented by -999
     5190!-------------------------------------------------------------------
     5191
     5192    LOGICAL function missing(x)
     5193
     5194      REAL(kind=wp), INTENT(in) :: x
     5195
     5196      ! bandwidth for checking (in)equalities of floats
     5197      REAL(kind=wp), PARAMETER :: EPS = 1.0e-5
     5198
     5199      missing = (abs(x + 999.) .le. EPS)
     5200
     5201    END function missing
     5202
     5203
     5204
     5205    ELEMENTAL FUNCTION sedimentation_velocity(rhopart, partsize, slipcor, visc) RESULT (vs)
     5206
     5207
     5208      IMPLICIT NONE
     5209
     5210      ! --- in/out ---------------------------------
     5211
     5212      REAL(kind=wp), INTENT(in)        ::  rhopart    ! kg/m3
     5213      REAL(kind=wp), INTENT(in)        ::  partsize   ! m
     5214      REAL(kind=wp), INTENT(in)        ::  slipcor    ! m
     5215      REAL(kind=wp), INTENT(in)        ::  visc       ! viscosity
     5216      REAL(kind=wp)                    ::  vs
     5217
     5218      ! acceleration of gravity:
     5219      REAL(kind=wp), PARAMETER         ::  grav = 9.80665    ! m/s2
     5220
     5221      ! --- begin ----------------------------------
     5222
     5223      !viscosity & sedimentation velocity
     5224      vs = rhopart * (partsize**2) * grav * slipcor / (18*visc)
     5225
     5226    END FUNCTION sedimentation_velocity
     5227
     5228!------------------------------------------------------------------------
     5229! Boundary-layer deposition resistance following Zhang (2001)
     5230!------------------------------------------------------------------------
     5231    SUBROUTINE drydepo_aero_zhang_vd( vd, Rs, vs1, partsize, slipcor, &
     5232         nwet, tsurf, dens1, viscos1, &
     5233         luc, ftop_lu, ustar_lu)
     5234
     5235
     5236      IMPLICIT NONE 
     5237
     5238      ! --- in/out ---------------------------------
     5239
     5240      REAL(kind=wp), INTENT(out)        ::  vd          ! deposition velocity (m/s)
     5241      REAL(kind=wp), INTENT(out)        ::  Rs          ! Sedimentaion resistance (s/m)
     5242      REAL(kind=wp), INTENT(in)         ::  vs1         ! sedimentation velocity in lowest layer
     5243      REAL(kind=wp), INTENT(in)         ::  partsize    ! particle diameter (m)
     5244      REAL(kind=wp), INTENT(in)         ::  slipcor     ! slip correction factor
     5245      REAL(kind=wp), INTENT(in)         ::  tsurf       ! surface temperature (K)
     5246      REAL(kind=wp), INTENT(in)         ::  dens1       ! air density (kg/m3) in lowest layer
     5247      REAL(kind=wp), INTENT(in)         ::  viscos1     ! air viscosity in lowest layer
     5248      REAL(kind=wp), INTENT(in)         ::  ftop_lu     ! atmospheric resistnace Ra
     5249      REAL(kind=wp), INTENT(in)         ::  ustar_lu    ! friction velocity u*   
     5250
     5251      INTEGER(iwp),  INTENT(in)         ::  nwet        ! 1=rain, 9=snowcover
     5252      INTEGER(iwp),  INTENT(in)         ::  luc         ! DEPAC LU
     5253
     5254
     5255      ! --- const ----------------------------------
     5256
     5257
     5258
     5259      ! acceleration of gravity:
     5260      REAL(kind=wp), PARAMETER         ::  grav     = 9.80665    ! m/s2
     5261
     5262      REAL(kind=wp), PARAMETER         ::  beta     = 2.0
     5263      REAL(kind=wp), PARAMETER         ::  epsilon0 = 3.0
     5264      REAL(kind=wp), PARAMETER         ::  kb       = 1.38066e-23
     5265      REAL(kind=wp), PARAMETER         ::  pi = 3.141592654_wp  !< PI
     5266
     5267      REAL, PARAMETER :: alfa_lu(nlu_dep) = &
     5268           (/1.2,  1.2,   1.2,  1.0,  1.0,   100.0, 1.5,  1.2, 50.0, 100.0, 1.2, 1.0, 100.0, 1.2, 50.0/)   
     5269      !           grass arabl  crops conif decid  water  urba  othr  desr  ice   sav   trf  wai   med  sem
     5270      REAL, PARAMETER :: gamma_lu(nlu_dep) = &
     5271           (/0.54, 0.54,  0.54, 0.56, 0.56,  0.50,  0.56, 0.54, 0.58, 0.50, 0.54, 0.56, 0.50, 0.54, 0.54/)
     5272      !           grass arabl  crops conif decid  water  urba  othr  desr  ice   sav   trf   wai   med   sem   
     5273      REAL, PARAMETER ::A_lu(nlu_dep) = &   
     5274           (/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./)
     5275      !           grass arabl  crops conif decid  water  urba  othr  desr  ice   sav  trf   wai  med   sem     
     5276
     5277
     5278      ! --- local ----------------------------------
     5279
     5280      REAL(kind=wp)              ::  kinvisc, diff_part
     5281      REAL(kind=wp)              ::  schmidt,stokes, Ebrown, Eimpac, Einterc, Reffic
     5282
     5283      ! --- begin ----------------------------------
     5284
     5285      ! kinetic viscosity & diffusivity
     5286      kinvisc = viscos1 / dens1 ! only needed at surface
     5287
     5288      diff_part = kb * tsurf * slipcor / (3*pi*viscos1*partsize)
     5289
     5290      ! Schmidt number
     5291      schmidt = kinvisc / diff_part
     5292
     5293      ! calculate collection efficiencie E
     5294      Ebrown = Schmidt**(-gamma_lu(luc)) ! Brownian diffusion
     5295      ! determine Stokes number, interception efficiency
     5296      ! and sticking efficiency R (1 = no rebound)
     5297      IF ( luc == ilu_ice .or. nwet.eq.9 .or. luc == ilu_water_sea .or. luc == ilu_water_inland ) THEN
     5298         stokes=vs1*ustar_lu**2/(grav*kinvisc)
     5299         Einterc=0.0
     5300         Reffic=1.0
     5301      ELSE IF ( luc == ilu_other .or. luc == ilu_desert ) THEN !tundra of desert
     5302         stokes=vs1*ustar_lu**2/(grav*kinvisc)
     5303         Einterc=0.0
     5304         Reffic=exp(-Stokes**0.5)
     5305      ELSE
     5306         stokes=vs1*ustar_lu/(grav*A_lu(luc)*1.e-3)
     5307         Einterc=0.5*(partsize/(A_lu(luc)*1e-3))**2
     5308         Reffic=exp(-Stokes**0.5)
     5309      END IF
     5310      ! when surface is wet all particles DO not rebound:
     5311      IF(nwet.eq.1) Reffic = 1.0
     5312      ! determine impaction efficiency:
     5313      Eimpac = ( stokes / (alfa_lu(luc)+stokes) )**beta
     5314
     5315      ! sedimentation resistance:
     5316      Rs = 1.0 / ( epsilon0 * ustar_lu * (Ebrown+Eimpac+Einterc) * Reffic )
     5317
     5318      ! deposition velocity according to Seinfeld and Pandis (= SP, 2006; eq 19.7):
     5319      !
     5320      !            1
     5321      !    vd = ------------------ + vs
     5322      !         Ra + Rs + Ra*Rs*vs
     5323      !
     5324      !    where: Rs = Rb (in SP, 2006)
     5325      !
     5326      !
     5327      vd = 1.0 / ( ftop_lu + Rs + ftop_lu*Rs*vs1) + vs1
     5328
     5329
     5330
     5331    END SUBROUTINE drydepo_aero_zhang_vd
     5332
     5333
     5334
     5335    SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffc, Rb )   
     5336
     5337      ! Compute quasi-laminar boundary layer resistance as a function of landuse and tracer.
     5338      !
     5339      ! Original EMEP formulation by (Simpson et al, 2003) is used.
     5340      !
     5341
     5342
     5343      IMPLICIT NONE
     5344
     5345      ! --- in/out ---------------------------------
     5346
     5347      LOGICAL      , INTENT(in)      ::  is_water
     5348      REAL(kind=wp), INTENT(in)      ::  z0h           ! roughness length for heat
     5349      REAL(kind=wp), INTENT(in)      ::  ustar         ! friction velocity
     5350      REAL(kind=wp), INTENT(in)      ::  diffc         ! coefficient of diffusivity
     5351      REAL(kind=wp), INTENT(out)     ::  Rb            ! boundary layer resistance
     5352
     5353      ! --- const ----------------------------------
     5354
     5355      ! thermal diffusivity of dry air 20 C
     5356      REAL(kind=wp), PARAMETER       ::  thk = 0.19e-4         ! http://en.wikipedia.org/wiki/Thermal_diffusivity  (@ T=300K)
     5357      REAL(kind=wp), PARAMETER       ::  kappa_stab = 0.35     ! von Karman constant
     5358
     5359      !  --- begin ---------------------------------
     5360
     5361
     5362      ! Use Simpson et al. (2003)
     5363      !IF ( is_water ) THEN
     5364      !Rb = 1.0_wp / (kappa_stab*max(0.01_wp,ustar)) * log(z0h/diffc*kappa_stab*max(0.01_wp,ustar))
     5365      !TODO: Check Rb over water calculation!!!
     5366      !   Rb = 1.0_wp / (kappa_stab*max(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*max(0.1_wp,ustar))
     5367      !ELSE
     5368      Rb = 5.0_wp / max(0.01_wp,ustar) * (thk/diffc)**0.67_wp
     5369      !END IF
     5370
     5371    END SUBROUTINE get_rb_cell
     5372
     5373
     5374!-----------------------------------------------------------------
     5375!  Compute water vapor partial pressure (e_w)
     5376!  given specific humidity Q [(kg water)/(kg air)].
     5377!
     5378!  Use that gas law for volume V with temperature T
     5379!  holds for the total mixture as well as the water part:
     5380!
     5381!    R T / V = p_air / n_air = p_water / n_water
     5382!
     5383!  thus:
     5384!
     5385!    p_water = p_air n_water / n_air
     5386!
     5387!  Use:
     5388!    n_air =   m_air   /        xm_air
     5389!            [kg air]  /  [(kg air)/(mole air)]
     5390!  and:
     5391!    n_water =  m_air * Q  /     xm_water
     5392!              [kg water]  /  [(kg water)/(mole water)]
     5393!  thus:
     5394!    p_water = p_air Q / (xm_water/xm_air)
     5395!
     5396
     5397    ELEMENTAL FUNCTION watervaporpartialpressure( Q, p ) result ( p_w )
     5398
     5399
     5400      ! --- in/out ---------------------------------
     5401
     5402      REAL(kind=wp), INTENT(in)      ::  Q   ! specific humidity [(kg water)/(kg air)]
     5403      REAL(kind=wp), INTENT(in)      ::  p   ! air pressure [Pa]
     5404      REAL(kind=wp)                  ::  p_w  ! water vapor partial pressure [Pa]
     5405
     5406      ! --- const ----------------------------------
     5407
     5408      ! mole mass ratio:
     5409      REAL(kind=wp), PARAMETER  ::  eps = xm_H2O / xm_air  !  ~ 0.622
     5410
     5411      ! --- begin ----------------------------------
     5412
     5413      ! partial pressure of water vapor:
     5414      p_w = p * Q / eps
     5415
     5416    END function watervaporpartialpressure
     5417
     5418
     5419
     5420!------------------------------------------------------------------   
     5421!  Saturation vapor pressure.
     5422!  From (Stull 1988, eq. 7.5.2d):
     5423!
     5424!      e_sat = p0 exp( 17.67 * (T-273.16) / (T-29.66) )     [Pa]
     5425!
     5426!  where:
     5427!      p0 = 611.2 [Pa]   : reference pressure
     5428!
     5429!  Arguments:
     5430!      T  [K]  : air temperature
     5431!  Result:
     5432!      e_sat_w  [Pa]  : saturation vapor pressure
     5433!
     5434!  References:
     5435!      Roland B. Stull, 1988
     5436!      An introduction to boundary layer meteorology.
     5437!
     5438
     5439    ELEMENTAL FUNCTION saturationvaporpressure( T ) result( e_sat_w )
     5440
     5441
     5442      ! --- in/out ---------------------------------
     5443
     5444      REAL(kind=wp), INTENT(in)      ::  T        ! temperature [K]
     5445      REAL(kind=wp)                  ::  e_sat_w  ! saturation vapor pressure  [Pa]
     5446
     5447      ! --- const ----------------------------------
     5448
     5449      ! base pressure:
     5450      REAL(kind=wp), PARAMETER  ::  p0 = 611.2   ! [Pa]
     5451
     5452      ! --- begin ----------------------------------
     5453
     5454      ! saturation vapor pressure:
     5455      e_sat_w = p0 * exp( 17.67 * (T-273.16) / (T-29.66) )   ! [Pa]
     5456
     5457    END FUNCTION saturationvaporpressure
     5458
     5459
     5460
     5461!------------------------------------------------------------------------
     5462!  Relative humidity RH [%] is by definition:
     5463!
     5464!           e_w            # water vapor partial pressure
     5465!    Rh = -------- * 100
     5466!         e_sat_w          # saturation vapor pressure
     5467!
     5468!  Compute here from:
     5469!    Q specific humidity [(kg water)/(kg air)]
     5470!    T temperature [K]
     5471!    P pressure [Pa]
     5472!
     5473   
     5474    ELEMENTAL FUNCTION relativehumidity_from_specifichumidity( Q, T, p ) result( Rh )
     5475
     5476
     5477      ! --- in/out ---------------------------------
     5478
     5479      REAL(kind=wp), INTENT(in)      ::  Q   ! specific humidity [(kg water)/(kg air)]
     5480      REAL(kind=wp), INTENT(in)      ::  T   ! temperature [K]
     5481      REAL(kind=wp), INTENT(in)      ::  p   ! air pressure [Pa]
     5482      REAL(kind=wp)                  ::  Rh  ! relative humidity [%]
     5483
     5484      ! --- begin ----------------------------------
     5485
     5486      ! relative humidity:
     5487      Rh = watervaporpartialpressure( Q, p ) / saturationvaporpressure( T ) * 100.0
     5488
     5489    END FUNCTION relativehumidity_from_specifichumidity
     5490
     5491
     5492     
    24345493 END MODULE chemistry_model_mod
    24355494

  • palm/trunk/SOURCE/date_and_time_mod.f90

    r3298 r3458  
    1919!
    2020! Current revisions:
    21 ! -----------------
     21! ------------------
    2222!
    2323!
    … …  
    2525! -----------------
    2626! $Id$
     27! from chemistry branch r3443, banzhafs:
     28! Added initial hour_of_day, hour_of_year, day_of_year and month_of_year to
     29! init_date_and_time
     30!
     31! 3298 2018-10-02 12:21:11Z kanani
    2732! - Minor formatting (kanani)
    2833! - Added Routines for DEFAULT mode of chemistry emissions (Russo)
    … …  
    6166!>       already implemented changes for calculating it from date_init in
    6267!>       calc_date_and_time
     68!> @todo time_utc during spin-up 
    6369!------------------------------------------------------------------------------!
    6470 MODULE date_and_time_mod
    … …  
    9197    INTEGER(iwp)        ::  index_mm                                !< index months of the default emission mode
    9298    INTEGER(iwp)        ::  index_dd                                !< index days of the default emission mode
    93     INTEGER(iwp)        ::  index_hh                                !< index hours of the default emission mode
     99    INTEGER(iwp)        ::  index_hh                                !< index hours of the emission mode
     100    INTEGER(iwp)        ::  hours_since_reference_point             !< hours of current simulation
    94101
    95102    REAL(wp)            ::  time_utc                     !< current model time in UTC
    … …  
    102109    REAL(wp), PARAMETER ::  d_seconds_year = 1.0_wp / 31536000.0_wp !< inverse of the seconds per year (1/(365*86400))
    103110   
    104     CHARACTER(len=8)    ::  date_init = "31122018"                  !< Starting date of simulation: We selected this because it was a monday
     111    CHARACTER(len=8)    ::  date_init = "21062017"                  !< Starting date of simulation: We selected this because it was a monday
    105112 
    106113    !> --- Parameters
    … …  
    121128    END INTERFACE time_default_indices
    122129
     130    !-- Get hour index in the PRE-PROCESSED case of chemistry emissions :
     131    INTERFACE time_preprocessed_indices
     132       MODULE PROCEDURE time_preprocessed_indices
     133    END INTERFACE time_preprocessed_indices
     134
     135
    123136    !-- Calculate current date and time
    124137    INTERFACE calc_date_and_time
    … …  
    128141
    129142    !-- Public Interfaces
    130     PUBLIC calc_date_and_time, time_default_indices, init_date_and_time
     143    PUBLIC calc_date_and_time, time_default_indices, init_date_and_time, time_preprocessed_indices
    131144
    132145    !-- Public Variables
    … …  
    147160       IMPLICIT NONE
    148161
     162       !--    Variables Definition
     163       INTEGER ::  i_mon       !< Index for going through the different months
     164
    149165       IF  (day_of_year_init == 0) THEN
    150166          ! Day of the month at starting time
    … …  
    159175       ENDIF
    160176
    161     END SUBROUTINE init_date_and_time
    162 
    163 !------------------------------------------------------------------------------!
    164 ! Description:
    165 ! ------------
    166 !> Calculate current date and time of the simulation
    167 !------------------------------------------------------------------------------!
    168  
    169     SUBROUTINE calc_date_and_time
    170 
    171        IMPLICIT NONE
    172 
    173 !--    Variables Definition
    174        INTEGER                          :: i_mon       !< Index for going through the different months
    175 
    176        !> Update simulation time in seconds
    177        time_update = simulated_time-coupling_start_time
    178 
    179 !--    Calculate current day of the simulated time
    180        days_since_reference_point=INT(FLOOR( (time_utc_init + time_update) &
    181                                / 86400.0_wp ) )
    182 
    183 !--    Calculate actual UTC time                       
    184        time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
    185        
    186 !sB    PRILIMINARY workaround for time_utc bug concerning time_since_reference_point:
    187        time_utc_emis = MOD((time_utc_init + time_update), 86400.0_wp)     
    188 
    189 !--    Calculate initial day of the year: it is calculated only once. In fact, day_of_year_init is initialized to 0 and then a positive value is passed. This condition is also called only when day_of_year_init is not given in the namelist.
    190 
     177
     178       !-- Calculate initial hour of the day: the first hour of the day is from 00:00:00 to 00:59:59.
     179
     180       hour_of_day = INT( FLOOR( time_utc_init/3600.0_wp ) ) + 1
     181
     182       !-- Calculate initial day day_of_year_init in case date_init is given or day_of_year_init is given
    191183       IF ( day_of_year_init == 0 ) THEN
    192184
    … …  
    214206       ENDIF
    215207
     208       
     209       !-- Initial day of the year
     210       day_of_year = day_of_year_init
     211
     212       !-- Initial hour of the year
     213       hour_of_year = ( (day_of_year-1) * 24 ) + hour_of_day
     214
     215       !--Initial day of the month and month of the year
     216       !> --------------------------------------------------------------------------------
     217       !> The first case is when date_init is not provided: we only know day_of_year_init     
     218       IF ( month_of_year == 0 .AND. day_of_month == 0) THEN
     219
     220         
     221          IF ( day_of_year .LE. 31 ) THEN
     222
     223             month_of_year=1
     224             day_of_month=day_of_year
     225
     226          ELSE
     227
     228             DO i_mon=2,12   !january is considered in the first case
     229                IF ( day_of_year .LE. SUM(days(1:i_mon)) .AND. day_of_year .GT. SUM(days(1:(i_mon-1))) ) THEN
     230           
     231                   month_of_year=i_mon
     232
     233                   day_of_month=INT(MOD(day_of_year, SUM(days(1:(i_mon-1)))))
     234
     235                   GOTO 38
     236
     237                ENDIF
     238
     239             38 ENDDO
     240          ENDIF
     241       !> --------------------------------------------------------------------------------
     242       !> in the second condition both day of month and month_of_year are either given in input (passed to date_init) or we are in some day successive to the initial one, so that day_of_month has already be computed in previous step
     243       !>TBD: something to calculate the current year is missing
     244       ELSEIF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN
     245 
     246          !> calculate month_of_year. TBD: test the condition when day_of_year==31
     247 
     248          IF (day_of_year==1) THEN  !> this allows to turn from december to January when passing from a year to another
     249 
     250             month_of_year = 1
     251       
     252          ELSE IF (day_of_year .GT. 1 .AND. day_of_year .GT. SUM(days(1:month_of_year))) THEN
     253
     254             month_of_year = month_of_year + 1
     255
     256          ENDIF
     257
     258          !> calculate day_of_month
     259          IF ( month_of_year == 1 ) THEN
     260           
     261            day_of_month=day_of_year
     262
     263          ELSE
     264
     265            day_of_month=INT(MOD(day_of_year, SUM(days(1:(month_of_year-1)))))
     266
     267          ENDIF
     268
     269
     270       ELSE
     271
     272          !> Condition when date_init is provided but it is given in the wrong format
     273          message_string = 'date_init not provided in the namelist or'            //          &
     274                              ' given in the wrong format: MUST BE DDMMYYYY'                 
     275          CALL message( 'init_date_and_time', 'DT0102', 2, 2, 0, 6, 0 )
     276
     277       ENDIF
     278
     279
     280    END SUBROUTINE init_date_and_time
     281
     282!------------------------------------------------------------------------------!
     283! Description:
     284! ------------
     285!> Calculate current date and time of the simulation
     286!------------------------------------------------------------------------------!
     287 
     288    SUBROUTINE calc_date_and_time
     289
     290       IMPLICIT NONE
     291
     292!--    Variables Definition
     293       INTEGER                          :: i_mon       !< Index for going through the different months
     294
     295       !> Update simulation time in seconds
     296       time_update = simulated_time-coupling_start_time
     297
     298!--    Calculate current day of the simulated time
     299       days_since_reference_point=INT(FLOOR( (time_utc_init + time_update) &
     300                               / 86400.0_wp ) )
     301
     302!--    Calculate actual UTC time                       
     303       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
     304       
     305!sB    PRILIMINARY workaround for time_utc changes due to changes in time_since_reference_point in
     306!sB    radiation_model_mod during runtime:
     307       time_utc_emis = MOD((time_utc_init + time_update), 86400.0_wp)     
     308
     309!--    Calculate initial day of the year: it is calculated only once. In fact, day_of_year_init is initialized to 0 and then a positive value is passed. This condition is also called only when day_of_year_init is not given in the namelist.
     310
     311       IF ( day_of_year_init == 0 ) THEN
     312
     313          !> Condition for printing an error when date_init is not provided when day_of_year_init is not given in the namelist or when the format of the date is not the one required by PALM.
     314          IF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN
     315       
     316             IF ( month_of_year == 1 ) THEN  !!month of year is read in input
     317
     318                day_of_year_init = day_of_month
     319
     320             ELSE
     321
     322                day_of_year_init= SUM(days( 1:(month_of_year-1) )) + day_of_month  !day_of_month is read in input in this case
     323
     324             ENDIF
     325!kanani: Revise, we cannot force users to provide date_init, maybe set a default value?
     326!           ELSE
     327!
     328!              message_string = 'date_init not provided in the namelist or'             //          &
     329!                               ' given in the wrong format: MUST BE DDMMYYYY'                       
     330!              CALL message( 'calc_date_and_time', 'DT0100', 2, 2, 0, 6, 0 )
     331     
     332          ENDIF
     333
     334       ENDIF
     335
    216336      !-- Calculate actual hour of the day: the first hour of the day is from 00:00:00 to 00:59:59.
    217337
    … …  
    302422
    303423      ENDIF       
    304      
     424
    305425    END SUBROUTINE calc_date_and_time
    306426
    307 !TBD: check the routines used for update of emission indices in the DEFAULT MODE
     427
     428!------------------------------------------------------------------------------!
     429! Description:
     430! ------------
     431!> This routine determines the time factor index in the PRE-PROCESSED emissions mode.
     432!------------------------------------------------------------------------------!
     433
     434 SUBROUTINE time_preprocessed_indices(index_hh)
     435
     436    USE indices
     437
     438    IMPLICIT NONE
     439
     440!
     441!-- In/output
     442    INTEGER, INTENT(INOUT) ::  index_hh    !> Index Hour
     443!
     444!-- Additional Variables for calculateing indices
     445!-- Constants
     446    INTEGER, PARAMETER ::  nhour = 24
     447
     448    IF (days_since_reference_point == 0) THEN
     449
     450       index_hh=hour_of_day
     451
     452    ELSE
     453
     454       index_hh=(days_since_reference_point*nhour)+(hour_of_day)
     455
     456    ENDIF
     457
     458
     459 END SUBROUTINE time_preprocessed_indices
     460
     461
    308462!------------------------------------------------------------------------------!
    309463! Description:

  • palm/trunk/SOURCE/flow_statistics.f90

    r3452 r3458  
    2525! -----------------
    2626! $Id$
     27! from chemistry branch r3443, basit:
     28! bug fixed in chemistry profiles
     29!
     30! 3452 2018-10-30 13:13:34Z schwenkel
    2731! Bugfix for profiles output
    2832!
    … …  
    18431847       ENDIF
    18441848
    1845        IF ( air_chemistry .AND. max_pr_cs > 0 )  THEN
     1849       IF ( air_chemistry  .AND. max_pr_cs > 0 )  THEN
    18461850          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1847              CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+1,0), sums(nzb,pr_palm+1), &
    1848                                  nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
     1851             DO  i = 1, max_pr_cs
     1852                CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+i,0),       &
     1853                                    sums(nzb,pr_palm+max_pr_user+i),           &
     1854                                    nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
     1855             ENDDO
    18491856       ENDIF
    18501857

  • palm/trunk/SOURCE/init_3d_model.f90

    r3448 r3458  
    2525! -----------------
    2626! $Id$
     27! from chemistry branch r3443, basit:
     28! bug fixed in sums and sums_l for chemistry profile output
     29!
     30! 3448 2018-10-29 18:14:31Z kanani
    2731! Add biometeorology
    2832!
    … …  
    711715              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                &
    712716              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                  &
    713               sums(nzb:nzt+1,pr_palm+max_pr_user),                             &
    714               sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),      &
     717              sums(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs),                   &
     718              sums_l(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs,0:threads_per_task-1),      &
    715719              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),    &
    716720              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                   &

  • palm/trunk/SOURCE/netcdf_data_input_mod.f90

    r3429 r3458  
    2525! -----------------
    2626! $Id$
     27! from chemistry branch r3443, banzhafs, Russo:
     28! Uncommented lines on dimension of surface_fractions
     29! Removed par_emis_time_factor variable, moved to chem_emissions_mod
     30! Initialized nspec and other emission variables at time of declaration
     31! Modified EXPERT mode to PRE-PROCESSED mode
     32! Introduced Chemistry static netcdf file
     33! Added the routine for reading-in netcdf data for chemistry
     34! Added routines to get_variable interface specific for chemistry files
     35!
     36! 3429 2018-10-25 13:04:23Z knoop
    2737! add default values of origin_x/y/z
    2838!
    … …  
    368378
    369379       !-DIMENSIONS
    370        INTEGER(iwp)                                 :: nspec                     !< number of chem species for which emission values are provided
    371        INTEGER(iwp)                                 :: ncat                      !< number of emission categories
    372        INTEGER(iwp)                                 :: nvoc                      !< number of VOCs components
    373        INTEGER(iwp)                                 :: npm                       !< number of PMs components
     380       INTEGER(iwp)                                 :: nspec=0                   !< number of chem species for which emission values are provided
     381       INTEGER(iwp)                                 :: ncat=0                    !< number of emission categories
     382       INTEGER(iwp)                                 :: nvoc=0                    !< number of VOCs components
     383       INTEGER(iwp)                                 :: npm=0                     !< number of PMs components
    374384       INTEGER(iwp)                                 :: nnox=2                    !< number of NOx components: NO and NO2
    375385       INTEGER(iwp)                                 :: nsox=2                    !< number of SOx components: SO and SO4
    … …  
    391401       INTEGER(iwp),ALLOCATABLE, DIMENSION(:)       :: species_index             !< Index of emission chem species
    392402
    393        REAL(wp), DIMENSION(24)                      :: par_emis_time_factor      !< time factors
    394                                                                                  !  for the parameterized mode: these are fixed for each hour
    395                                                                                  !  of a single day.
    396403       REAL(wp),ALLOCATABLE, DIMENSION(:)           :: xm                        !< Molecular masses of emission chem species
    397404
    … …  
    928935         ENDDO
    929936
    930           !Assign Constant Values of time factors:
    931           emt_att%par_emis_time_factor( : ) = (/ 0.01, 0.01, 0.01, 0.02, 0.03, 0.07, 0.09, 0.09, 0.05, 0.03, 0.03, 0.03, &
    932                                                  0.03, 0.03, 0.03, 0.04, 0.05, 0.09, 0.09, 0.09, 0.04, 0.02, 0.01, 0.01 /)
    933937
    934938       !> DEFAULT AND PRE-PROCESSED MODE

  • palm/trunk/SOURCE/palm.f90

    r3337 r3458  
    2525! -----------------
    2626! $Id$
     27! from chemistry branch r3443, forkel:
     28! removed double do_emis check around CALL chem_init
     29! replaced call to calc_date_and_time to init_date_and_time
     30!
     31! 3337 2018-10-12 15:17:09Z kanani
    2732! (from branch resler)
    2833! Fix chemistry call
    … …  
    436441
    437442       CALL chem_init
     443
    438444!       CALL photolysis_init   ! probably also required for restart
    439 
    440        CALL init_date_and_time     !initialize the time of chemistry emissions
    441445
    442446    ENDIF
    … …  
    448452!
    449453!-- Initialize all necessary variables
    450     CALL calc_date_and_time !this is required for chemistry emissions
     454!
     455!-- Initial time for chem_emissions_mod
     456    CALL init_date_and_time
    451457
    452458    CALL init_3d_model

  • palm/trunk/SOURCE/prognostic_equations.f90

    r3386 r3458  
    2525! -----------------
    2626! $Id$
     27! remove duplicate USE chem_modules
     28! from chemistry branch r3443, banzhafs, basit:
     29! chem_depo call introduced
     30! code added for decycling chemistry
     31!
     32! 3386 2018-10-19 16:28:22Z gronemeier
    2733! Renamed tcm_prognostic to tcm_prognostic_equations
    2834!
    … …  
    340346    USE buoyancy_mod,                                                          &
    341347        ONLY:  buoyancy
    342 
    343     USE chemistry_model_mod,                                                   &
    344         ONLY:  chem_integrate, chem_species,  chem_prognostic_equations,       &
    345                nspec, nvar, spc_names, chem_boundary_conds
    346348           
    347349    USE chem_modules,                                                          &
    348         ONLY:  call_chem_at_all_substeps, chem_gasphase_on
     350        ONLY:  call_chem_at_all_substeps, chem_gasphase_on, cs_name, do_depo
    349351
    350352    USE chem_photolysis_mod,                                                   &
    351353        ONLY:  photolysis_control
    352354
    353     USE chem_modules,                                                          &
    354         ONLY:  call_chem_at_all_substeps, chem_gasphase_on, cs_name
     355    USE chemistry_model_mod,                                                   &
     356        ONLY:  chem_boundary_conds, chem_depo, chem_integrate,                 &
     357               chem_prognostic_equations, chem_species,                        &
     358               nspec, nvar, spc_names
    355359
    356360    USE control_parameters,                                                    &
    … …  
    495499
    496500                IF ( intermediate_timestep_count == 1 .OR.                        &
    497                      call_chem_at_all_substeps ) THEN
    498                    CALL chem_integrate (i,j)                                               
     501                     call_chem_at_all_substeps )  THEN
     502                   CALL chem_integrate (i,j)
     503                   IF ( do_depo )  THEN
     504                      CALL chem_depo(i,j)
     505                   ENDIF
    499506                ENDIF
    500507             ENDDO

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/Makefile

    r2696 r3458  
    2828
    2929distclean: clean
    30         rm -f $(bindir)/kp4.exe
     30        rm -f $(bindir)/kpp4palm.exe
    3131
    3232# ----------------------------------------------

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/bin/kpp4palm.ksh

    r3298 r3458  
    4040# -----------------
    4141# $Id$
     42# forkel   30. Oktober 2018: Integrating contents of kp4_compress into chem_gasphase_mod.f90
    4243# forkel   25. September 2018: Added cat for $MECH to pass mechanism name to kpp4palm
    4344# ketelsen 18. September 2018: Added cat for '#INLINE F90_GLOBAL'
    … …  
    121122      k)   KEEP="YES";;              # keep Working directory
    122123
    123       o)   OUTDIR=$OPTARG;;          # Output directory of Generated Code
    124 
    125       p)   PREFIX=$OPTARG;;          # Name Prefix
    126 
    127       s)   KPP_SOLVER=$OPTARG;;      # Update sample f90 code in the def_MECH directory
    128 
    129       u)   UPDT="YES";;              # keep Working directory
     124      o)   OUTDIR=$OPTARG;;          # Output directory of Generated Code =
     125                                     # SOURCE directory (do not change)
     126
     127      p)   PREFIX=$OPTARG;;          # Name Prefix (chem_gasphase_mod, do not change)
     128
     129      s)   KPP_SOLVER=$OPTARG;;      # Chosen solver (only Rosebrock solvers work for vector mode)
     130
     131      u)   UPDT="YES";;              # update mechanisms/def_$MECH/chem_gasphase_mod.f90
    130132
    131133      v)   MODE="vector";;           # Set to vector Mode
    … …  
    163165
    164166KPP_FILE_LIST="Initialize Integrator LinearAlgebra Jacobian Function Rates Util"
    165 
     167if [[ $MODE = "vector" ]]
     168 then
     169KPP_FILE_LIST="$KPP_FILE_LIST kp4_compress_subroutines"
     170fi
    166171
    167172KPP_SUBROUTINE_LIST="Initialize"
    … …  
    189194   # get vector Solver
    190195   cp $BASE/templates/${KPP_SOLVER}_vec.f90 ${KPP_SOLVER}.f90
     196   cp $BASE/templates/kp4_compress_header ${PREFIX}_kp4_compress_header.f90
     197   cp $BASE/templates/kp4_compress_subroutines ${PREFIX}_kp4_compress_subroutines.f90
    191198fi
    192199
    … …  
    207214         KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST Rosenbrock  FunTemplate JacTemplate"
    208215        KPP_INTERFACE_IGNORE="WAXPY"
     216
    209217      else
    210218        KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST FunTemplate JacTemplate"
     219        KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST kco_initialize kco_compress kco_finalize"
    211220      fi;;
    212221
    … …  
    246255
    247256KPP_INCLUDE_LIST="Parameters Global JacobianSP Monitor"
     257if [[ $MODE = "vector" ]]
     258 then
     259KPP_INCLUDE_LIST="$KPP_INCLUDE_LIST kp4_compress_header"
     260fi
    248261
    249262#Get definition Files
    … …  
    309322done
    310323
    311 echo start kp4.exe with arguments
     324echo start kpp4palm.exe with arguments
    312325echo $PREFIX $MODE $VLEN $DE_INDEX $DE_INDEX_FAST
    313326

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.C

    r3298 r3458  
    1212//Current revisions:
    1313//------------------
    14 / 
    15 / 
     14//
     15//
    1616//Former revisions:
    1717//-----------------
    18 //$Id: create_kpp_module.C 2470 2017-09-14 13:56:42Z forkel $
    19 // forkel 25.09.2018:  added  automatic line with meechanism name (read mech_list)
    20 // forkel 20.09.2018:  added  vl_glo = size(tempi,1)
    21 // ketelsen 18.09.2018: Removed creation of fill_ Subroutine and creation of calls thereof
    22 //
    23 // forkel 12.09.2018: Fix in order not to loose the values of qvap and fakt
    24 //
    25 // forkel 03.09.2018: Bug fix: moved kppi.add_line("    CALL initialize after fakt = fakti(is)
    26 // forkel June 2018: added qvap and fakt, re-established original uppercase/lowercase
    27 //                    Deleted definition of qvap,fakt in create_kpp_integrate again
    28 //
    29 // ketelsen July 2018: Changes for vector mode (edit_WAXPY, edit_FunTemplate, edit_JacTemplate,
    30 //                            some cleanup of comments, various changes in create_kpp_integrate)
    31 //
    32 // forkel June 2018: Added qvap and fakt, re-established original uppercase/lowercase
    33 //
    34 //
    35 / 2017-09-14 13:56:42Z forkel $
    36 // Removed preprocessor directive __chem again
    37 //
    38 // 2017-09-14 13:56:42Z forkel $
    39 //
     18//$Id: create_kpp_module.C 3453 2018-10-30 13:21:51Z forkel $
     19// exclude kco_compress from handling by global_variables2vector (30.10.2018, forkel)
     20//
     21// Added  automatic line with mechanism name (read mech_list) (25.09.2018, forkel)
     22//
     23// Added  vl_glo = size(tempi,1) (20.09.2018, forkel)
     24//
     25// Removed creation of fill_ Subroutine and creation of calls thereof (18.09.2018, ketelsen)
     26//
     27// Fix in order not to loose the values of qvap and fakt (12.09.2018, forkel)
     28//
     29// Bug fixes: moved kppi.add_line("    CALL initialize after fakt = fakti(is)
     30// Deleted definition of qvap,fakt in create_kpp_integrate again (03.09.2018, forkel)
     31//
     32// Changes for vector mode (edit_WAXPY, edit_FunTemplate, edit_JacTemplate,
     33// some cleanup of comments, various changes in create_kpp_integrate) (July 2018, ketelsen)
     34//
     35//
     36// Added qvap and fakt                            (June 2018, forkel)
     37// --> Change in module_header: qvap, fakt added  (June 2018, forkel)
     38//
     39// re-established original uppercase/lowercase     (June 2018, forkel)
     40// --> Change in module_header: reset case in  Initialize, Integrate, and
     41//     Update_rconst                               (June 2018, forkel)
     42//
     43// Removed preprocessor directive __chem again (2017-09-14, forkel)
    4044//
    41 //added phot
    42 //change of some output to lowercase with uppercase Fortran
    43 //
    44 // Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
    45 //
    46 
     45// Added phot                                                 (2017-09-14, forkel)
     46// --> Change in module_header: Variables for photolyis added (2017-09-14, forkel)
     47//
     48// change of some output to lowercase with uppercase Fortran (2017, forkel)
     49//
     50// Intial version of KP4 adapted to PALM                      (Nov. 2016, ketelsen)
     51//
    4752#include <stdio.h>
    48 //  mz_rs_20090111+
    4953// stdlib is necessary to define getenv:
    5054#include <stdlib.h>
    51 //  mz_rs_20090111-
    5255
    5356#include "create_kpp_module.h"
    … …  
    131134     }
    132135
    133 //   Change global variables to vector
     136//   Change global variables to vector (except for kp4_compress, which has already the right form)
    134137     
    135138     for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) {
    136        it->global_variables2vector (global_variable_list);
     139       if(it->get_name() != "kco_compress" ) {
     140         it->global_variables2vector (global_variable_list);
     141       }
    137142     }
    138143
    … …  
    162167       }
    163168     }
    164 // cout << "## after Edit individual subroutines    " <<endl;
    165 
    166169   }
    167170

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.h

    r3298 r3458  
    11#ifndef mecca
     2#define mecca 1
     3//
     4// ############################################################################
     5//
     6//     create_mz_kpp_module
     7//
     8//     create scalar code from .f90 sources created by KPP to be used in MESSy
     9//
     10//     COPYRIGHT Klaus Ketelsen and MPI-CH   April 2007
     11//
     12// ############################################################################
     13//
     14//Current revisions:
     15//------------------
     16//
     17//
     18// Former revisions:
     19// -----------------------
     20// $Id$
    221// ketelsen 18.09.2018: Removed create_fill_routine
    3 
    4 #include <iostream>
    5 #include <fstream>
    6 
    7 #include <string>
    8 #include <list>
    9 #include <vector>
    10 
    11 #include "fortran_file.h"
    12 #include "expand_decomp.h"
    13 
    14 // Class to create module which contains code generated by kpp.
    15 
    16 /
    17 /
    18 // ketelsen 18.09.2018: Removed create_fill_routine
     22//
     23// initial version                                  (Nov. 2016, ketelsen)
     24//
    1925
    2026#include <iostream>

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/expand_decomp.C

    r3298 r3458  
    99//
    1010// ############################################################################
     11//Current revisions:
     12//------------------
     13//
     14//
     15// Former revisions:
     16// -----------------------
     17// $Id: expand_decomp.C 3327 2018-10-09 19:55:00Z forkel $
     18//
     19// initial version                                  (Nov. 2016, ketelsen)
     20
    1121
    1222#include <fstream>

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/expand_decomp.h

    r2696 r3458  
    1111//
    1212// ############################################################################
     13//Current revisions:
     14//------------------
     15//
     16//
     17// Former revisions:
     18// -----------------------
     19// $Id: expand_decomp.h 3327 2018-10-09 19:55:00Z forkel $
     20//
     21// initial version                                  (Nov. 2016, ketelsen)
     22
    1323
    1424#include <iostream>

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.C

    r3298 r3458  
    1313//Current revisions:
    1414//------------------
    15 /
    16 /
     15//
     16//
    1717//Former revisions:
    1818//-----------------
    19 //$Id: fortran_file.C 2470 2017-09-14 13:56:42Z forkel $
    20 //  ketelsen 18.09.2018: (Rev 3260) added vector switch and creation of dimension statement
    21 //  forkel Sept.2018:  'fill_' routine applied to qvap and fakt
    22 //                      removal of unnecessary variables (LOOKAT, monitor etc. from Fortran code)
    23 //  forkel June 2018:  Moved adaption to PALM conventions to the end of the process
    24 //                     in order to make future use of vector code adaptations possible
    25 //  forkel 20.04.2018: Replace Roundoff = WLAMCH('E') by Roundoff = epsilon(one)
     19//$Id:
     20//
     21// Added vector switch and creation of dimension statement (rev. 3260, 18.09.2018, ketelsen)
    2622//
     23// removal of unnecessary variables (LOOKAT, monitor etc. from Fortran code) (Sept.2018, forkel)
    2724//
     25// Replace Roundoff = WLAMCH('E') by Roundoff = epsilon(one)   (20.04.2018, forkel)
    2826//
    29 / 2017-09-14 13:56:42Z forkel $
     27// Moved adaption to PALM conventions to the end of the processing kpp output
     28//  in order to make future use of vector code adaptations possible (June 2018, forkel)
     29//
     30// changed KPP-generated code to lowercase with uppercase Fortran  expressions
     31// added photolysis variables                               (2017-09-14, forkel)
    3032//
    31 //
    32 //changed KPP-generated code to lowercase with uppercase Fortran  expressions
    33 //added photolysis variables
    34 //
    35 //
    36 //
    37 //
    38 // Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
     33// Initial version of KP4 adapted to PALM                   (Nov 2016, ketelsen)
    3934//
    4035

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.h

    r3298 r3458  
    1414//Current revisions:
    1515//------------------
    16 /
    17 /
    18 //Former revisions:
    19 //-----------------
    20 // Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
    21 //  forkel Sept.2018:  added edit_Initialize
    22 //  forkel          :  added global_subtolower
    2316//
    24 / Intial version of KP4
    25 
     17//
     18// Former revisions:
     19// -----------------------
     20// $Id: fortran_file.h 3327 2018-10-09 19:55:00Z forkel $
     21//
     22// added edit_Initialize                            (Sept.2018, forkel)
     23//
     24// added global_subtolower                          (forkel)
     25// initial version                                  (Nov. 2016, ketelsen)
    2626
    2727#include <iostream>

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file_vec.C

    r3298 r3458  
    99//
    1010// ############################################################################
     11//
    1112//Current revisions:
    1213//------------------
    13 / 
    14 / 
     14//
     15//
    1516//Former revisions:
    16 //-----------------
    17 //$Id: fortran_file.C 2470 2017-09-14 13:56:42Z forkel $
    18 //  ketelsen 18.09.2018: Line 112: do k=is,ie bydo k=1,vl; line 156 ff: replaced index k by j
    19 //
    20 //  forkel Sept.2018:  added edit_Initialize
    21 //                     changed loop direction for update_rconst in edit_Update_RCONST
    22 / 2017-09-14 13:56:42Z forkel $
    23 //
    24 //
    25 // Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
    26 //
    27 
     17//-----------------------
     18//$Id: fortran_file_vec.C 3327 2018-10-09 19:55:00Z forkel $
     19// Line 112: do k=is,ie bydo k=1,vl; line 156 ff: replaced index k by j (18.09.2018, ketelsen)
     20//
     21// added edit_Initialize, changed loop direction for update_rconst in edit_Update_RCONST (Sept. 2018, forkel)
     22//
     23// initial version       (Nov. 2016, ketelsen)
     24//
    2825
    2926#include <fstream>

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/main.C

    r2696 r3458  
    99//
    1010// ############################################################################
     11//
     12//Current revisions:
     13//------------------
     14//
     15//
     16// Former revisions:
     17// -----------------------
     18// $Id: main.C 3327 2018-10-09 19:55:00Z forkel $
     19//
     20// initial version       (Nov. 2016, ketelsen)
     21//
    1122
    1223//  mz_rs_20090111+

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.C

    r3298 r3458  
    99//
    1010// ############################################################################
     11//
     12//Current revisions:
     13//------------------
     14//
     15//
     16// Former revisions:
     17// -----------------------
     18// $Id: program_line.C 3327 2018-10-09 19:55:00Z forkel $
     19//
     20// global_subtolower     (June 2018, forkel)
     21//
     22// initial version       (Nov. 2016, ketelsen)
     23//
     24
    1125
    1226#include "program_line.h"

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.h

    r3298 r3458  
    1212// ############################################################################
    1313//
    14 // forkel June 2018: added    void   global_subtolower(string &line);
     14//Current revisions:
     15//------------------
     16//
     17//
     18// Former revisions:
     19// -----------------------
     20// $Id: program_line.h 3327 2018-10-09 19:55:00Z forkel $
     21//
     22// added    void   global_subtolower(string &line); (June 2018, forkel)
     23//
     24// initial version                                  (Nov. 2016, ketelsen)
    1525
    1626#include <iostream>

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/utils.C

    r2696 r3458  
    99//
    1010// ############################################################################
     11//
     12//Current revisions:
     13//------------------
     14//
     15//
     16// Former revisions:
     17// -----------------------
     18// $Id: utils.C 3327 2018-10-09 19:55:00Z forkel $
     19//
     20// initial version                                  (Nov. 2016, ketelsen)
     21
    1122
    1223//  mz_rs_20090111+

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/utils.h

    r2696 r3458  
    1111//
    1212// ############################################################################
     13//
     14//Current revisions:
     15//------------------
     16//
     17//
     18// Former revisions:
     19// -----------------------
     20// $Id: utils.h 3327 2018-10-09 19:55:00Z forkel $
     21//
     22// initial version                                  (Nov. 2016, ketelsen)
     23//
    1324
    1425#include <string>

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/Rosenbrock_vec.f90

    r3298 r3458  
    22!       and can only be used with KP4.
    33!
     4!Current revisions:
     5!------------------
     6!
     7!
     8! Former revisions:
     9! -----------------------
     10! $Id$
     11! commented USE kp4_compress                       (30.10.2018, forkel)
     12!
     13! initial version (Rev. 3185)                      (June 2018, ketelsen)
     14!
     15
    416SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, &
    517           AbsTol,RelTol,              &
    … …  
    358370!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    359371
    360   USE kp4_compress,    ONLY: kco_initialize, kco_compress, &
    361                              kco_finalize, cell_done
     372! USE kp4_compress,    ONLY: kco_initialize, kco_compress, &
     373!                            kco_finalize, cell_done
    362374  IMPLICIT NONE
    363375

  • palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header

    r3298 r3458  
    4040!--------------------------------------------------------------------------------!
    4141!
    42 ! Current revisions:
    43 ! ------------------
    44 !
    45 !
    46 ! Former revisions:
    47 ! -----------------
    48 ! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
    49 ! forkel June 2018: qvap, fakt added
    50 ! forkel June 2018: reset case in  Initialize, Integrate, Update_rconst
    51 !
    52 !
    53 ! 2460 2017-09-13 14:47:48Z forkel
    5442!
    55 ! forkel Sept. 2017: Variables for photolyis added
     43! MODULE HEADER TEMPLATE
    5644!
    57 !
    58 ! Nov. 2016: Intial version (Klaus Ketelsen)
    59 !
    60 !------------------------------------------------------------------------------!
    61 !
    62 
     45!  Initial version (Nov. 2016, ketelsen), for later modifications of module_header
     46!  see comments in kpp4palm/src/create_kpp_module.C
    6347
    6448! Set kpp Double Precision to PALM Default Precision

  • palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_cbm4/chem_gasphase_mod.f90

    r3298 r3458  
    4444!--------------------------------------------------------------------------------!
    4545!
    46 ! Current revisions:
    47 ! ------------------
    48 !
    49 !
    50 ! Former revisions:
    51 ! -----------------
    52 ! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
    53 ! forkel June 2018: qvap,fakt added
    54 ! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    55 !
    56 !
    57 ! 2460 2017-09-13 14:47:48Z forkel
    5846!
    59 ! forkel Sept. 2017: Variables for photolyis added
     47! MODULE HEADER TEMPLATE
    6048!
    61 !
    62 ! Nov. 2016: Intial version (Klaus Ketelsen)
    63 !
    64 !------------------------------------------------------------------------------!
    65 !
    66 
     49!  Initial version (Nov. 2016,ketelsen),for later modifications of module_header
     50!  see comments in kpp4palm/src/create_kpp_module.C
    6751
    6852! Set kpp Double Precision to PALM Default Precision
    … …  
    129113!
    130114! File                 : chem_gasphase_mod_Parameters.f90
    131 ! Time                 : Tue Sep 25 18:35:13 2018
    132 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
     115! Time                 : Tue Oct 30 14:03:29 2018
     116! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    133117! Equation file        : chem_gasphase_mod.kpp
    134118! Output root filename : chem_gasphase_mod
    … …  
    236220!
    237221! File                 : chem_gasphase_mod_Global.f90
    238 ! Time                 : Tue Sep 25 18:35:13 2018
    239 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
     222! Time                 : Tue Oct 30 14:03:29 2018
     223! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    240224! Equation file        : chem_gasphase_mod.kpp
    241225! Output root filename : chem_gasphase_mod
    … …  
    310294!
    311295! File                 : chem_gasphase_mod_JacobianSP.f90
    312 ! Time                 : Tue Sep 25 18:35:13 2018
    313 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
     296! Time                 : Tue Oct 30 14:03:29 2018
     297! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    314298! Equation file        : chem_gasphase_mod.kpp
    315299! Output root filename : chem_gasphase_mod
    … …  
    406390!
    407391! File                 : chem_gasphase_mod_Monitor.f90
    408 ! Time                 : Tue Sep 25 18:35:13 2018
    409 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
     392! Time                 : Tue Oct 30 14:03:29 2018
     393! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    410394! Equation file        : chem_gasphase_mod.kpp
    411395! Output root filename : chem_gasphase_mod
    … …  
    569553!
    570554! File                 : chem_gasphase_mod_Initialize.f90
    571 ! Time                 : Tue Sep 25 18:35:13 2018
    572 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
     555! Time                 : Tue Oct 30 14:03:29 2018
     556! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    573557! Equation file        : chem_gasphase_mod.kpp
    574558! Output root filename : chem_gasphase_mod
    … …  
    595579!
    596580! File                 : chem_gasphase_mod_Integrator.f90
    597 ! Time                 : Tue Sep 25 18:35:13 2018
    598 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
     581! Time                 : Tue Oct 30 14:03:29 2018
     582! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    599583! Equation file        : chem_gasphase_mod.kpp
    600584! Output root filename : chem_gasphase_mod
    … …  
    653637!
    654638! File                 : chem_gasphase_mod_LinearAlgebra.f90
    655 ! Time                 : Tue Sep 25 18:35:13 2018
    656 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
     639! Time                 : Tue Oct 30 14:03:29 2018
     640! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    657641! Equation file        : chem_gasphase_mod.kpp
    658642! Output root filename : chem_gasphase_mod
    … …  
    680664!
    681665! File                 : chem_gasphase_mod_Jacobian.f90
    682 ! Time                 : Tue Sep 25 18:35:13 2018
    683 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
     666! Time                 : Tue Oct 30 14:03:29 2018
     667! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    684668! Equation file        : chem_gasphase_mod.kpp
    685669! Output root filename : chem_gasphase_mod
    … …  
    707691!
    708692! File                 : chem_gasphase_mod_Function.f90
    709 ! Time                 : Tue Sep 25 18:35:13 2018
    710 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
     693! Time                 : Tue Oct 30 14:03:29 2018
     694! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    711695! Equation file        : chem_gasphase_mod.kpp
    712696! Output root filename : chem_gasphase_mod
    … …  
    736720!
    737721! File                 : chem_gasphase_mod_Rates.f90
    738 ! Time                 : Tue Sep 25 18:35:13 2018
    739 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
     722! Time                 : Tue Oct 30 14:03:29 2018
     723! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    740724! Equation file        : chem_gasphase_mod.kpp
    741725! Output root filename : chem_gasphase_mod
    … …  
    762746!
    763747! File                 : chem_gasphase_mod_Util.f90
    764 ! Time                 : Tue Sep 25 18:35:13 2018
    765 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
     748! Time                 : Tue Oct 30 14:03:29 2018
     749! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    766750! Equation file        : chem_gasphase_mod.kpp
    767751! Output root filename : chem_gasphase_mod

  • palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive/chem_gasphase_mod.f90

    r3298 r3458  
    4646! Current revisions:
    4747! ------------------
     48! forkel June 2018: qvap,fakt added
     49! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    4850!
    4951!
    … …  
    5153! -----------------
    5254! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
    53 ! forkel June 2018: qvap,fakt added
    54 ! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    55 !
    56 !
    57 ! 2460 2017-09-13 14:47:48Z forkel
    5855!
    5956! forkel Sept. 2017: Variables for photolyis added

  • palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive1/chem_gasphase_mod.f90

    r3298 r3458  
    4646! Current revisions:
    4747! ------------------
     48! forkel June 2018: qvap,fakt added
     49! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    4850!
    4951!
    … …  
    5153! -----------------
    5254! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
    53 ! forkel June 2018: qvap,fakt added
    54 ! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    55 !
    56 !
    57 ! 2460 2017-09-13 14:47:48Z forkel
    5855!
    5956! forkel Sept. 2017: Variables for photolyis added

  • palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_phstat/chem_gasphase_mod.f90

    r3298 r3458  
    4646! Current revisions:
    4747! ------------------
     48! forkel June 2018: qvap,fakt added
     49! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    4850!
    4951!
    … …  
    5153! -----------------
    5254! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
    53 ! forkel June 2018: qvap,fakt added
    54 ! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    55 !
    56 !
    57 ! 2460 2017-09-13 14:47:48Z forkel
    5855!
    5956! forkel Sept. 2017: Variables for photolyis added

  • palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.f90

    r3298 r3458  
    4646! Current revisions:
    4747! ------------------
     48! forkel June 2018: qvap,fakt added
     49! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    4850!
    4951!
    … …  
    5153! -----------------
    5254! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
    53 ! forkel June 2018: qvap,fakt added
    54 ! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    55 !
    56 !
    57 ! 2460 2017-09-13 14:47:48Z forkel
    5855!
    5956! forkel Sept. 2017: Variables for photolyis added

  • palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simplep/chem_gasphase_mod.f90

    r3298 r3458  
    4646! Current revisions:
    4747! ------------------
     48! forkel June 2018: qvap,fakt added
     49! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    4850!
    4951!
    … …  
    5153! -----------------
    5254! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
    53 ! forkel June 2018: qvap,fakt added
    54 ! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    55 !
    56 !
    57 ! 2460 2017-09-13 14:47:48Z forkel
    5855!
    5956! forkel Sept. 2017: Variables for photolyis added

  • palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_smog/chem_gasphase_mod.f90

    r3298 r3458  
    4646! Current revisions:
    4747! ------------------
     48! forkel June 2018: qvap,fakt added
     49! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    4850!
    4951!
    … …  
    5153! -----------------
    5254! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
    53 ! forkel June 2018: qvap,fakt added
    54 ! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
    55 !
    56 !
    57 ! 2460 2017-09-13 14:47:48Z forkel
    5855!
    5956! forkel Sept. 2017: Variables for photolyis added

  • palm/trunk/UTIL/chemistry/gasphase_preproc/run_kpp4palm.ksh

    r3298 r3458  
    2424# Current revisions:
    2525# ------------------
    26 #
     26# Some correctetion of comments
    2727#
    2828# Former revisions:
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |