!> @file chemistry_model_mod.f90 !------------------------------------------------------------------------------! ! This file is part of the PALM model system. ! ! PALM is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 2017-2018 Leibniz Universitaet Hannover ! Copyright 2017-2018 Karlsruhe Institute of Technology ! Copyright 2017-2018 Freie Universitaet Berlin !------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: chemistry_model_mod.f90 3543 2018-11-20 17:06:15Z suehring $ ! Remove tabs ! ! 3542 2018-11-20 17:04:13Z suehring ! working precision added to make code Fortran 2008 conform ! ! 3458 2018-10-30 14:51:23Z kanani ! from chemistry branch r3443, banzhafs, basit: ! replace surf_lsm_h%qv1(m) by q(k,j,i) for mixing ratio in chem_depo ! bug fix in chem_depo: allow different surface fractions for one ! surface element and set lai to zero for non vegetated surfaces ! bug fixed in chem_data_output_2d ! bug fix in chem_depo subroutine ! added code on deposition of gases and particles ! removed cs_profile_name from chem_parin ! bug fixed in output profiles and code cleaned ! ! 3449 2018-10-29 19:36:56Z suehring ! additional output - merged from branch resler ! ! 3438 2018-10-28 19:31:42Z pavelkrc ! Add terrain-following masked output ! ! 3373 2018-10-18 15:25:56Z kanani ! Remove MPI_Abort, replace by message ! ! 3318 2018-10-08 11:43:01Z sward ! Fixed faulty syntax of message string ! ! 3298 2018-10-02 12:21:11Z kanani ! Add remarks (kanani) ! Merge with trunk, replaced cloud_physics by bulk_cloud_model 28.09.2018 forkel ! Subroutines header and chem_check_parameters added 25.09.2018 basit ! Removed chem_emission routine now declared in chem_emissions.f90 30.07.2018 ERUSSO ! Introduced emissions namelist parameters 30.07.2018 ERUSSO ! ! Timestep steering added in subroutine chem_integrate_ij and ! output of chosen solver in chem_parin added 30.07.2018 ketelsen ! ! chem_check_data_output_pr: added unit for PM compounds 20.07.2018 forkel ! replaced : by nzb+1:nzt for pt,q,ql (found by kk) 18.07.2018 forkel ! debugged restart run for chem species 06.07.2018 basit ! reorganized subroutines in alphabetical order. 27.06.2018 basit ! subroutine chem_parin updated for profile output 27.06.2018 basit ! Added humidity arrays to USE section and tmp_qvap in chem_integrate 26.6.2018 forkel ! Merged chemistry with with trunk (nzb_do and nzt_do in 3d output) 26.6.2018 forkel ! ! reorganized subroutines in alphabetical order. basit 22.06.2018 ! subroutine chem_parin updated for profile output basit 22.06.2018 ! subroutine chem_statistics added ! subroutine chem_check_data_output_pr add 21.06.2018 basit ! subroutine chem_data_output_mask added 20.05.2018 basit ! subroutine chem_data_output_2d added 20.05.2018 basit ! subroutine chem_statistics added 04.06.2018 basit ! subroutine chem_emissions: Set cssws to zero before setting values 20.03.2018 forkel ! subroutine chem_emissions: Introduced different conversion factors ! for PM and gaseous compounds 15.03.2018 forkel ! subroutine chem_emissions updated to take variable number of chem_spcs and ! emission factors. 13.03.2018 basit ! chem_boundary_conds_decycle improved. 05.03.2018 basit ! chem_boundary_conds_decycle subroutine added 21.02.2018 basit ! chem_init_profiles subroutines re-activated after correction 21.02.2018 basit ! ! ! 3293 2018-09-28 12:45:20Z forkel ! Modularization of all bulk cloud physics code components ! ! 3248 2018-09-14 09:42:06Z sward ! Minor formating changes ! ! 3246 2018-09-13 15:14:50Z sward ! Added error handling for input namelist via parin_fail_message ! ! 3241 2018-09-12 15:02:00Z raasch ! +nest_chemistry ! ! 3209 2018-08-27 16:58:37Z suehring ! Rename flags indicating outflow boundary conditions ! ! 3182 2018-07-27 13:36:03Z suehring ! Revise output of surface quantities in case of overhanging structures ! ! 3045 2018-05-28 07:55:41Z Giersch ! error messages revised ! ! 3014 2018-05-09 08:42:38Z maronga ! Bugfix: nzb_do and nzt_do were not used for 3d data output ! ! 3004 2018-04-27 12:33:25Z Giersch ! Comment concerning averaged data output added ! ! 2932 2018-03-26 09:39:22Z maronga ! renamed chemistry_par to chemistry_parameters ! ! 2894 2018-03-15 09:17:58Z Giersch ! Calculations of the index range of the subdomain on file which overlaps with ! the current subdomain are already done in read_restart_data_mod, ! chem_last_actions was renamed to chem_wrd_local, chem_read_restart_data was ! renamed to chem_rrd_local, chem_write_var_list was renamed to ! chem_wrd_global, chem_read_var_list was renamed to chem_rrd_global, ! chem_skip_var_list has been removed, variable named found has been ! introduced for checking if restart data was found, reading of restart strings ! has been moved completely to read_restart_data_mod, chem_rrd_local is already ! inside the overlap loop programmed in read_restart_data_mod, todo list has ! bees extended, redundant characters in chem_wrd_local have been removed, ! the marker *** end chemistry *** is not necessary anymore, strings and their ! respective lengths are written out and read now in case of restart runs to ! get rid of prescribed character lengths ! ! 2815 2018-02-19 11:29:57Z suehring ! Bugfix in restart mechanism, ! rename chem_tendency to chem_prognostic_equations, ! implement vector-optimized version of chem_prognostic_equations, ! some clean up (incl. todo list) ! ! 2773 2018-01-30 14:12:54Z suehring ! Declare variables required for nesting as public ! ! 2772 2018-01-29 13:10:35Z suehring ! Bugfix in string handling ! ! 2768 2018-01-24 15:38:29Z kanani ! Shorten lines to maximum length of 132 characters ! ! 2766 2018-01-22 17:17:47Z kanani ! Removed preprocessor directive __chem ! ! 2756 2018-01-16 18:11:14Z suehring ! Fill values in 3D output introduced. ! ! 2718 2018-01-02 08:49:38Z maronga ! Initial revision ! ! ! ! ! Authors: ! -------- ! @author Renate Forkel ! @author Farah Kanani-Suehring ! @author Klaus Ketelsen ! @author Basit Khan ! @author Sabine Banzhaf ! ! !------------------------------------------------------------------------------! ! Description: ! ------------ !> Chemistry model for PALM-4U !> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not !> allowed to use the chemistry model in a precursor run and additionally !> not using it in a main run !> @todo Update/clean-up todo list! (FK) !> @todo Set proper fill values (/= 0) for chem output arrays! (FK) !> @todo Add routine chem_check_parameters, add checks for inconsistent or !> unallowed parameter settings. !> CALL of chem_check_parameters from check_parameters. (FK) !> @todo Make routine chem_header available, CALL from header.f90 !> (see e.g. how it is done in routine lsm_header in !> land_surface_model_mod.f90). chem_header should include all setup !> info about chemistry parameter settings. (FK) !> @todo Implement turbulent inflow of chem spcs in inflow_turbulence. !> @todo Separate boundary conditions for each chem spcs to be implemented !> @todo Currently only total concentration are calculated. Resolved, parameterized !> and chemistry fluxes although partially and some completely coded but !> are not operational/activated in this version. bK. !> @todo slight differences in passive scalar and chem spcs when chem reactions !> turned off. Need to be fixed. bK !> @todo test nesting for chem spcs, was implemented by suehring (kanani) !> @todo chemistry error messages !> @todo Format this module according to PALM coding standard (see e.g. module !> template under http://palm.muk.uni-hannover.de/mosaik/downloads/8 or !> D3_coding_standard.pdf under https://palm.muk.uni-hannover.de/trac/downloads/16) ! !------------------------------------------------------------------------------! MODULE chemistry_model_mod USE kinds, ONLY: wp, iwp USE indices, ONLY: nz, nzb,nzt,nysg,nyng,nxlg,nxrg, & nys,nyn,nx,nxl,nxr,ny,wall_flags_0 USE pegrid, ONLY: myid, threads_per_task USE bulk_cloud_model_mod, & ONLY: bulk_cloud_model USE control_parameters, ONLY: bc_lr_cyc, bc_ns_cyc, dt_3d, humidity, & initializing_actions, message_string, & omega, tsc, intermediate_timestep_count, & intermediate_timestep_count_max, max_pr_user, & timestep_scheme, use_prescribed_profile_data, ws_scheme_sca USE arrays_3d, ONLY: exner, hyp, pt, q, ql, rdf_sc, tend, zu USE chem_gasphase_mod, ONLY: nspec, spc_names, nkppctrl, nmaxfixsteps, & t_steps, chem_gasphase_integrate, vl_dim, & nvar, nreact, atol, rtol, nphot, phot_names USE cpulog, ONLY: cpu_log, log_point USE chem_modules USE statistics IMPLICIT NONE PRIVATE SAVE LOGICAL :: nest_chemistry = .TRUE. !< flag for nesting mode of chemical species, independent on parent or not ! !- Define chemical variables TYPE species_def CHARACTER(LEN=8) :: name CHARACTER(LEN=16) :: unit REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc_av REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc_p REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: tconc_m REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: cssws_av REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: flux_s_cs REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: diss_s_cs REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:) :: flux_l_cs REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:) :: diss_l_cs REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: conc_pr_init END TYPE species_def TYPE photols_def CHARACTER(LEN=8) :: name CHARACTER(LEN=16) :: unit REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: freq END TYPE photols_def PUBLIC species_def PUBLIC photols_def TYPE(species_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC :: chem_species TYPE(photols_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC :: phot_frequen REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_1 REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_2 REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_3 REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_av REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: freq_1 INTEGER,DIMENSION(nkppctrl) :: icntrl ! Fine tuning kpp REAL(kind=wp),DIMENSION(nkppctrl) :: rcntrl ! Fine tuning kpp LOGICAL :: decycle_chem_lr = .FALSE. LOGICAL :: decycle_chem_ns = .FALSE. CHARACTER (LEN=20), DIMENSION(4) :: decycle_method = & (/'dirichlet','dirichlet','dirichlet','dirichlet'/) !< Decycling method at horizontal boundaries, !< 1=left, 2=right, 3=south, 4=north !< dirichlet = initial size distribution and !< chemical composition set for the ghost and !< first three layers !< neumann = zero gradient REAL(kind=wp), PUBLIC :: cs_time_step = 0._wp CHARACTER(10), PUBLIC :: photolysis_scheme ! 'constant', ! 'simple' (Simple parameterisation from MCM, Saunders et al., 2003, Atmos. Chem. Phys., 3, 161-180 ! 'fastj' (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED !----------------------------------------------------------------------- ! Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010) ! INTEGER(iwp), PARAMETER :: nlu_dep = 15 ! Number of DEPAC landuse classes (lu's) INTEGER(iwp), PARAMETER :: ncmp = 10 ! Number of DEPAC gas components !------------------------------------------------------------------------ ! DEPAC landuse classes as defined in LOTOS-EUROS model v2.1 ! INTEGER(iwp) :: ilu_grass = 1 INTEGER(iwp) :: ilu_arable = 2 INTEGER(iwp) :: ilu_permanent_crops = 3 INTEGER(iwp) :: ilu_coniferous_forest = 4 INTEGER(iwp) :: ilu_deciduous_forest = 5 INTEGER(iwp) :: ilu_water_sea = 6 INTEGER(iwp) :: ilu_urban = 7 INTEGER(iwp) :: ilu_other = 8 INTEGER(iwp) :: ilu_desert = 9 INTEGER(iwp) :: ilu_ice = 10 INTEGER(iwp) :: ilu_savanna = 11 INTEGER(iwp) :: ilu_tropical_forest = 12 INTEGER(iwp) :: ilu_water_inland = 13 INTEGER(iwp) :: ilu_mediterrean_scrub = 14 INTEGER(iwp) :: ilu_semi_natural_veg = 15 ! !-------------------------------------------------------------------------- ! NH3/SO2 ratio regimes: INTEGER, PARAMETER :: iratns_low = 1 ! low ratio NH3/SO2 INTEGER, PARAMETER :: iratns_high = 2 ! high ratio NH3/SO2 INTEGER, PARAMETER :: iratns_very_low = 3 ! very low ratio NH3/SO2 ! Default: INTEGER, PARAMETER :: iratns_default = iratns_low ! ! Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) naar W m-2 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 ! ! Set temperatures per land use for F_temp 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/) 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 /) 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 /) ! ! Set F_min: 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/) ! Set maximal conductance (m/s) ! (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from mmol O3/m2/s to m/s ! Could be refined to a function of T and P. in Jones REAL(wp), DIMENSION(nlu_dep), PARAMETER :: g_max =(/270., 300., 300., 140., 150., -999., -999.,270., -999.,-999.,270., 150.,-999.,300., 422./)/41000 ! ! Set max, min for vapour pressure deficit vpd; 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/) 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/) ! ! !------------------------------------------------------------------------ PUBLIC nest_chemistry PUBLIC nspec PUBLIC nreact PUBLIC nvar PUBLIC spc_names PUBLIC spec_conc_2 !- Interface section INTERFACE chem_3d_data_averaging MODULE PROCEDURE chem_3d_data_averaging END INTERFACE chem_3d_data_averaging INTERFACE chem_boundary_conds MODULE PROCEDURE chem_boundary_conds MODULE PROCEDURE chem_boundary_conds_decycle END INTERFACE chem_boundary_conds INTERFACE chem_check_data_output MODULE PROCEDURE chem_check_data_output END INTERFACE chem_check_data_output INTERFACE chem_data_output_2d MODULE PROCEDURE chem_data_output_2d END INTERFACE chem_data_output_2d INTERFACE chem_data_output_3d MODULE PROCEDURE chem_data_output_3d END INTERFACE chem_data_output_3d INTERFACE chem_data_output_mask MODULE PROCEDURE chem_data_output_mask END INTERFACE chem_data_output_mask INTERFACE chem_check_data_output_pr MODULE PROCEDURE chem_check_data_output_pr END INTERFACE chem_check_data_output_pr INTERFACE chem_check_parameters MODULE PROCEDURE chem_check_parameters END INTERFACE chem_check_parameters INTERFACE chem_define_netcdf_grid MODULE PROCEDURE chem_define_netcdf_grid END INTERFACE chem_define_netcdf_grid INTERFACE chem_header MODULE PROCEDURE chem_header END INTERFACE chem_header INTERFACE chem_init MODULE PROCEDURE chem_init END INTERFACE chem_init INTERFACE chem_init_profiles MODULE PROCEDURE chem_init_profiles END INTERFACE chem_init_profiles INTERFACE chem_integrate MODULE PROCEDURE chem_integrate_ij END INTERFACE chem_integrate INTERFACE chem_parin MODULE PROCEDURE chem_parin END INTERFACE chem_parin INTERFACE chem_prognostic_equations MODULE PROCEDURE chem_prognostic_equations MODULE PROCEDURE chem_prognostic_equations_ij END INTERFACE chem_prognostic_equations INTERFACE chem_rrd_local MODULE PROCEDURE chem_rrd_local END INTERFACE chem_rrd_local INTERFACE chem_statistics MODULE PROCEDURE chem_statistics END INTERFACE chem_statistics INTERFACE chem_swap_timelevel MODULE PROCEDURE chem_swap_timelevel END INTERFACE chem_swap_timelevel INTERFACE chem_wrd_local MODULE PROCEDURE chem_wrd_local END INTERFACE chem_wrd_local INTERFACE chem_depo MODULE PROCEDURE chem_depo END INTERFACE chem_depo INTERFACE drydepos_gas_depac MODULE PROCEDURE drydepos_gas_depac END INTERFACE drydepos_gas_depac INTERFACE rc_special MODULE PROCEDURE rc_special END INTERFACE rc_special INTERFACE rc_gw MODULE PROCEDURE rc_gw END INTERFACE rc_gw INTERFACE rw_so2 MODULE PROCEDURE rw_so2 END INTERFACE rw_so2 INTERFACE rw_nh3_sutton MODULE PROCEDURE rw_nh3_sutton END INTERFACE rw_nh3_sutton INTERFACE rw_constant MODULE PROCEDURE rw_constant END INTERFACE rw_constant INTERFACE rc_gstom MODULE PROCEDURE rc_gstom END INTERFACE rc_gstom INTERFACE rc_gstom_emb MODULE PROCEDURE rc_gstom_emb END INTERFACE rc_gstom_emb INTERFACE par_dir_diff MODULE PROCEDURE par_dir_diff END INTERFACE par_dir_diff INTERFACE rc_get_vpd MODULE PROCEDURE rc_get_vpd END INTERFACE rc_get_vpd INTERFACE rc_gsoil_eff MODULE PROCEDURE rc_gsoil_eff END INTERFACE rc_gsoil_eff INTERFACE rc_rinc MODULE PROCEDURE rc_rinc END INTERFACE rc_rinc INTERFACE rc_rctot MODULE PROCEDURE rc_rctot END INTERFACE rc_rctot INTERFACE rc_comp_point_rc_eff MODULE PROCEDURE rc_comp_point_rc_eff END INTERFACE rc_comp_point_rc_eff INTERFACE drydepo_aero_zhang_vd MODULE PROCEDURE drydepo_aero_zhang_vd END INTERFACE drydepo_aero_zhang_vd INTERFACE get_rb_cell MODULE PROCEDURE get_rb_cell END INTERFACE get_rb_cell PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, & chem_check_data_output_pr, chem_check_parameters, & chem_data_output_2d, chem_data_output_3d, chem_data_output_mask, & chem_define_netcdf_grid, chem_header, chem_init, & chem_init_profiles, chem_integrate, chem_parin, & chem_prognostic_equations, chem_rrd_local, & chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo CONTAINS ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine for averaging 3D data of chemical species. Due to the fact that !> the averaged chem arrays are allocated in chem_init, no if-query concerning !> the allocation is required (in any mode). Attention: If you just specify an !> averaged output quantity in the _p3dr file during restarts the first output !> includes the time between the beginning of the restart run and the first !> output time (not necessarily the whole averaging_interval you have !> specified in your _p3d/_p3dr file ) !------------------------------------------------------------------------------! SUBROUTINE chem_3d_data_averaging ( mode, variable ) USE control_parameters USE indices USE kinds USE surface_mod, & ONLY: surf_def_h, surf_lsm_h, surf_usm_h IMPLICIT NONE CHARACTER (LEN=*) :: mode !< CHARACTER (LEN=*) :: variable !< LOGICAL :: match_def !< flag indicating natural-type surface LOGICAL :: match_lsm !< flag indicating natural-type surface LOGICAL :: match_usm !< flag indicating urban-type surface INTEGER(iwp) :: i !< grid index x direction INTEGER(iwp) :: j !< grid index y direction INTEGER(iwp) :: k !< grid index z direction INTEGER(iwp) :: m !< running index surface type INTEGER(iwp) :: lsp !< running index for chem spcs IF ( mode == 'allocate' ) THEN DO lsp = 1, nspec IF ( TRIM(variable(1:3)) == 'kc_' .AND. & TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN chem_species(lsp)%conc_av = 0.0_wp ENDIF ENDDO ELSEIF ( mode == 'sum' ) THEN DO lsp = 1, nspec IF ( TRIM(variable(1:3)) == 'kc_' .AND. & TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) + & chem_species(lsp)%conc(k,j,i) ENDDO ENDDO ENDDO ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') ) THEN DO i = nxl, nxr DO j = nys, nyn match_def = surf_def_h(0)%start_index(j,i) <= & surf_def_h(0)%end_index(j,i) match_lsm = surf_lsm_h%start_index(j,i) <= & surf_lsm_h%end_index(j,i) match_usm = surf_usm_h%start_index(j,i) <= & surf_usm_h%end_index(j,i) IF ( match_def ) THEN m = surf_def_h(0)%end_index(j,i) chem_species(lsp)%cssws_av(j,i) = & chem_species(lsp)%cssws_av(j,i) + & surf_def_h(0)%cssws(lsp,m) ELSEIF ( match_lsm .AND. .NOT. match_usm ) THEN m = surf_lsm_h%end_index(j,i) chem_species(lsp)%cssws_av(j,i) = & chem_species(lsp)%cssws_av(j,i) + & surf_lsm_h%cssws(lsp,m) ELSEIF ( match_usm ) THEN m = surf_usm_h%end_index(j,i) chem_species(lsp)%cssws_av(j,i) = & chem_species(lsp)%cssws_av(j,i) + & surf_usm_h%cssws(lsp,m) ENDIF ENDDO ENDDO ENDIF ENDDO ELSEIF ( mode == 'average' ) THEN DO lsp = 1, nspec IF ( TRIM(variable(1:3)) == 'kc_' .AND. & TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') ) THEN DO i = nxlg, nxrg DO j = nysg, nyng chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp ) ENDIF ENDDO ENDIF END SUBROUTINE chem_3d_data_averaging ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine to initialize and set all boundary conditions for chemical species !------------------------------------------------------------------------------! SUBROUTINE chem_boundary_conds( mode ) USE control_parameters, & ONLY: air_chemistry, bc_radiation_l, bc_radiation_n, bc_radiation_r, & bc_radiation_s USE indices, & ONLY: nxl, nxr, nxlg, nxrg, nyng, nysg, nzt USE arrays_3d, & ONLY: dzu USE surface_mod, & ONLY: bc_h CHARACTER (len=*), INTENT(IN) :: mode INTEGER(iwp) :: i !< grid index x direction. INTEGER(iwp) :: j !< grid index y direction. INTEGER(iwp) :: k !< grid index z direction. INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing. INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls. INTEGER(iwp) :: m !< running index surface elements. INTEGER(iwp) :: lsp !< running index for chem spcs. INTEGER(iwp) :: lph !< running index for photolysis frequencies SELECT CASE ( TRIM( mode ) ) CASE ( 'init' ) IF ( bc_cs_b == 'dirichlet' ) THEN ibc_cs_b = 0 ELSEIF ( bc_cs_b == 'neumann' ) THEN ibc_cs_b = 1 ELSE message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 ) !< ENDIF ! !-- Set Integer flags and check for possible erroneous settings for top !-- boundary condition. IF ( bc_cs_t == 'dirichlet' ) THEN ibc_cs_t = 0 ELSEIF ( bc_cs_t == 'neumann' ) THEN ibc_cs_t = 1 ELSEIF ( bc_cs_t == 'initial_gradient' ) THEN ibc_cs_t = 2 ELSEIF ( bc_cs_t == 'nested' ) THEN ibc_cs_t = 3 ELSE message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"' CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 ) ENDIF CASE ( 'set_bc_bottomtop' ) !-- Bottom boundary condtions for chemical species DO lsp = 1, nspec IF ( ibc_cs_b == 0 ) THEN DO l = 0, 1 !-- Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e. !-- the chem_species(nspec)%conc_p value at the topography top (k-1) !-- is set; for downward-facing surfaces (l=1), kb=1, i.e. the !-- value at the topography bottom (k+1) is set. kb = MERGE( -1, 1, l == 0 ) !$OMP PARALLEL DO PRIVATE( i, j, k ) DO m = 1, bc_h(l)%ns i = bc_h(l)%i(m) j = bc_h(l)%j(m) k = bc_h(l)%k(m) chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i) ENDDO ENDDO ELSEIF ( ibc_cs_b == 1 ) THEN !-- in boundary_conds there is som extra loop over m here for passive tracer DO l = 0, 1 kb = MERGE( -1, 1, l == 0 ) !$OMP PARALLEL DO PRIVATE( i, j, k ) DO m = 1, bc_h(l)%ns i = bc_h(l)%i(m) j = bc_h(l)%j(m) k = bc_h(l)%k(m) chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc_p(k,j,i) ENDDO ENDDO ENDIF ENDDO ! end lsp loop !-- Top boundary conditions for chemical species - Should this not be done for all species? IF ( ibc_cs_t == 0 ) THEN DO lsp = 1, nspec chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:) ENDDO ELSEIF ( ibc_cs_t == 1 ) THEN DO lsp = 1, nspec chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) ENDDO ELSEIF ( ibc_cs_t == 2 ) THEN DO lsp = 1, nspec chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1) ENDDO ENDIF ! CASE ( 'set_bc_lateral' ) !-- Lateral boundary conditions for chem species at inflow boundary !-- are automatically set when chem_species concentration is !-- initialized. The initially set value at the inflow boundary is not !-- touched during time integration, hence, this boundary value remains !-- at a constant value, which is the concentration that flows into the !-- domain. !-- Lateral boundary conditions for chem species at outflow boundary IF ( bc_radiation_s ) THEN DO lsp = 1, nspec chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:) ENDDO ELSEIF ( bc_radiation_n ) THEN DO lsp = 1, nspec chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:) ENDDO ELSEIF ( bc_radiation_l ) THEN DO lsp = 1, nspec chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl) ENDDO ELSEIF ( bc_radiation_r ) THEN DO lsp = 1, nspec chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr) ENDDO ENDIF END SELECT END SUBROUTINE chem_boundary_conds ! !------------------------------------------------------------------------------! ! Description: ! ------------ !> Boundary conditions for prognostic variables in chemistry: decycling in the !> x-direction !----------------------------------------------------------------------------- SUBROUTINE chem_boundary_conds_decycle (cs_3d, cs_pr_init ) USE pegrid, & ONLY: myid IMPLICIT NONE INTEGER(iwp) :: boundary !< INTEGER(iwp) :: ee !< INTEGER(iwp) :: copied !< INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< INTEGER(iwp) :: ss !< REAL(wp), DIMENSION(nzb:nzt+1) :: cs_pr_init REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: cs_3d REAL(wp) :: flag !< flag to mask topography grid points flag = 0.0_wp !-- Left and right boundaries IF ( decycle_chem_lr .AND. bc_lr_cyc ) THEN DO boundary = 1, 2 IF ( decycle_method(boundary) == 'dirichlet' ) THEN ! !-- Initial profile is copied to ghost and first three layers ss = 1 ee = 0 IF ( boundary == 1 .AND. nxl == 0 ) THEN ss = nxlg ee = nxl+2 ELSEIF ( boundary == 2 .AND. nxr == nx ) THEN ss = nxr-2 ee = nxrg ENDIF DO i = ss, ee DO j = nysg, nyng DO k = nzb+1, nzt flag = MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) ) cs_3d(k,j,i) = cs_pr_init(k) * flag ENDDO ENDDO ENDDO ELSEIF ( decycle_method(boundary) == 'neumann' ) THEN ! !-- The value at the boundary is copied to the ghost layers to simulate !-- an outlet with zero gradient ss = 1 ee = 0 IF ( boundary == 1 .AND. nxl == 0 ) THEN ss = nxlg ee = nxl-1 copied = nxl ELSEIF ( boundary == 2 .AND. nxr == nx ) THEN ss = nxr+1 ee = nxrg copied = nxr ENDIF DO i = ss, ee DO j = nysg, nyng DO k = nzb+1, nzt flag = MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) ) cs_3d(k,j,i) = cs_3d(k,j,copied) * flag ENDDO ENDDO ENDDO ELSE WRITE(message_string,*) & 'unknown decycling method: decycle_method (', & boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"' CALL message( 'chem_boundary_conds_decycle', 'CM0431', & 1, 2, 0, 6, 0 ) ENDIF ENDDO ENDIF !-- South and north boundaries IF ( decycle_chem_ns .AND. bc_ns_cyc ) THEN DO boundary = 3, 4 IF ( decycle_method(boundary) == 'dirichlet' ) THEN ! !-- Initial profile is copied to ghost and first three layers ss = 1 ee = 0 IF ( boundary == 3 .AND. nys == 0 ) THEN ss = nysg ee = nys+2 ELSEIF ( boundary == 4 .AND. nyn == ny ) THEN ss = nyn-2 ee = nyng ENDIF DO i = nxlg, nxrg DO j = ss, ee DO k = nzb+1, nzt flag = MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) ) cs_3d(k,j,i) = cs_pr_init(k) * flag ENDDO ENDDO ENDDO ELSEIF ( decycle_method(boundary) == 'neumann' ) THEN ! !-- The value at the boundary is copied to the ghost layers to simulate !-- an outlet with zero gradient ss = 1 ee = 0 IF ( boundary == 3 .AND. nys == 0 ) THEN ss = nysg ee = nys-1 copied = nys ELSEIF ( boundary == 4 .AND. nyn == ny ) THEN ss = nyn+1 ee = nyng copied = nyn ENDIF DO i = nxlg, nxrg DO j = ss, ee DO k = nzb+1, nzt flag = MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) ) cs_3d(k,j,i) = cs_3d(k,copied,i) * flag ENDDO ENDDO ENDDO ELSE WRITE(message_string,*) & 'unknown decycling method: decycle_method (', & boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"' CALL message( 'chem_boundary_conds_decycle', 'CM0432', & 1, 2, 0, 6, 0 ) ENDIF ENDDO ENDIF END SUBROUTINE chem_boundary_conds_decycle ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine for checking data output for chemical species !------------------------------------------------------------------------------! SUBROUTINE chem_check_data_output( var, unit, i, ilen, k ) USE control_parameters, & ONLY: data_output, message_string IMPLICIT NONE CHARACTER (LEN=*) :: unit !< CHARACTER (LEN=*) :: var !< INTEGER(iwp) :: i, lsp INTEGER(iwp) :: ilen INTEGER(iwp) :: k CHARACTER(len=16) :: spec_name unit = 'illegal' spec_name = TRIM(var(4:)) !< var 1:3 is 'kc_' or 'em_'. IF ( TRIM(var(1:3)) == 'em_' ) THEN DO lsp=1,nspec IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN unit = 'mol m-2 s-1' ENDIF !-- It is possible to plant PM10 and PM25 into the gasphase chemistry code !-- as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): !-- set unit to micrograms per m**3 for PM10 and PM25 (PM2.5) IF (spec_name(1:2) == 'PM') THEN unit = 'kg m-2 s-1' ENDIF ENDDO ELSE DO lsp=1,nspec IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN unit = 'ppm' ENDIF ! It is possible to plant PM10 and PM25 into the gasphase chemistry code ! as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): ! set unit to kilograms per m**3 for PM10 and PM25 (PM2.5) IF (spec_name(1:2) == 'PM') THEN unit = 'kg m-3' ENDIF ENDDO DO lsp=1,nphot IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name)) THEN unit = 'sec-1' ENDIF ENDDO ENDIF RETURN END SUBROUTINE chem_check_data_output ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine for checking data output of profiles for chemistry model !------------------------------------------------------------------------------! SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit ) USE arrays_3d USE control_parameters, & ONLY: data_output_pr, message_string, air_chemistry USE indices USE profil_parameter USE statistics IMPLICIT NONE CHARACTER (LEN=*) :: unit !< CHARACTER (LEN=*) :: variable !< CHARACTER (LEN=*) :: dopr_unit CHARACTER(len=16) :: spec_name INTEGER(iwp) :: var_count, lsp !< spec_name = TRIM(variable(4:)) IF ( .NOT. air_chemistry ) THEN message_string = 'data_output_pr = ' // & TRIM( data_output_pr(var_count) ) // ' is not imp' // & 'lemented for air_chemistry = .FALSE.' CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 ) ELSE DO lsp = 1, nspec IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN cs_pr_count = cs_pr_count+1 cs_pr_index(cs_pr_count) = lsp dopr_index(var_count) = pr_palm+cs_pr_count dopr_unit = 'ppm' IF (spec_name(1:2) == 'PM') THEN dopr_unit = 'kg m-3' ENDIF hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 ) unit = dopr_unit ENDIF ENDDO ENDIF END SUBROUTINE chem_check_data_output_pr ! !------------------------------------------------------------------------------! ! Description: ! ------------ !> Check parameters routine for chemistry_model_mod !------------------------------------------------------------------------------! SUBROUTINE chem_check_parameters IMPLICIT NONE LOGICAL :: found INTEGER (iwp) :: lsp_usr !< running index for user defined chem spcs INTEGER (iwp) :: lsp !< running index for chem spcs. !!-- check for chemical reactions status IF ( chem_gasphase_on ) THEN message_string = 'Chemical reactions: ON' CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 ) ELSEIF ( .not. (chem_gasphase_on) ) THEN message_string = 'Chemical reactions: OFF' CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 ) ENDIF !-- check for chemistry time-step IF ( call_chem_at_all_substeps ) THEN message_string = 'Chemistry is calculated at all meteorology time-step' CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 ) ELSEIF ( .not. (call_chem_at_all_substeps) ) THEN message_string = 'Sub-time-steps are skipped for chemistry time-steps' CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 ) ENDIF !-- check for photolysis scheme IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant') ) THEN message_string = 'Incorrect photolysis scheme selected, please check spelling' CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 ) ENDIF !-- check for decycling of chem species IF ( (.not. any(decycle_method == 'neumann') ) .AND. (.not. any(decycle_method == 'dirichlet') ) ) THEN message_string = 'Incorrect boundary conditions. Only neumann or ' & // 'dirichlet &available for decycling chemical species ' CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 ) ENDIF !--------------------- !-- chem_check_parameters is called before the array chem_species is allocated! !-- temporary switch of this part of the check RETURN !--------------------- !-- check for initial chem species input lsp_usr = 1 lsp = 1 DO WHILE ( cs_name (lsp_usr) /= 'novalue') found = .FALSE. DO lsp = 1, nvar IF ( TRIM(cs_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN found = .TRUE. EXIT ENDIF ENDDO IF ( .not. found ) THEN message_string = 'Incorrect input for initial surface vaue: ' // TRIM(cs_name(lsp_usr)) CALL message( 'chem_check_parameters', 'CM0427', 0, 1, 0, 6, 0 ) ENDIF lsp_usr = lsp_usr + 1 ENDDO !-- check for surface emission flux chem species lsp_usr = 1 lsp = 1 DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue') found = .FALSE. DO lsp = 1, nvar IF ( TRIM(surface_csflux_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN found = .TRUE. EXIT ENDIF ENDDO IF ( .not. found ) THEN message_string = 'Incorrect input of chemical species for surface emission fluxes: ' & // TRIM(surface_csflux_name(lsp_usr)) CALL message( 'chem_check_parameters', 'CM0428', 0, 1, 0, 6, 0 ) ENDIF lsp_usr = lsp_usr + 1 ENDDO END SUBROUTINE chem_check_parameters ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine defining 2D output variables for chemical species !> @todo: Remove "mode" from argument list, not used. !------------------------------------------------------------------------------! SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf, & two_d, nzb_do, nzt_do, fill_value ) USE indices USE kinds USE pegrid, ONLY: myid, threads_per_task IMPLICIT NONE CHARACTER (LEN=*) :: grid !< CHARACTER (LEN=*) :: mode !< CHARACTER (LEN=*) :: variable !< INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data INTEGER(iwp) :: nzb_do !< lower limit of the domain (usually nzb) INTEGER(iwp) :: nzt_do !< upper limit of the domain (usually nzt+1) LOGICAL :: found !< LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections) REAL(wp) :: fill_value REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) :: local_pf !< !-- local variables. CHARACTER(len=16) :: spec_name INTEGER(iwp) :: lsp INTEGER(iwp) :: i !< grid index along x-direction INTEGER(iwp) :: j !< grid index along y-direction INTEGER(iwp) :: k !< grid index along z-direction INTEGER(iwp) :: m !< running index surface elements INTEGER(iwp) :: char_len !< length of a character string found = .TRUE. char_len = LEN_TRIM(variable) spec_name = TRIM( variable(4:char_len-3) ) DO lsp=1,nspec IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name) .AND. & ( (variable(char_len-2:) == '_xy') .OR. & (variable(char_len-2:) == '_xz') .OR. & (variable(char_len-2:) == '_yz') ) ) THEN ! !-- todo: remove or replace by "CALL message" mechanism (kanani) ! IF(myid == 0) WRITE(6,*) 'Output of species ' // TRIM(variable) // & ! TRIM(chem_species(lsp)%name) IF (av == 0) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = MERGE( & chem_species(lsp)%conc(k,j,i), & REAL( fill_value, KIND = wp ), & BTEST( wall_flags_0(k,j,i), 0 ) ) ENDDO ENDDO ENDDO ELSE DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = MERGE( & chem_species(lsp)%conc_av(k,j,i), & REAL( fill_value, KIND = wp ), & BTEST( wall_flags_0(k,j,i), 0 ) ) ENDDO ENDDO ENDDO ENDIF grid = 'zu' ENDIF ENDDO RETURN END SUBROUTINE chem_data_output_2d ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine defining 3D output variables for chemical species !------------------------------------------------------------------------------! SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do ) USE indices USE kinds USE surface_mod IMPLICIT NONE CHARACTER (LEN=*) :: variable !< INTEGER(iwp) :: av !< INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0) INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d) LOGICAL :: found !< REAL(wp) :: fill_value !< REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !-- local variables CHARACTER(len=16) :: spec_name INTEGER(iwp) :: i, j, k INTEGER(iwp) :: m, l !< running indices for surfaces INTEGER(iwp) :: lsp !< running index for chem spcs found = .FALSE. IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) ) THEN RETURN ENDIF spec_name = TRIM(variable(4:)) IF ( variable(1:3) == 'em_' ) THEN local_pf = 0.0_wp DO lsp = 1, nvar !!! cssws - nvar species, chem_species - nspec species !!! IF ( TRIM(spec_name) == TRIM(chem_species(lsp)%name) ) THEN ! no average for now DO m = 1, surf_usm_h%ns local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = & local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m) ENDDO DO m = 1, surf_lsm_h%ns local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = & local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m) ENDDO DO l = 0, 3 DO m = 1, surf_usm_v(l)%ns local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = & local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) + surf_usm_v(l)%cssws(lsp,m) ENDDO DO m = 1, surf_lsm_v(l)%ns local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = & local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) + surf_lsm_v(l)%cssws(lsp,m) ENDDO ENDDO found = .TRUE. ENDIF ENDDO ELSE DO lsp=1,nspec IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN ! !-- todo: remove or replace by "CALL message" mechanism (kanani) ! IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM(variable) // & ! TRIM(chem_species(lsp)%name) IF (av == 0) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = MERGE( & chem_species(lsp)%conc(k,j,i), & REAL( fill_value, KIND = wp ), & BTEST( wall_flags_0(k,j,i), 0 ) ) ENDDO ENDDO ENDDO ELSE DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = MERGE( & chem_species(lsp)%conc_av(k,j,i),& REAL( fill_value, KIND = wp ), & BTEST( wall_flags_0(k,j,i), 0 ) ) ENDDO ENDDO ENDDO ENDIF found = .TRUE. ENDIF ENDDO ENDIF RETURN END SUBROUTINE chem_data_output_3d ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine defining mask output variables for chemical species !------------------------------------------------------------------------------! SUBROUTINE chem_data_output_mask( av, variable, found, local_pf ) USE control_parameters USE indices USE kinds USE pegrid, ONLY: myid, threads_per_task USE surface_mod, ONLY: get_topography_top_index_ji IMPLICIT NONE CHARACTER(LEN=5) :: grid !< flag to distinquish between staggered grids CHARACTER (LEN=*):: variable !< INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data LOGICAL :: found !< REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: & local_pf !< !-- local variables. CHARACTER(len=16) :: spec_name INTEGER(iwp) :: lsp INTEGER(iwp) :: i !< grid index along x-direction INTEGER(iwp) :: j !< grid index along y-direction INTEGER(iwp) :: k !< grid index along z-direction INTEGER(iwp) :: topo_top_ind !< k index of highest horizontal surface found = .TRUE. grid = 's' spec_name = TRIM( variable(4:) ) DO lsp=1,nspec IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name) ) THEN ! !-- todo: remove or replace by "CALL message" mechanism (kanani) ! IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM(variable) // & ! TRIM(chem_species(lsp)%name) IF (av == 0) THEN IF ( .NOT. mask_surface(mid) ) THEN DO i = 1, mask_size_l(mid,1) DO j = 1, mask_size_l(mid,2) DO k = 1, mask_size(mid,3) local_pf(i,j,k) = chem_species(lsp)%conc( & mask_k(mid,k), & mask_j(mid,j), & mask_i(mid,i) ) ENDDO ENDDO ENDDO ELSE ! !-- Terrain-following masked output DO i = 1, mask_size_l(mid,1) DO j = 1, mask_size_l(mid,2) ! !-- Get k index of highest horizontal surface topo_top_ind = get_topography_top_index_ji( & mask_j(mid,j), & mask_i(mid,i), & grid ) ! !-- Save output array DO k = 1, mask_size_l(mid,3) local_pf(i,j,k) = chem_species(lsp)%conc( & MIN( topo_top_ind+mask_k(mid,k), & nzt+1 ), & mask_j(mid,j), & mask_i(mid,i) ) ENDDO ENDDO ENDDO ENDIF ELSE IF ( .NOT. mask_surface(mid) ) THEN DO i = 1, mask_size_l(mid,1) DO j = 1, mask_size_l(mid,2) DO k = 1, mask_size_l(mid,3) local_pf(i,j,k) = chem_species(lsp)%conc_av( & mask_k(mid,k), & mask_j(mid,j), & mask_i(mid,i) ) ENDDO ENDDO ENDDO ELSE ! !-- Terrain-following masked output DO i = 1, mask_size_l(mid,1) DO j = 1, mask_size_l(mid,2) ! !-- Get k index of highest horizontal surface topo_top_ind = get_topography_top_index_ji( & mask_j(mid,j), & mask_i(mid,i), & grid ) ! !-- Save output array DO k = 1, mask_size_l(mid,3) local_pf(i,j,k) = chem_species(lsp)%conc_av( & MIN( topo_top_ind+mask_k(mid,k), & nzt+1 ), & mask_j(mid,j), & mask_i(mid,i) ) ENDDO ENDDO ENDDO ENDIF ENDIF found = .FALSE. ENDIF ENDDO RETURN END SUBROUTINE chem_data_output_mask ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine defining appropriate grid for netcdf variables. !> It is called out from subroutine netcdf. !------------------------------------------------------------------------------! SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) IMPLICIT NONE CHARACTER (LEN=*), INTENT(IN) :: var !< LOGICAL, INTENT(OUT) :: found !< CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< found = .TRUE. IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' ) THEN !< always the same grid for chemistry variables grid_x = 'x' grid_y = 'y' grid_z = 'zu' ELSE found = .FALSE. grid_x = 'none' grid_y = 'none' grid_z = 'none' ENDIF END SUBROUTINE chem_define_netcdf_grid ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine defining header output for chemistry model !------------------------------------------------------------------------------! SUBROUTINE chem_header ( io ) IMPLICIT NONE INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file INTEGER(iwp) :: lsp !< running index for chem spcs INTEGER(iwp) :: cs_fixed CHARACTER (LEN=80) :: docsflux_chr CHARACTER (LEN=80) :: docsinit_chr ! !-- Write chemistry model header WRITE( io, 1 ) !-- Gasphase reaction status IF ( chem_gasphase_on ) THEN WRITE( io, 2 ) ELSE WRITE( io, 3 ) ENDIF ! Chemistry time-step WRITE ( io, 4 ) cs_time_step !-- Emission mode info IF ( mode_emis == "DEFAULT" ) THEN WRITE( io, 5 ) ELSEIF ( mode_emis == "PARAMETERIZED" ) THEN WRITE( io, 6 ) ELSEIF ( mode_emis == "PRE-PROCESSED" ) THEN WRITE( io, 7 ) ENDIF !-- Photolysis scheme info IF ( photolysis_scheme == "simple" ) THEN WRITE( io, 8 ) ELSEIF (photolysis_scheme == "conastant" ) THEN WRITE( io, 9 ) ENDIF !-- Emission flux info lsp = 1 docsflux_chr ='Chemical species for surface emission flux: ' DO WHILE ( surface_csflux_name(lsp) /= 'novalue' ) docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ',' IF ( LEN_TRIM( docsflux_chr ) >= 75 ) THEN WRITE ( io, 10 ) docsflux_chr docsflux_chr = ' ' ENDIF lsp = lsp + 1 ENDDO IF ( docsflux_chr /= '' ) THEN WRITE ( io, 10 ) docsflux_chr ENDIF !-- initializatoin of Surface and profile chemical species lsp = 1 docsinit_chr ='Chemical species for initial surface and profile emissions: ' DO WHILE ( cs_name(lsp) /= 'novalue' ) docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ',' IF ( LEN_TRIM( docsinit_chr ) >= 75 ) THEN WRITE ( io, 11 ) docsinit_chr docsinit_chr = ' ' ENDIF lsp = lsp + 1 ENDDO IF ( docsinit_chr /= '' ) THEN WRITE ( io, 11 ) docsinit_chr ENDIF !-- number of variable and fix chemical species and number of reactions cs_fixed = nspec - nvar WRITE ( io, * ) ' --> Chemical species, variable: ', nvar WRITE ( io, * ) ' --> Chemical species, fixed : ', cs_fixed WRITE ( io, * ) ' --> Total number of reactions : ', nreact 1 FORMAT (//' Chemistry model information:'/ & ' ----------------------------'/) 2 FORMAT (' --> Chemical reactions are turned on') 3 FORMAT (' --> Chemical reactions are turned off') 4 FORMAT (' --> Time-step for chemical species: ',F6.2, ' s') 5 FORMAT (' --> Emission mode = DEFAULT ') 6 FORMAT (' --> Emission mode = PARAMETERIZED ') 7 FORMAT (' --> Emission mode = PRE-PROCESSED ') 8 FORMAT (' --> Photolysis scheme used = simple ') 9 FORMAT (' --> Photolysis scheme used = constant ') 10 FORMAT (/' ',A) 11 FORMAT (/' ',A) ! ! END SUBROUTINE chem_header ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine initializating chemistry_model_mod !------------------------------------------------------------------------------! SUBROUTINE chem_init USE control_parameters, & ONLY: message_string, io_blocks, io_group, turbulent_inflow USE arrays_3d, & ONLY: mean_inflow_profiles USE pegrid IMPLICIT none !-- local variables INTEGER :: i,j !< running index for for horiz numerical grid points INTEGER :: lsp !< running index for chem spcs INTEGER :: lpr_lev !< running index for chem spcs profile level ! !-- NOPOINTER version not implemented yet ! #if defined( __nopointer ) ! message_string = 'The chemistry module only runs with POINTER version' ! CALL message( 'chemistry_model_mod', 'CM0001', 1, 2, 0, 6, 0 ) ! #endif ! !-- Allocate memory for chemical species ALLOCATE( chem_species(nspec) ) ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) ALLOCATE( phot_frequen(nphot) ) ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) ) ALLOCATE( bc_cs_t_val(nspec) ) ! !-- Initialize arrays spec_conc_1 (:,:,:,:) = 0.0_wp spec_conc_2 (:,:,:,:) = 0.0_wp spec_conc_3 (:,:,:,:) = 0.0_wp spec_conc_av(:,:,:,:) = 0.0_wp DO lsp = 1, nspec chem_species(lsp)%name = spc_names(lsp) chem_species(lsp)%conc (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1 (:,:,:,lsp) chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2 (:,:,:,lsp) chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_3 (:,:,:,lsp) chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_av(:,:,:,lsp) ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg)) chem_species(lsp)%cssws_av = 0.0_wp ! !-- The following block can be useful when emission module is not applied. & !-- if emission module is applied the following block will be overwritten. ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1)) ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1)) ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1)) ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1)) chem_species(lsp)%flux_s_cs = 0.0_wp chem_species(lsp)%flux_l_cs = 0.0_wp chem_species(lsp)%diss_s_cs = 0.0_wp chem_species(lsp)%diss_l_cs = 0.0_wp ! !-- Allocate memory for initial concentration profiles !-- (concentration values come from namelist) !-- (@todo (FK): Because of this, chem_init is called in palm before !-- check_parameters, since conc_pr_init is used there. !-- We have to find another solution since chem_init should !-- eventually be called from init_3d_model!!) ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) ) chem_species(lsp)%conc_pr_init(:) = 0.0_wp ENDDO ! !-- Initial concentration of profiles is prescribed by parameters cs_profile !-- and cs_heights in the namelist &chemistry_parameters CALL chem_init_profiles ! !-- Initialize model variables IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN !-- First model run of a possible job queue. !-- Initial profiles of the variables must be computed. IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN CALL location_message( 'initializing with 1D chemistry model profiles', .FALSE. ) ! !-- Transfer initial profiles to the arrays of the 3D model DO lsp = 1, nspec DO i = nxlg, nxrg DO j = nysg, nyng DO lpr_lev = 1, nz + 1 chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev) ENDDO ENDDO ENDDO ENDDO ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) & THEN CALL location_message( 'initializing with constant chemistry profiles', .FALSE. ) DO lsp = 1, nspec DO i = nxlg, nxrg DO j = nysg, nyng chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init ENDDO ENDDO ENDDO ENDIF ! !-- If required, change the surface chem spcs at the start of the 3D run IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN DO lsp = 1, nspec chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) + & cs_surface_initial_change(lsp) ENDDO ENDIF ! !-- Initiale old and new time levels. DO lsp = 1, nvar chem_species(lsp)%tconc_m = 0.0_wp chem_species(lsp)%conc_p = chem_species(lsp)%conc ENDDO ENDIF !--- new code add above this line DO lsp = 1, nphot phot_frequen(lsp)%name = phot_names(lsp) ! !-- todo: remove or replace by "CALL message" mechanism (kanani) ! IF( myid == 0 ) THEN ! WRITE(6,'(a,i4,3x,a)') 'Photolysis: ',lsp,trim(phot_names(lsp)) ! ENDIF phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => freq_1(:,:,:,lsp) ENDDO RETURN END SUBROUTINE chem_init ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine defining initial vertical profiles of chemical species (given by !> namelist parameters chem_profiles and chem_heights) --> which should work !> analogue to parameters u_profile, v_profile and uv_heights) !------------------------------------------------------------------------------! SUBROUTINE chem_init_profiles !< SUBROUTINE is called from chem_init in case of !< TRIM( initializing_actions ) /= 'read_restart_data' !< We still need to see what has to be done in case of restart run USE chem_modules IMPLICIT NONE !-- Local variables INTEGER :: lsp !< running index for number of species in derived data type species_def INTEGER :: lsp_usr !< running index for number of species (user defined) in cs_names, cs_profiles etc INTEGER :: lpr_lev !< running index for profile level for each chem spcs. INTEGER :: npr_lev !< the next available profile lev !----------------- !-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles !-- and heights. If parameter "cs_profile" is not prescribed then initial surface values !-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and !-- "cs_heights" are prescribed, their values will!override the constant profile given by !-- "cs_surface". ! IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN lsp_usr = 1 DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' ) !'novalue' is the default DO lsp = 1, nspec ! !-- create initial profile (conc_pr_init) for each chemical species IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) ) THEN ! IF ( cs_profile(lsp_usr,1) == 9999999.9_wp ) THEN !-- set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species DO lpr_lev = 0, nzt+1 chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr) ENDDO ELSE IF ( cs_heights(1,1) /= 0.0_wp ) THEN message_string = 'The surface value of cs_heights must be 0.0' CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 ) ENDIF use_prescribed_profile_data = .TRUE. npr_lev = 1 ! chem_species(lsp)%conc_pr_init(0) = 0.0_wp DO lpr_lev = 1, nz+1 IF ( npr_lev < 100 ) THEN DO WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) ) npr_lev = npr_lev + 1 IF ( npr_lev == 100 ) THEN message_string = 'number of chem spcs exceeding the limit' CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 ) EXIT ENDIF ENDDO ENDIF IF ( npr_lev < 100 .AND. cs_heights(lsp_usr, npr_lev + 1) /= 9999999.9_wp ) THEN chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) + & ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) / & ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) * & ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) ) ELSE chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) ENDIF ENDDO ENDIF !-- If a profile is prescribed explicity using cs_profiles and cs_heights, then !-- chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based !-- on the cs_profiles(lsp_usr,:) and cs_heights(lsp_usr,:). ENDIF ENDDO lsp_usr = lsp_usr + 1 ENDDO ENDIF END SUBROUTINE chem_init_profiles ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine to integrate chemical species in the given chemical mechanism !------------------------------------------------------------------------------! SUBROUTINE chem_integrate_ij (i, j) USE cpulog, & ONLY: cpu_log, log_point USE statistics, & ONLY: weight_pres USE control_parameters, & ONLY: dt_3d, intermediate_timestep_count,simulated_time IMPLICIT none INTEGER,INTENT(IN) :: i,j !-- local variables INTEGER(iwp) :: lsp !< running index for chem spcs. INTEGER(iwp) :: lph !< running index for photolysis frequencies INTEGER :: k,m,istatf INTEGER,DIMENSION(20) :: istatus REAL(kind=wp),DIMENSION(nzb+1:nzt,nspec) :: tmp_conc REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_temp REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_qvap REAL(kind=wp),DIMENSION(nzb+1:nzt,nphot) :: tmp_phot REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_fact REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_fact_i !< conversion factor between !< molecules cm^{-3} and ppm INTEGER,DIMENSION(nzb+1:nzt) :: nacc !< Number of accepted steps INTEGER,DIMENSION(nzb+1:nzt) :: nrej !< Number of rejected steps REAL(wp) :: conv !< conversion factor REAL(wp), PARAMETER :: ppm2fr = 1.0e-6_wp !< Conversion factor ppm to fraction REAL(wp), PARAMETER :: xm_air = 28.96_wp !< Mole mass of dry air REAL(wp), PARAMETER :: xm_h2o = 18.01528_wp !< Mole mass of water vapor REAL(wp), PARAMETER :: pref_i = 1._wp / 100000.0_wp !< inverse reference pressure (1/Pa) REAL(wp), PARAMETER :: t_std = 273.15_wp !< standard pressure (Pa) REAL(wp), PARAMETER :: p_std = 101325.0_wp !< standard pressure (Pa) REAL(wp), PARAMETER :: vmolcm = 22.414e3_wp !< Mole volume (22.414 l) in cm^{-3} REAL(wp), PARAMETER :: xna = 6.022e23_wp !< Avogadro number (molecules/mol) REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local REAL(kind=wp) :: dt_chem CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'start' ) !< set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry IF (chem_gasphase_on) THEN nacc = 0 nrej = 0 tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt) ! ppm to molecules/cm**3 ! tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp * hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp ) conv = ppm2fr * xna / vmolcm tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std) tmp_fact_i = 1.0_wp/tmp_fact IF ( humidity ) THEN IF ( bulk_cloud_model ) THEN tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) * xm_air/xm_h2o * tmp_fact(:) ELSE tmp_qvap(:) = q(nzb+1:nzt,j,i) * xm_air/xm_h2o * tmp_fact(:) ENDIF ELSE tmp_qvap(:) = 0.01 * tmp_fact(:) !< Constant value for q if water vapor is not computed ENDIF DO lsp = 1,nspec tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:) ENDDO DO lph = 1,nphot tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i) ENDDO ! !-- todo: remove (kanani) ! IF(myid == 0 .AND. chem_debug0 ) THEN ! IF (i == 10 .and. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d ! ENDIF !-- Compute length of time step IF ( call_chem_at_all_substeps ) THEN dt_chem = dt_3d * weight_pres(intermediate_timestep_count) ELSE dt_chem = dt_3d ENDIF cs_time_step = dt_chem CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'start' ) IF(maxval(rcntrl) > 0.0) THEN ! Only if rcntrl is set IF( simulated_time <= 2*dt_3d) THEN rcntrl_local = 0 ! !-- todo: remove (kanani) ! WRITE(9,'(a,2f10.3)') 'USE Default rcntrl in the first steps ',simulated_time,dt_3d ELSE rcntrl_local = rcntrl ENDIF ELSE rcntrl_local = 0 END IF CALL chem_gasphase_integrate (dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot, & icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus) CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'stop' ) DO lsp = 1,nspec chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:) ENDDO ENDIF CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'stop' ) RETURN END SUBROUTINE chem_integrate_ij ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine defining parin for &chemistry_parameters for chemistry model !------------------------------------------------------------------------------! SUBROUTINE chem_parin USE chem_modules USE control_parameters USE kinds USE pegrid USE statistics IMPLICIT NONE CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file CHARACTER (LEN=3) :: cs_prefix REAL(wp), DIMENSION(nmaxfixsteps) :: my_steps !< List of fixed timesteps my_step(1) = 0.0 automatic stepping INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: max_pr_cs_tmp !< NAMELIST /chemistry_parameters/ bc_cs_b, & bc_cs_t, & call_chem_at_all_substeps, & chem_debug0, & chem_debug1, & chem_debug2, & chem_gasphase_on, & cs_heights, & cs_name, & cs_profile, & cs_surface, & decycle_chem_lr, & decycle_chem_ns, & decycle_method, & do_depo, & emiss_factor_main, & emiss_factor_side, & icntrl, & main_street_id, & max_street_id, & my_steps, & nest_chemistry, & rcntrl, & side_street_id, & photolysis_scheme, & wall_csflux, & cs_vertical_gradient, & top_csflux, & surface_csflux, & surface_csflux_name, & cs_surface_initial_change, & cs_vertical_gradient_level, & ! namelist parameters for emissions mode_emis, & time_fac_type, & daytype_mdh, & do_emis !-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj) !-- so this way we could prescribe a specific flux value for each species !> chemistry_parameters for initial profiles !> cs_names = 'O3', 'NO2', 'NO', ... to set initial profiles) !> cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3) !> cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.) !> cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, ..... (chem spcs conc at height lvls chem_heights(1,:)) etc. !> If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)" !> then write these cs_surface values to chem_species(lsp)%conc_pr_init(:) ! !-- Read chem namelist INTEGER :: ier CHARACTER(LEN=64) :: text CHARACTER(LEN=8) :: solver_type icntrl = 0 rcntrl = 0.0_wp my_steps = 0.0_wp photolysis_scheme = 'simple' atol = 1.0_wp rtol = 0.01_wp ! !-- Try to find chemistry package REWIND ( 11 ) line = ' ' DO WHILE ( INDEX( line, '&chemistry_parameters' ) == 0 ) READ ( 11, '(A)', END=20 ) line ENDDO BACKSPACE ( 11 ) ! !-- Read chemistry namelist READ ( 11, chemistry_parameters, ERR = 10, END = 20 ) ! !-- Enable chemistry model air_chemistry = .TRUE. GOTO 20 10 BACKSPACE( 11 ) READ( 11 , '(A)') line CALL parin_fail_message( 'chemistry_parameters', line ) 20 CONTINUE ! !-- check for emission mode for chem species IF ( (mode_emis /= 'PARAMETERIZED') .AND. ( mode_emis /= 'DEFAULT') .AND. (mode_emis /= 'PRE-PROCESSED' ) ) THEN message_string = 'Incorrect mode_emiss option select. Please check spelling' CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 ) ENDIF t_steps = my_steps !-- Determine the number of user-defined profiles and append them to the !-- standard data output (data_output_pr) max_pr_cs_tmp = 0 i = 1 DO WHILE ( data_output_pr(i) /= ' ' .AND. i <= 100 ) IF ( TRIM( data_output_pr(i)(1:3)) == 'kc_' ) THEN max_pr_cs_tmp = max_pr_cs_tmp+1 ENDIF i = i +1 ENDDO IF ( max_pr_cs_tmp > 0 ) THEN cs_pr_namelist_found = .TRUE. max_pr_cs = max_pr_cs_tmp ENDIF ! Set Solver Type IF(icntrl(3) == 0) THEN solver_type = 'rodas3' !Default ELSE IF(icntrl(3) == 1) THEN solver_type = 'ros2' ELSE IF(icntrl(3) == 2) THEN solver_type = 'ros3' ELSE IF(icntrl(3) == 3) THEN solver_type = 'ro4' ELSE IF(icntrl(3) == 4) THEN solver_type = 'rodas3' ELSE IF(icntrl(3) == 5) THEN solver_type = 'rodas4' ELSE IF(icntrl(3) == 6) THEN solver_type = 'Rang3' ELSE message_string = 'illegal solver type' CALL message( 'chem_parin', 'PA0506', 1, 2, 0, 6, 0 ) END IF ! !-- todo: remove or replace by "CALL message" mechanism (kanani) ! write(text,*) 'gas_phase chemistry: solver_type = ',trim(solver_type) !kk Has to be changed to right calling sequence !kk CALL location_message( trim(text), .FALSE. ) ! IF(myid == 0) THEN ! write(9,*) ' ' ! write(9,*) 'kpp setup ' ! write(9,*) ' ' ! write(9,*) ' gas_phase chemistry: solver_type = ',trim(solver_type) ! write(9,*) ' ' ! write(9,*) ' Hstart = ',rcntrl(3) ! write(9,*) ' FacMin = ',rcntrl(4) ! write(9,*) ' FacMax = ',rcntrl(5) ! write(9,*) ' ' ! IF(vl_dim > 1) THEN ! write(9,*) ' Vector mode vektor length = ',vl_dim ! ELSE ! write(9,*) ' Scalar mode' ! ENDIF ! write(9,*) ' ' ! END IF RETURN END SUBROUTINE chem_parin ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine calculating prognostic equations for chemical species !> (vector-optimized). !> Routine is called separately for each chemical species over a loop from !> prognostic_equations. !------------------------------------------------------------------------------! SUBROUTINE chem_prognostic_equations ( cs_scalar_p, cs_scalar, tcs_scalar_m, & pr_init_cs, ilsp ) USE advec_s_pw_mod, & ONLY: advec_s_pw USE advec_s_up_mod, & ONLY: advec_s_up USE advec_ws, & ONLY: advec_s_ws USE diffusion_s_mod, & ONLY: diffusion_s USE indices, & ONLY: nxl, nxr, nyn, nys, wall_flags_0 USE pegrid USE surface_mod, & ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & surf_usm_v IMPLICIT NONE INTEGER :: i !< running index INTEGER :: j !< running index INTEGER :: k !< running index INTEGER(iwp),INTENT(IN) :: ilsp !< REAL(wp), DIMENSION(0:nz+1) :: pr_init_cs !< REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar !< REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar_p !< REAL(wp), DIMENSION(:,:,:), POINTER :: tcs_scalar_m !< ! !-- Tendency terms for chemical species tend = 0.0_wp ! !-- Advection terms IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( cs_scalar, 'kc' ) ELSE CALL advec_s_pw( cs_scalar ) ENDIF ELSE CALL advec_s_up( cs_scalar ) ENDIF ! !-- Diffusion terms (the last three arguments are zero) CALL diffusion_s( cs_scalar, & surf_def_h(0)%cssws(ilsp,:), & surf_def_h(1)%cssws(ilsp,:), & surf_def_h(2)%cssws(ilsp,:), & surf_lsm_h%cssws(ilsp,:), & surf_usm_h%cssws(ilsp,:), & surf_def_v(0)%cssws(ilsp,:), & surf_def_v(1)%cssws(ilsp,:), & surf_def_v(2)%cssws(ilsp,:), & surf_def_v(3)%cssws(ilsp,:), & surf_lsm_v(0)%cssws(ilsp,:), & surf_lsm_v(1)%cssws(ilsp,:), & surf_lsm_v(2)%cssws(ilsp,:), & surf_lsm_v(3)%cssws(ilsp,:), & surf_usm_v(0)%cssws(ilsp,:), & surf_usm_v(1)%cssws(ilsp,:), & surf_usm_v(2)%cssws(ilsp,:), & surf_usm_v(3)%cssws(ilsp,:) ) ! !-- Prognostic equation for chemical species DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt cs_scalar_p(k,j,i) = cs_scalar(k,j,i) & + ( dt_3d * & ( tsc(2) * tend(k,j,i) & + tsc(3) * tcs_scalar_m(k,j,i) & ) & - tsc(5) * rdf_sc(k) & * ( cs_scalar(k,j,i) - pr_init_cs(k) ) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( cs_scalar_p(k,j,i) < 0.0_wp ) cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i) ENDDO ENDDO ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tcs_scalar_m(k,j,i) = tend(k,j,i) ENDDO ENDDO ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tcs_scalar_m(k,j,i) = - 9.5625_wp * tend(k,j,i) & + 5.3125_wp * tcs_scalar_m(k,j,i) ENDDO ENDDO ENDDO ENDIF ENDIF END SUBROUTINE chem_prognostic_equations !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine calculating prognostic equations for chemical species !> (cache-optimized). !> Routine is called separately for each chemical species over a loop from !> prognostic_equations. !------------------------------------------------------------------------------! SUBROUTINE chem_prognostic_equations_ij ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs, & i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs, & flux_l_cs, diss_l_cs ) USE pegrid USE advec_ws, ONLY: advec_s_ws USE advec_s_pw_mod, ONLY: advec_s_pw USE advec_s_up_mod, ONLY: advec_s_up USE diffusion_s_mod, ONLY: diffusion_s USE indices, ONLY: wall_flags_0 USE surface_mod, ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & surf_usm_v IMPLICIT NONE REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar_p, cs_scalar, tcs_scalar_m INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn, ilsp REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: flux_s_cs !< REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: diss_s_cs !< REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: flux_l_cs !< REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: diss_l_cs !< REAL(wp), DIMENSION(0:nz+1) :: pr_init_cs !< !-- local variables INTEGER :: k ! !-- Tendency-terms for chem spcs. tend(:,j,i) = 0.0_wp ! !-- Advection terms IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( i, j, cs_scalar, 'kc', flux_s_cs, diss_s_cs, & flux_l_cs, diss_l_cs, i_omp_start, tn ) ELSE CALL advec_s_pw( i, j, cs_scalar ) ENDIF ELSE CALL advec_s_up( i, j, cs_scalar ) ENDIF ! !-- Diffusion terms (the last three arguments are zero) CALL diffusion_s( i, j, cs_scalar, & surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:), & surf_def_h(2)%cssws(ilsp,:), & surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:), & surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:), & surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:), & surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:), & surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:), & surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:), & surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) ) ! !-- Prognostic equation for chem spcs DO k = nzb+1, nzt cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d * & ( tsc(2) * tend(k,j,i) + & tsc(3) * tcs_scalar_m(k,j,i) ) & - tsc(5) * rdf_sc(k) & * ( cs_scalar(k,j,i) - pr_init_cs(k) ) & ) & * MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) & ) IF ( cs_scalar_p(k,j,i) < 0.0_wp ) cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i) !FKS6 ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = nzb+1, nzt tcs_scalar_m(k,j,i) = tend(k,j,i) ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = nzb+1, nzt tcs_scalar_m(k,j,i) = -9.5625_wp * tend(k,j,i) + & 5.3125_wp * tcs_scalar_m(k,j,i) ENDDO ENDIF ENDIF END SUBROUTINE chem_prognostic_equations_ij ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine to read restart data of chemical species !------------------------------------------------------------------------------! SUBROUTINE chem_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, & nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc, & nys_on_file, tmp_3d, found ) USE control_parameters USE indices USE pegrid IMPLICIT NONE CHARACTER (LEN=20) :: spc_name_av !< INTEGER(iwp) :: i, lsp !< INTEGER(iwp) :: k !< INTEGER(iwp) :: nxlc !< INTEGER(iwp) :: nxlf !< INTEGER(iwp) :: nxl_on_file !< INTEGER(iwp) :: nxrc !< INTEGER(iwp) :: nxrf !< INTEGER(iwp) :: nxr_on_file !< INTEGER(iwp) :: nync !< INTEGER(iwp) :: nynf !< INTEGER(iwp) :: nyn_on_file !< INTEGER(iwp) :: nysc !< INTEGER(iwp) :: nysf !< INTEGER(iwp) :: nys_on_file !< LOGICAL, INTENT(OUT) :: found REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !< 3D array to temp store data found = .FALSE. IF ( ALLOCATED(chem_species) ) THEN DO lsp = 1, nspec !< for time-averaged chemical conc. spc_name_av = TRIM(chem_species(lsp)%name)//'_av' IF (restart_string(1:length) == TRIM(chem_species(lsp)%name) ) & THEN !< read data into tmp_3d IF ( k == 1 ) READ ( 13 ) tmp_3d !< fill ..%conc in the restart run chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp, & nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) found = .TRUE. ELSEIF (restart_string(1:length) == spc_name_av ) THEN IF ( k == 1 ) READ ( 13 ) tmp_3d chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp, & nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) found = .TRUE. ENDIF ENDDO ENDIF END SUBROUTINE chem_rrd_local ! !-------------------------------------------------------------------------------! !> Description: !> Calculation of horizontally averaged profiles !> This routine is called for every statistic region (sr) defined by the user, !> but at least for the region "total domain" (sr=0). !> quantities. !------------------------------------------------------------------------------! SUBROUTINE chem_statistics( mode, sr, tn ) USE arrays_3d USE indices USE kinds USE pegrid USE statistics USE user IMPLICIT NONE CHARACTER (LEN=*) :: mode !< INTEGER(iwp) :: i !< running index on x-axis INTEGER(iwp) :: j !< running index on y-axis INTEGER(iwp) :: k !< vertical index counter INTEGER(iwp) :: sr !< statistical region INTEGER(iwp) :: tn !< thread number INTEGER(iwp) :: n !< INTEGER(iwp) :: m !< INTEGER(iwp) :: lpr !< running index chem spcs ! REAL(wp), & ! DIMENSION(dots_num_palm+1:dots_max) :: & ! ts_value_l !< IF ( mode == 'profiles' ) THEN ! !-- Sample on how to calculate horizontally averaged profiles of user- !-- defined quantities. Each quantity is identified by the index !-- "pr_palm+#" where "#" is an integer starting from 1. These !-- user-profile-numbers must also be assigned to the respective strings !-- given by data_output_pr_cs in routine user_check_data_output_pr. !-- hom(:,:,:,:) = dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g. ! w*pt*), dim-4 = statistical region. !$OMP DO DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 DO lpr = 1, cs_pr_count sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) + & chem_species(cs_pr_index(lpr))%conc(k,j,i) * & rmask(j,i,sr) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 22 ) ) ENDDO ENDDO ENDDO ENDDO ELSEIF ( mode == 'time_series' ) THEN CALL location_message( 'Time series not calculated for chemistry', .TRUE. ) ENDIF END SUBROUTINE chem_statistics ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine for swapping of timelevels for chemical species !> called out from subroutine swap_timelevel !------------------------------------------------------------------------------! SUBROUTINE chem_swap_timelevel (level) IMPLICIT none INTEGER,INTENT(IN) :: level !-- local variables INTEGER :: lsp IF ( level == 0 ) THEN DO lsp=1, nvar chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp) chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp) ENDDO ELSE DO lsp=1, nvar chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp) chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp) ENDDO ENDIF RETURN END SUBROUTINE chem_swap_timelevel ! !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine to write restart data for chemistry model !------------------------------------------------------------------------------! SUBROUTINE chem_wrd_local IMPLICIT NONE INTEGER(iwp) :: lsp !< DO lsp = 1, nspec CALL wrd_write_string( TRIM( chem_species(lsp)%name )) WRITE ( 14 ) chem_species(lsp)%conc CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' ) WRITE ( 14 ) chem_species(lsp)%conc_av ENDDO END SUBROUTINE chem_wrd_local !-------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine to calculate the deposition of gases and PMs. For now deposition !> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT !> considered. The deposition of particles is derived following Zhang et al., !> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010). !> !> @TODO: Consider deposition on vertical surfaces !> @TODO: Consider overlaying horizontal surfaces !> @TODO: Check error messages !-------------------------------------------------------------------------------! SUBROUTINE chem_depo (i,j) USE control_parameters, & ONLY: dt_3d, intermediate_timestep_count, latitude USE arrays_3d, & ONLY: dzw, rho_air_zw USE date_and_time_mod, & ONLY: day_of_year USE surface_mod, & ONLY: surf_lsm_h, surf_usm_h, surf_type, ind_veg_wall, & ind_pav_green, ind_wat_win USE radiation_model_mod, & ONLY: zenith IMPLICIT NONE INTEGER,INTENT(IN) :: i,j INTEGER(iwp) :: k ! matching k to surface m at i,j INTEGER(iwp) :: lsp ! running index for chem spcs. INTEGER(iwp) :: lu ! running index for landuse classes INTEGER(iwp) :: lu_palm ! index of PALM LSM vegetation_type at current surface element INTEGER(iwp) :: lup_palm ! index of PALM LSM pavement_type at current surface element INTEGER(iwp) :: luw_palm ! index of PALM LSM water_type at current surface element INTEGER(iwp) :: luu_palm ! index of PALM USM walls/roofs at current surface element INTEGER(iwp) :: lug_palm ! index of PALM USM green walls/roofs at current surface element INTEGER(iwp) :: lud_palm ! index of PALM USM windows at current surface element INTEGER(iwp) :: lu_dep ! matching DEPAC LU to lu_palm INTEGER(iwp) :: lup_dep ! matching DEPAC LU to lup_palm INTEGER(iwp) :: luw_dep ! matching DEPAC LU to luw_palm INTEGER(iwp) :: luu_dep ! matching DEPAC LU to luu_palm INTEGER(iwp) :: lug_dep ! matching DEPAC LU to lug_palm INTEGER(iwp) :: lud_dep ! matching DEPAC LU to lud_palm INTEGER(iwp) :: m ! index for horizontal surfaces INTEGER(iwp) :: pspec ! running index INTEGER(iwp) :: i_pspec ! index for matching depac gas component ! Vegetation !Match to DEPAC INTEGER(iwp) :: ind_lu_user = 0 ! --> ERROR INTEGER(iwp) :: ind_lu_b_soil = 1 ! --> ilu_desert INTEGER(iwp) :: ind_lu_mixed_crops = 2 ! --> ilu_arable INTEGER(iwp) :: ind_lu_s_grass = 3 ! --> ilu_grass INTEGER(iwp) :: ind_lu_ev_needle_trees = 4 ! --> ilu_coniferous_forest INTEGER(iwp) :: ind_lu_de_needle_trees = 5 ! --> ilu_coniferous_forest INTEGER(iwp) :: ind_lu_ev_broad_trees = 6 ! --> ilu_tropical_forest INTEGER(iwp) :: ind_lu_de_broad_trees = 7 ! --> ilu_deciduous_forest INTEGER(iwp) :: ind_lu_t_grass = 8 ! --> ilu_grass INTEGER(iwp) :: ind_lu_desert = 9 ! --> ilu_desert INTEGER(iwp) :: ind_lu_tundra = 10 ! --> ilu_other INTEGER(iwp) :: ind_lu_irr_crops = 11 ! --> ilu_arable INTEGER(iwp) :: ind_lu_semidesert = 12 ! --> ilu_other INTEGER(iwp) :: ind_lu_ice = 13 ! --> ilu_ice INTEGER(iwp) :: ind_lu_marsh = 14 ! --> ilu_other INTEGER(iwp) :: ind_lu_ev_shrubs = 15 ! --> ilu_mediterrean_scrub INTEGER(iwp) :: ind_lu_de_shrubs = 16 ! --> ilu_mediterrean_scrub INTEGER(iwp) :: ind_lu_mixed_forest = 17 ! --> ilu_coniferous_forest (ave(decid+conif)) INTEGER(iwp) :: ind_lu_intrup_forest = 18 ! --> ilu_other (ave(other+decid)) ! Water INTEGER(iwp) :: ind_luw_user = 0 ! --> ERROR INTEGER(iwp) :: ind_luw_lake = 1 ! --> ilu_water_inland INTEGER(iwp) :: ind_luw_river = 2 ! --> ilu_water_inland INTEGER(iwp) :: ind_luw_ocean = 3 ! --> ilu_water_sea INTEGER(iwp) :: ind_luw_pond = 4 ! --> ilu_water_inland INTEGER(iwp) :: ind_luw_fountain = 5 ! --> ilu_water_inland ! Pavement INTEGER(iwp) :: ind_lup_user = 0 ! --> ERROR INTEGER(iwp) :: ind_lup_asph_conc = 1 ! --> ilu_desert INTEGER(iwp) :: ind_lup_asph = 2 ! --> ilu_desert INTEGER(iwp) :: ind_lup_conc = 3 ! --> ilu_desert INTEGER(iwp) :: ind_lup_sett = 4 ! --> ilu_desert INTEGER(iwp) :: ind_lup_pav_stones = 5 ! --> ilu_desert INTEGER(iwp) :: ind_lup_cobblest = 6 ! --> ilu_desert INTEGER(iwp) :: ind_lup_metal = 7 ! --> ilu_desert INTEGER(iwp) :: ind_lup_wood = 8 ! --> ilu_desert INTEGER(iwp) :: ind_lup_gravel = 9 ! --> ilu_desert INTEGER(iwp) :: ind_lup_f_gravel = 10 ! --> ilu_desert INTEGER(iwp) :: ind_lup_pebblest = 11 ! --> ilu_desert INTEGER(iwp) :: ind_lup_woodchips = 12 ! --> ilu_desert INTEGER(iwp) :: ind_lup_tartan = 13 ! --> ilu_desert INTEGER(iwp) :: ind_lup_art_turf = 14 ! --> ilu_desert INTEGER(iwp) :: ind_lup_clay = 15 ! --> ilu_desert !-- Particle parameters according to the respective aerosol classes (PM25, PM10) INTEGER(iwp) :: ind_p_size = 1 !< index for partsize in particle_pars INTEGER(iwp) :: ind_p_dens = 2 !< index for rhopart in particle_pars INTEGER(iwp) :: ind_p_slip = 3 !< index for slipcor in particle_pars INTEGER(iwp) :: part_type INTEGER(iwp) :: nwet ! wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow REAL(kind=wp) :: dt_chem REAL(kind=wp) :: dh REAL(kind=wp) :: inv_dh REAL(kind=wp) :: dt_dh REAL(kind=wp) :: dens ! Density at lowest layer at i,j REAL(kind=wp) :: r_aero_lu ! aerodynamic resistance (s/m) at current surface element REAL(kind=wp) :: ustar_lu ! ustar at current surface element REAL(kind=wp) :: z0h_lu ! roughness length for heat at current surface element REAL(kind=wp) :: glrad ! rad_sw_in at current surface element REAL(kind=wp) :: ppm_to_ugm3 ! conversion factor REAL(kind=wp) :: relh ! relative humidity at current surface element REAL(kind=wp) :: lai ! leaf area index at current surface element REAL(kind=wp) :: sai ! surface area index at current surface element assumed to be lai + 1 REAL(kind=wp) :: slinnfac REAL(kind=wp) :: visc ! Viscosity REAL(kind=wp) :: vs ! Sedimentation velocity REAL(kind=wp) :: vd_lu ! deposition velocity (m/s) REAL(kind=wp) :: Rs ! Sedimentaion resistance (s/m) REAL(kind=wp) :: Rb ! quasi-laminar boundary layer resistance (s/m) REAL(kind=wp) :: Rc_tot ! total canopy resistance Rc (s/m) REAL(kind=wp) :: catm ! gasconc. at i,j,k=1 in ug/m3 REAL(kind=wp) :: diffc ! Diffusivity REAL(kind=wp),DIMENSION(nspec) :: bud_now_lu ! budget for LSM vegetation type at current surface element REAL(kind=wp),DIMENSION(nspec) :: bud_now_lup ! budget for LSM pavement type at current surface element REAL(kind=wp),DIMENSION(nspec) :: bud_now_luw ! budget for LSM water type at current surface element REAL(kind=wp),DIMENSION(nspec) :: bud_now_luu ! budget for USM walls/roofs at current surface element REAL(kind=wp),DIMENSION(nspec) :: bud_now_lug ! budget for USM green surfaces at current surface element REAL(kind=wp),DIMENSION(nspec) :: bud_now_lud ! budget for USM windows at current surface element REAL(kind=wp),DIMENSION(nspec) :: bud_now ! overall budget at current surface element REAL(kind=wp),DIMENSION(nspec) :: cc ! concentration i,j,k REAL(kind=wp),DIMENSION(nspec) :: ccomp_tot ! total compensation point (ug/m3) (= 0 for species that don't have a compensation) ! For now kept to zero for all species! REAL(kind=wp) :: ttemp ! temperatur at i,j,k REAL(kind=wp) :: ts ! surface temperatur in degrees celsius REAL(kind=wp) :: qv_tmp ! surface mixing ratio at current surface element !-- Particle parameters (PM10 (1), PM25 (2)) !-- partsize (diameter in m) , rhopart (density in kg/m3), slipcor (slip correction factor DIMENSIONless, Seinfeld and Pandis 2006, Table 9.3) REAL(wp), DIMENSION(1:3,1:2), PARAMETER :: particle_pars = RESHAPE( (/ & 8.0e-6_wp, 1.14e3_wp, 1.016_wp, & ! 1 0.7e-6_wp, 1.14e3_wp, 1.082_wp & ! 2 /), (/ 3, 2 /) ) LOGICAL :: match_lsm ! flag indicating natural-type surface LOGICAL :: match_usm ! flag indicating urban-type surface !-- List of names of possible tracers CHARACTER(len=*), PARAMETER :: pspecnames(69) = (/ & 'NO2 ', & ! NO2 'NO ', & ! NO 'O3 ', & ! O3 'CO ', & ! CO 'form ', & ! FORM 'ald ', & ! ALD 'pan ', & ! PAN 'mgly ', & ! MGLY 'par ', & ! PAR 'ole ', & ! OLE 'eth ', & ! ETH 'tol ', & ! TOL 'cres ', & ! CRES 'xyl ', & ! XYL 'SO4a_f ', & ! SO4a_f 'SO2 ', & ! SO2 'HNO2 ', & ! HNO2 'CH4 ', & ! CH4 'NH3 ', & ! NH3 'NO3 ', & ! NO3 'OH ', & ! OH 'HO2 ', & ! HO2 'N2O5 ', & ! N2O5 'SO4a_c ', & ! SO4a_c 'NH4a_f ', & ! NH4a_f 'NO3a_f ', & ! NO3a_f 'NO3a_c ', & ! NO3a_c 'C2O3 ', & ! C2O3 'XO2 ', & ! XO2 'XO2N ', & ! XO2N 'cro ', & ! CRO 'HNO3 ', & ! HNO3 'H2O2 ', & ! H2O2 'iso ', & ! ISO 'ispd ', & ! ISPD 'to2 ', & ! TO2 'open ', & ! OPEN 'terp ', & ! TERP 'ec_f ', & ! EC_f 'ec_c ', & ! EC_c 'pom_f ', & ! POM_f 'pom_c ', & ! POM_c 'ppm_f ', & ! PPM_f 'ppm_c ', & ! PPM_c 'na_ff ', & ! Na_ff 'na_f ', & ! Na_f 'na_c ', & ! Na_c 'na_cc ', & ! Na_cc 'na_ccc ', & ! Na_ccc 'dust_ff ', & ! dust_ff 'dust_f ', & ! dust_f 'dust_c ', & ! dust_c 'dust_cc ', & ! dust_cc 'dust_ccc ', & ! dust_ccc 'tpm10 ', & ! tpm10 'tpm25 ', & ! tpm25 'tss ', & ! tss 'tdust ', & ! tdust 'tc ', & ! tc 'tcg ', & ! tcg 'tsoa ', & ! tsoa 'tnmvoc ', & ! tnmvoc 'SOxa ', & ! SOxa 'NOya ', & ! NOya 'NHxa ', & ! NHxa 'NO2_obs ', & ! NO2_obs 'tpm10_biascorr', & ! tpm10_biascorr 'tpm25_biascorr', & ! tpm25_biascorr 'O3_biascorr ' /) ! o3_biascorr ! tracer mole mass: REAL, PARAMETER :: specmolm(69) = (/ & xm_O * 2 + xm_N, & ! NO2 xm_O + xm_N, & ! NO xm_O * 3, & ! O3 xm_C + xm_O, & ! CO xm_H * 2 + xm_C + xm_O, & ! FORM xm_H * 3 + xm_C * 2 + xm_O, & ! ALD xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, & ! PAN xm_H * 4 + xm_C * 3 + xm_O * 2, & ! MGLY xm_H * 3 + xm_C, & ! PAR xm_H * 3 + xm_C * 2, & ! OLE xm_H * 4 + xm_C * 2, & ! ETH xm_H * 8 + xm_C * 7, & ! TOL xm_H * 8 + xm_C * 7 + xm_O, & ! CRES xm_H * 10 + xm_C * 8, & ! XYL xm_S + xm_O * 4, & ! SO4a_f xm_S + xm_O * 2, & ! SO2 xm_H + xm_O * 2 + xm_N, & ! HNO2 xm_H * 4 + xm_C, & ! CH4 xm_H * 3 + xm_N, & ! NH3 xm_O * 3 + xm_N, & ! NO3 xm_H + xm_O, & ! OH xm_H + xm_O * 2, & ! HO2 xm_O * 5 + xm_N * 2, & ! N2O5 xm_S + xm_O * 4, & ! SO4a_c xm_H * 4 + xm_N, & ! NH4a_f xm_O * 3 + xm_N, & ! NO3a_f xm_O * 3 + xm_N, & ! NO3a_c xm_C * 2 + xm_O * 3, & ! C2O3 xm_dummy, & ! XO2 xm_dummy, & ! XO2N xm_dummy, & ! CRO xm_H + xm_O * 3 + xm_N, & ! HNO3 xm_H * 2 + xm_O * 2, & ! H2O2 xm_H * 8 + xm_C * 5, & ! ISO xm_dummy, & ! ISPD xm_dummy, & ! TO2 xm_dummy, & ! OPEN xm_H * 16 + xm_C * 10, & ! TERP xm_dummy, & ! EC_f xm_dummy, & ! EC_c xm_dummy, & ! POM_f xm_dummy, & ! POM_c xm_dummy, & ! PPM_f xm_dummy, & ! PPM_c xm_Na, & ! Na_ff xm_Na, & ! Na_f xm_Na, & ! Na_c xm_Na, & ! Na_cc xm_Na, & ! Na_ccc xm_dummy, & ! dust_ff xm_dummy, & ! dust_f xm_dummy, & ! dust_c xm_dummy, & ! dust_cc xm_dummy, & ! dust_ccc xm_dummy, & ! tpm10 xm_dummy, & ! tpm25 xm_dummy, & ! tss xm_dummy, & ! tdust xm_dummy, & ! tc xm_dummy, & ! tcg xm_dummy, & ! tsoa xm_dummy, & ! tnmvoc xm_dummy, & ! SOxa xm_dummy, & ! NOya xm_dummy, & ! NHxa xm_O * 2 + xm_N, & ! NO2_obs xm_dummy, & ! tpm10_biascorr xm_dummy, & ! tpm25_biascorr xm_O * 3 /) ! o3_biascorr ! Initialize surface element m m = 0 k = 0 ! LSM or USM surface present at i,j: ! Default surfaces are NOT considered for deposition match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i) match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i) !!!!!!! ! LSM ! !!!!!!! IF ( match_lsm ) THEN ! Get surface element information at i,j: m = surf_lsm_h%start_index(j,i) ! Get corresponding k k = surf_lsm_h%k(m) ! Get needed variables for surface element m ustar_lu = surf_lsm_h%us(m) z0h_lu = surf_lsm_h%z0h(m) r_aero_lu = surf_lsm_h%r_a(m) glrad = surf_lsm_h%rad_sw_in(m) lai = surf_lsm_h%lai(m) sai = lai + 1 ! For small grid spacing neglect R_a IF (dzw(k) <= 1.0) THEN r_aero_lu = 0.0_wp ENDIF !Initialize lu's lu_palm = 0 lu_dep = 0 lup_palm = 0 lup_dep = 0 luw_palm = 0 luw_dep = 0 !Initialize budgets bud_now_lu = 0.0_wp bud_now_lup = 0.0_wp bud_now_luw = 0.0_wp ! Get land use for i,j and assign to DEPAC lu IF (surf_lsm_h%frac(ind_veg_wall,m) > 0) THEN lu_palm = surf_lsm_h%vegetation_type(m) IF (lu_palm == ind_lu_user) THEN message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation' CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 ) ELSEIF (lu_palm == ind_lu_b_soil) THEN lu_dep = 9 ELSEIF (lu_palm == ind_lu_mixed_crops) THEN lu_dep = 2 ELSEIF (lu_palm == ind_lu_s_grass) THEN lu_dep = 1 ELSEIF (lu_palm == ind_lu_ev_needle_trees) THEN lu_dep = 4 ELSEIF (lu_palm == ind_lu_de_needle_trees) THEN lu_dep = 4 ELSEIF (lu_palm == ind_lu_ev_broad_trees) THEN lu_dep = 12 ELSEIF (lu_palm == ind_lu_de_broad_trees) THEN lu_dep = 5 ELSEIF (lu_palm == ind_lu_t_grass) THEN lu_dep = 1 ELSEIF (lu_palm == ind_lu_desert) THEN lu_dep = 9 ELSEIF (lu_palm == ind_lu_tundra ) THEN lu_dep = 8 ELSEIF (lu_palm == ind_lu_irr_crops) THEN lu_dep = 2 ELSEIF (lu_palm == ind_lu_semidesert) THEN lu_dep = 8 ELSEIF (lu_palm == ind_lu_ice) THEN lu_dep = 10 ELSEIF (lu_palm == ind_lu_marsh) THEN lu_dep = 8 ELSEIF (lu_palm == ind_lu_ev_shrubs) THEN lu_dep = 14 ELSEIF (lu_palm == ind_lu_de_shrubs ) THEN lu_dep = 14 ELSEIF (lu_palm == ind_lu_mixed_forest) THEN lu_dep = 4 ELSEIF (lu_palm == ind_lu_intrup_forest) THEN lu_dep = 8 ENDIF ENDIF IF (surf_lsm_h%frac(ind_pav_green,m) > 0) THEN lup_palm = surf_lsm_h%pavement_type(m) IF (lup_palm == ind_lup_user) THEN message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation' CALL message( 'chem_depo', 'CM0452', 1, 2, 0, 6, 0 ) ELSEIF (lup_palm == ind_lup_asph_conc) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_asph) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_conc) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_sett) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_pav_stones) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_cobblest) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_metal) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_wood) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_gravel) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_f_gravel) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_pebblest) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_woodchips) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_tartan) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_art_turf) THEN lup_dep = 9 ELSEIF (lup_palm == ind_lup_clay) THEN lup_dep = 9 ENDIF ENDIF IF (surf_lsm_h%frac(ind_wat_win,m) > 0) THEN luw_palm = surf_lsm_h%water_type(m) IF (luw_palm == ind_luw_user) THEN message_string = 'No water type defined. Please define water type to enable deposition calculation' CALL message( 'chem_depo', 'CM0453', 1, 2, 0, 6, 0 ) ELSEIF (luw_palm == ind_luw_lake) THEN luw_dep = 13 ELSEIF (luw_palm == ind_luw_river) THEN luw_dep = 13 ELSEIF (luw_palm == ind_luw_ocean) THEN luw_dep = 6 ELSEIF (luw_palm == ind_luw_pond) THEN luw_dep = 13 ELSEIF (luw_palm == ind_luw_fountain) THEN luw_dep = 13 ENDIF ENDIF ! Set wetness indicator to dry or wet for lsm vegetation or pavement IF (surf_lsm_h%c_liq(m) > 0) THEN nwet = 1 ELSE nwet = 0 ENDIF ! Compute length of time step IF ( call_chem_at_all_substeps ) THEN dt_chem = dt_3d * weight_pres(intermediate_timestep_count) ELSE dt_chem = dt_3d ENDIF dh = dzw(k) inv_dh = 1.0_wp / dh dt_dh = dt_chem/dh ! Concentration at i,j,k DO lsp = 1, nspec cc(lsp) = chem_species(lsp)%conc(k,j,i) ENDDO ! Temperature column at i,j ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp ! Surface temperature in degrees Celcius: ts = ttemp - 273.15 ! in degrees celcius ! Viscosity of air in lowest layer visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0) ! Air density in lowest layer dens = rho_air_zw(k) ! Calculate relative humidity from specific humidity for DEPAC qv_tmp = max(q(k,j,i),0.0_wp) relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(k) ) ! Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget ! for each surface fraction. Then derive overall budget taking into account the surface fractions. ! Vegetation IF (surf_lsm_h%frac(ind_veg_wall,m) > 0) THEN slinnfac = 1.0_wp ! Get vd DO lsp = 1, nvar !Initialize vs = 0.0_wp vd_lu = 0.0_wp Rs = 0.0_wp Rb = 0.0_wp Rc_tot = 0.0_wp IF (spc_names(lsp) == 'PM10') THEN part_type = 1 ! sedimentation velocity in lowest layer vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & lu_dep, & r_aero_lu, ustar_lu) bud_now_lu(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSEIF (spc_names(lsp) == 'PM25') THEN part_type = 2 ! sedimentation velocity in lowest layer: vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & lu_dep , & r_aero_lu, ustar_lu) bud_now_lu(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSE ! GASES ! Read spc_name of current species for gas parameter IF (ANY(pspecnames(:) == spc_names(lsp))) THEN i_pspec = 0 DO pspec = 1, 69 IF (pspecnames(pspec) == spc_names(lsp)) THEN i_pspec = pspec END IF ENDDO ELSE ! Species for now not deposited CYCLE ENDIF ! Factor used for conversion from ppb to ug/m3 : ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) ! c 1e-9 xm_tracer 1e9 / xm_air dens ! thus: ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 ! Use density at lowest layer: ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 ! atmospheric concentration in DEPAC is requested in ug/m3: ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 ! Diffusivity for DEPAC relevant gases ! Use default value diffc = 0.11e-4 ! overwrite with known coefficients of diffusivity from Massman (1998) IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 ! Get quasi-laminar boundary layer resistance Rb: CALL get_rb_cell( (lu_dep == ilu_water_sea) .or. (lu_dep == ilu_water_inland), & z0h_lu, ustar_lu, diffc, & Rb) ! Get Rc_tot CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & r_aero_lu , Rb) ! Calculate budget IF (Rc_tot .le. 0.0) THEN bud_now_lu(lsp) = 0.0_wp ELSE ! Compute exchange velocity for current lu: vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) bud_now_lu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ENDIF ENDIF ENDDO ENDIF ! Pavement IF (surf_lsm_h%frac(ind_pav_green,m) > 0) THEN ! No vegetation on pavements: lai = 0.0_wp sai = 0.0_wp slinnfac = 1.0_wp ! Get vd DO lsp = 1, nvar !Initialize vs = 0.0_wp vd_lu = 0.0_wp Rs = 0.0_wp Rb = 0.0_wp Rc_tot = 0.0_wp IF (spc_names(lsp) == 'PM10') THEN part_type = 1 ! sedimentation velocity in lowest layer: vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & lup_dep, & r_aero_lu, ustar_lu) bud_now_lup(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSEIF (spc_names(lsp) == 'PM25') THEN part_type = 2 ! sedimentation velocity in lowest layer: vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & lup_dep, & r_aero_lu, ustar_lu) bud_now_lup(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSE ! GASES ! Read spc_name of current species for gas parameter IF (ANY(pspecnames(:) == spc_names(lsp))) THEN i_pspec = 0 DO pspec = 1, 69 IF (pspecnames(pspec) == spc_names(lsp)) THEN i_pspec = pspec END IF ENDDO ELSE ! Species for now not deposited CYCLE ENDIF ! Factor used for conversion from ppb to ug/m3 : ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) ! c 1e-9 xm_tracer 1e9 / xm_air dens ! thus: ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 ! Use density at lowest layer: ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 ! atmospheric concentration in DEPAC is requested in ug/m3: ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 ! Diffusivity for DEPAC relevant gases ! Use default value diffc = 0.11e-4 ! overwrite with known coefficients of diffusivity from Massman (1998) IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 ! Get quasi-laminar boundary layer resistance Rb: CALL get_rb_cell( (lup_dep == ilu_water_sea) .or. (lup_dep == ilu_water_inland), & z0h_lu, ustar_lu, diffc, & Rb) ! Get Rc_tot CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & r_aero_lu , Rb) ! Calculate budget IF (Rc_tot .le. 0.0) THEN bud_now_lup(lsp) = 0.0_wp ELSE ! Compute exchange velocity for current lu: vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ENDIF ENDIF ENDDO ENDIF ! Water IF (surf_lsm_h%frac(ind_wat_win,m) > 0) THEN ! No vegetation on water: lai = 0.0_wp sai = 0.0_wp slinnfac = 1.0_wp ! Get vd DO lsp = 1, nvar !Initialize vs = 0.0_wp vd_lu = 0.0_wp Rs = 0.0_wp Rb = 0.0_wp Rc_tot = 0.0_wp IF (spc_names(lsp) == 'PM10') THEN part_type = 1 ! sedimentation velocity in lowest layer: vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & luw_dep, & r_aero_lu, ustar_lu) bud_now_luw(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSEIF (spc_names(lsp) == 'PM25') THEN part_type = 2 ! sedimentation velocity in lowest layer: vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & luw_dep, & r_aero_lu, ustar_lu) bud_now_luw(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSE ! GASES ! Read spc_name of current species for gas PARAMETER IF (ANY(pspecnames(:) == spc_names(lsp))) THEN i_pspec = 0 DO pspec = 1, 69 IF (pspecnames(pspec) == spc_names(lsp)) THEN i_pspec = pspec END IF ENDDO ELSE ! Species for now not deposited CYCLE ENDIF ! Factor used for conversion from ppb to ug/m3 : ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) ! c 1e-9 xm_tracer 1e9 / xm_air dens ! thus: ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 ! Use density at lowest layer: ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 ! atmospheric concentration in DEPAC is requested in ug/m3: ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 ! Diffusivity for DEPAC relevant gases ! Use default value diffc = 0.11e-4 ! overwrite with known coefficients of diffusivity from Massman (1998) IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 ! Get quasi-laminar boundary layer resistance Rb: CALL get_rb_cell( (luw_dep == ilu_water_sea) .or. (luw_dep == ilu_water_inland), & z0h_lu, ustar_lu, diffc, & Rb) ! Get Rc_tot CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & r_aero_lu , Rb) ! Calculate budget IF (Rc_tot .le. 0.0) THEN bud_now_luw(lsp) = 0.0_wp ELSE ! Compute exchange velocity for current lu: vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) bud_now_luw(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ENDIF ENDIF ENDDO ENDIF bud_now = 0.0_wp ! Calculate budget for surface m and adapt concentration DO lsp = 1, nspec bud_now(lsp) = surf_lsm_h%frac(ind_veg_wall,m)*bud_now_lu(lsp) + & surf_lsm_h%frac(ind_pav_green,m)*bud_now_lup(lsp) + & surf_lsm_h%frac(ind_wat_win,m)*bud_now_luw(lsp) ! Compute new concentration, add contribution of all exchange processes: cc(lsp) = cc(lsp) + bud_now(lsp)*inv_dh chem_species(lsp)%conc(k,j,i) = max(0.0_wp, cc(lsp)) ENDDO ENDIF !!!!!!! ! USM ! !!!!!!! IF ( match_usm ) THEN ! Get surface element information at i,j: m = surf_usm_h%start_index(j,i) k = surf_usm_h%k(m) ! Get needed variables for surface element m ustar_lu = surf_usm_h%us(m) z0h_lu = surf_usm_h%z0h(m) r_aero_lu = surf_usm_h%r_a(m) glrad = surf_usm_h%rad_sw_in(m) lai = surf_usm_h%lai(m) sai = lai + 1 ! For small grid spacing neglect R_a IF (dzw(k) <= 1.0) THEN r_aero_lu = 0.0_wp ENDIF !Initialize lu's luu_palm = 0 luu_dep = 0 lug_palm = 0 lug_dep = 0 lud_palm = 0 lud_dep = 0 !Initialize budgets bud_now_luu = 0.0_wp bud_now_lug = 0.0_wp bud_now_lud = 0.0_wp ! Get land use for i,j and assign to DEPAC lu IF (surf_usm_h%frac(ind_pav_green,m) > 0) THEN ! For green urban surfaces (e.g. green roofs ! assume LU short grass lug_palm = ind_lu_s_grass IF (lug_palm == ind_lu_user) THEN message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation' CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 ) ELSEIF (lug_palm == ind_lu_b_soil) THEN lug_dep = 9 ELSEIF (lug_palm == ind_lu_mixed_crops) THEN lug_dep = 2 ELSEIF (lug_palm == ind_lu_s_grass) THEN lug_dep = 1 ELSEIF (lug_palm == ind_lu_ev_needle_trees) THEN lug_dep = 4 ELSEIF (lug_palm == ind_lu_de_needle_trees) THEN lug_dep = 4 ELSEIF (lug_palm == ind_lu_ev_broad_trees) THEN lug_dep = 12 ELSEIF (lug_palm == ind_lu_de_broad_trees) THEN lug_dep = 5 ELSEIF (lug_palm == ind_lu_t_grass) THEN lug_dep = 1 ELSEIF (lug_palm == ind_lu_desert) THEN lug_dep = 9 ELSEIF (lug_palm == ind_lu_tundra ) THEN lug_dep = 8 ELSEIF (lug_palm == ind_lu_irr_crops) THEN lug_dep = 2 ELSEIF (lug_palm == ind_lu_semidesert) THEN lug_dep = 8 ELSEIF (lug_palm == ind_lu_ice) THEN lug_dep = 10 ELSEIF (lug_palm == ind_lu_marsh) THEN lug_dep = 8 ELSEIF (lug_palm == ind_lu_ev_shrubs) THEN lug_dep = 14 ELSEIF (lug_palm == ind_lu_de_shrubs ) THEN lug_dep = 14 ELSEIF (lug_palm == ind_lu_mixed_forest) THEN lug_dep = 4 ELSEIF (lug_palm == ind_lu_intrup_forest) THEN lug_dep = 8 ENDIF ENDIF IF (surf_usm_h%frac(ind_veg_wall,m) > 0) THEN ! For walls in USM assume concrete walls/roofs, ! assumed LU class desert as also assumed for ! pavements in LSM luu_palm = ind_lup_conc IF (luu_palm == ind_lup_user) THEN message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation' CALL message( 'chem_depo', 'CM0455', 1, 2, 0, 6, 0 ) ELSEIF (luu_palm == ind_lup_asph_conc) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_asph) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_conc) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_sett) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_pav_stones) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_cobblest) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_metal) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_wood) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_gravel) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_f_gravel) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_pebblest) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_woodchips) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_tartan) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_art_turf) THEN luu_dep = 9 ELSEIF (luu_palm == ind_lup_clay) THEN luu_dep = 9 ENDIF ENDIF IF (surf_usm_h%frac(ind_wat_win,m) > 0) THEN ! For windows in USM assume metal as this is ! as close as we get, assumed LU class desert ! as also assumed for pavements in LSM lud_palm = ind_lup_metal IF (lud_palm == ind_lup_user) THEN message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation' CALL message( 'chem_depo', 'CM0456', 1, 2, 0, 6, 0 ) ELSEIF (lud_palm == ind_lup_asph_conc) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_asph) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_conc) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_sett) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_pav_stones) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_cobblest) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_metal) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_wood) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_gravel) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_f_gravel) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_pebblest) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_woodchips) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_tartan) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_art_turf) THEN lud_dep = 9 ELSEIF (lud_palm == ind_lup_clay) THEN lud_dep = 9 ENDIF ENDIF ! TODO: Activate these lines as soon as new ebsolver branch is merged: ! Set wetness indicator to dry or wet for usm vegetation or pavement !IF (surf_usm_h%c_liq(m) > 0) THEN ! nwet = 1 !ELSE nwet = 0 !ENDIF ! Compute length of time step IF ( call_chem_at_all_substeps ) THEN dt_chem = dt_3d * weight_pres(intermediate_timestep_count) ELSE dt_chem = dt_3d ENDIF dh = dzw(k) inv_dh = 1.0_wp / dh dt_dh = dt_chem/dh ! Concentration at i,j,k DO lsp = 1, nspec cc(lsp) = chem_species(lsp)%conc(k,j,i) ENDDO ! Temperature at i,j,k ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp ! Surface temperature in degrees Celcius: ts = ttemp - 273.15 ! in degrees celcius ! Viscosity of air in lowest layer visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0) ! Air density in lowest layer dens = rho_air_zw(k) ! Calculate relative humidity from specific humidity for DEPAC qv_tmp = max(q(k,j,i),0.0_wp) relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(nzb) ) ! Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget ! for each surface fraction. Then derive overall budget taking into account the surface fractions. ! Walls/roofs IF (surf_usm_h%frac(ind_veg_wall,m) > 0) THEN ! No vegetation on non-green walls: lai = 0.0_wp sai = 0.0_wp slinnfac = 1.0_wp ! Get vd DO lsp = 1, nvar !Initialize vs = 0.0_wp vd_lu = 0.0_wp Rs = 0.0_wp Rb = 0.0_wp Rc_tot = 0.0_wp IF (spc_names(lsp) == 'PM10') THEN part_type = 1 ! sedimentation velocity in lowest layer vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & luu_dep, & r_aero_lu, ustar_lu) bud_now_luu(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSEIF (spc_names(lsp) == 'PM25') THEN part_type = 2 ! sedimentation velocity in lowest layer: vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & luu_dep , & r_aero_lu, ustar_lu) bud_now_luu(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSE ! GASES ! Read spc_name of current species for gas parameter IF (ANY(pspecnames(:) == spc_names(lsp))) THEN i_pspec = 0 DO pspec = 1, 69 IF (pspecnames(pspec) == spc_names(lsp)) THEN i_pspec = pspec END IF ENDDO ELSE ! Species for now not deposited CYCLE ENDIF ! Factor used for conversion from ppb to ug/m3 : ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) ! c 1e-9 xm_tracer 1e9 / xm_air dens ! thus: ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 ! Use density at lowest layer: ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 ! atmospheric concentration in DEPAC is requested in ug/m3: ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 ! Diffusivity for DEPAC relevant gases ! Use default value diffc = 0.11e-4 ! overwrite with known coefficients of diffusivity from Massman (1998) IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 ! Get quasi-laminar boundary layer resistance Rb: CALL get_rb_cell( (luu_dep == ilu_water_sea) .or. (luu_dep == ilu_water_inland), & z0h_lu, ustar_lu, diffc, & Rb) ! Get Rc_tot CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & r_aero_lu , Rb) ! Calculate budget IF (Rc_tot .le. 0.0) THEN bud_now_luu(lsp) = 0.0_wp ELSE ! Compute exchange velocity for current lu: vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) bud_now_luu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ENDIF ENDIF ENDDO ENDIF ! Green usm surfaces IF (surf_usm_h%frac(ind_pav_green,m) > 0) THEN slinnfac = 1.0_wp ! Get vd DO lsp = 1, nvar !Initialize vs = 0.0_wp vd_lu = 0.0_wp Rs = 0.0_wp Rb = 0.0_wp Rc_tot = 0.0_wp IF (spc_names(lsp) == 'PM10') THEN part_type = 1 ! sedimentation velocity in lowest layer: vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & lug_dep, & r_aero_lu, ustar_lu) bud_now_lug(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSEIF (spc_names(lsp) == 'PM25') THEN part_type = 2 ! sedimentation velocity in lowest layer: vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & lug_dep, & r_aero_lu, ustar_lu) bud_now_lug(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSE ! GASES ! Read spc_name of current species for gas parameter IF (ANY(pspecnames(:) == spc_names(lsp))) THEN i_pspec = 0 DO pspec = 1, 69 IF (pspecnames(pspec) == spc_names(lsp)) THEN i_pspec = pspec END IF ENDDO ELSE ! Species for now not deposited CYCLE ENDIF ! Factor used for conversion from ppb to ug/m3 : ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) ! c 1e-9 xm_tracer 1e9 / xm_air dens ! thus: ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 ! Use density at lowest layer: ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 ! atmospheric concentration in DEPAC is requested in ug/m3: ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 ! Diffusivity for DEPAC relevant gases ! Use default value diffc = 0.11e-4 ! overwrite with known coefficients of diffusivity from Massman (1998) IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 ! Get quasi-laminar boundary layer resistance Rb: CALL get_rb_cell( (lug_dep == ilu_water_sea) .or. (lug_dep == ilu_water_inland), & z0h_lu, ustar_lu, diffc, & Rb) ! Get Rc_tot CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & r_aero_lu , Rb) ! Calculate budget IF (Rc_tot .le. 0.0) THEN bud_now_lug(lsp) = 0.0_wp ELSE ! Compute exchange velocity for current lu: vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) bud_now_lug(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ENDIF ENDIF ENDDO ENDIF ! Windows IF (surf_usm_h%frac(ind_wat_win,m) > 0) THEN ! No vegetation on windows: lai = 0.0_wp sai = 0.0_wp slinnfac = 1.0_wp ! Get vd DO lsp = 1, nvar !Initialize vs = 0.0_wp vd_lu = 0.0_wp Rs = 0.0_wp Rb = 0.0_wp Rc_tot = 0.0_wp IF (spc_names(lsp) == 'PM10') THEN part_type = 1 ! sedimentation velocity in lowest layer: vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & lud_dep, & r_aero_lu, ustar_lu) bud_now_lud(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSEIF (spc_names(lsp) == 'PM25') THEN part_type = 2 ! sedimentation velocity in lowest layer: vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & visc) CALL drydepo_aero_zhang_vd( vd_lu, Rs, & vs, & particle_pars(ind_p_size, part_type), & particle_pars(ind_p_slip, part_type), & nwet, ttemp, dens, visc, & lud_dep, & r_aero_lu, ustar_lu) bud_now_lud(lsp) = - cc(lsp) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ELSE ! GASES ! Read spc_name of current species for gas PARAMETER IF (ANY(pspecnames(:) == spc_names(lsp))) THEN i_pspec = 0 DO pspec = 1, 69 IF (pspecnames(pspec) == spc_names(lsp)) THEN i_pspec = pspec END IF ENDDO ELSE ! Species for now not deposited CYCLE ENDIF ! Factor used for conversion from ppb to ug/m3 : ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) ! c 1e-9 xm_tracer 1e9 / xm_air dens ! thus: ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 ! Use density at lowest layer: ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 ! atmospheric concentration in DEPAC is requested in ug/m3: ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 ! Diffusivity for DEPAC relevant gases ! Use default value diffc = 0.11e-4 ! overwrite with known coefficients of diffusivity from Massman (1998) IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 ! Get quasi-laminar boundary layer resistance Rb: CALL get_rb_cell( (lud_dep == ilu_water_sea) .or. (lud_dep == ilu_water_inland), & z0h_lu, ustar_lu, diffc, & Rb) ! Get Rc_tot CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & r_aero_lu , Rb) ! Calculate budget IF (Rc_tot .le. 0.0) THEN bud_now_lud(lsp) = 0.0_wp ELSE ! Compute exchange velocity for current lu: vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) bud_now_lud(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & (1.0 - exp(-vd_lu*dt_dh ))*dh ENDIF ENDIF ENDDO ENDIF bud_now = 0.0_wp ! Calculate budget for surface m and adapt concentration DO lsp = 1, nspec bud_now(lsp) = surf_usm_h%frac(ind_veg_wall,m)*bud_now_luu(lsp) + & surf_usm_h%frac(ind_pav_green,m)*bud_now_lug(lsp) + & surf_usm_h%frac(ind_wat_win,m)*bud_now_lud(lsp) ! Compute new concentration, add contribution of all exchange processes: cc(lsp) = cc(lsp) + bud_now(lsp)*inv_dh chem_species(lsp)%conc(k,j,i) = max(0.0_wp, cc(lsp)) ENDDO ENDIF END SUBROUTINE chem_depo !---------------------------------------------------------------------------------- ! ! DEPAC: ! Code of the DEPAC routine and corresponding subroutines below from the DEPAC ! module of the LOTOS-EUROS model (Manders et al., 2017) ! ! Original DEPAC routines by RIVM and TNO (2015), for Documentation see ! van Zanten et al., 2010. !--------------------------------------------------------------------------------- ! !---------------------------------------------------------------------------------- ! ! depac: compute total canopy (or surface) resistance Rc for gases !---------------------------------------------------------------------------------- SUBROUTINE drydepos_gas_depac(compnam, day_of_year, lat, t, ust, glrad, sinphi, & rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, catm, diffc, & ra, rb) ! The last four rows of depac arguments are OPTIONAL: ! ! A. compute Rc_tot without compensation points (ccomp_tot will be zero): ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi]) ! ! B. compute Rc_tot with compensation points (used for LOTOS-EUROS): ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water) ! ! C. compute effective Rc based on compensation points (used for OPS): ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & ! ra, rb, rc_eff) ! X1. Extra (OPTIONAL) output variables: ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & ! ra, rb, rc_eff, & ! gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out) ! X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves): ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & ! ra, rb, rc_eff, & ! gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, & ! calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA) IMPLICIT NONE CHARACTER(len=*) , INTENT(in) :: compnam ! component name ! 'HNO3','NO','NO2','O3','SO2','NH3' INTEGER(iwp) , INTENT(in) :: day_of_year ! day of year, 1 ... 365 (366) REAL(kind=wp) , INTENT(in) :: lat ! latitude Northern hemisphere (degrees) (DEPAC cannot be used for S. hemisphere) REAL(kind=wp) , INTENT(in) :: t ! temperature (C) ! NB discussion issue is temp T_2m or T_surf or T_leaf? REAL(kind=wp) , INTENT(in) :: ust ! friction velocity (m/s) REAL(kind=wp) , INTENT(in) :: glrad ! global radiation (W/m2) REAL(kind=wp) , INTENT(in) :: sinphi ! sin of solar elevation angle REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) REAL(kind=wp) , INTENT(in) :: lai ! one-sidedleaf area index (-) REAL(kind=wp) , INTENT(in) :: sai ! surface area index (-) (lai + branches and stems) REAL(kind=wp) , INTENT(in) :: diffc ! Diffusivity INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow INTEGER(iwp) , INTENT(in) :: lu ! land use type, lu = 1,...,nlu INTEGER(iwp) , INTENT(in) :: iratns ! index for NH3/SO2 ratio used for SO2; ! iratns = 1: low NH3/SO2 ! iratns = 2: high NH3/SO2 ! iratns = 3: very low NH3/SO2 REAL(kind=wp) , INTENT(out) :: rc_tot ! total canopy resistance Rc (s/m) REAL(kind=wp) , INTENT(out) :: ccomp_tot ! total compensation point (ug/m3) [= 0 for species that don't have a compensation REAL(kind=wp) ,INTENT(in) :: p ! pressure (Pa) REAL(kind=wp), INTENT(in) :: catm ! actual atmospheric concentration (ug/m3) REAL(kind=wp), INTENT(in) :: ra ! aerodynamic resistance (s/m) REAL(kind=wp), INTENT(in) :: rb ! boundary layer resistance (s/m) ! Local variables: REAL(kind=wp) :: laimax ! maximum leaf area index (-) LOGICAL :: ready ! Rc has been set ! = 1 -> constant Rc ! = 2 -> temperature dependent Rc REAL(kind=wp) :: gw ! external leaf conductance (m/s) REAL(kind=wp) :: gstom ! stomatal conductance (m/s) REAL(kind=wp) :: gsoil_eff ! effective soil conductance (m/s) REAL(kind=wp) :: gc_tot ! total canopy conductance (m/s) REAL(kind=wp) :: cw ! external leaf surface compensation point (ug/m3) REAL(kind=wp) :: cstom ! stomatal compensation point (ug/m3) REAL(kind=wp) :: csoil ! soil compensation point (ug/m3) ! Vegetation indicators: LOGICAL :: LAI_present ! leaves are present for current land use type LOGICAL :: SAI_present ! vegetation is present for current land use type ! Component number taken from component name, paramteres matched with include files INTEGER(iwp) :: icmp ! component numbers: INTEGER(iwp), PARAMETER :: icmp_o3 = 1 INTEGER(iwp), PARAMETER :: icmp_so2 = 2 INTEGER(iwp), PARAMETER :: icmp_no2 = 3 INTEGER(iwp), PARAMETER :: icmp_no = 4 INTEGER(iwp), PARAMETER :: icmp_nh3 = 5 INTEGER(iwp), PARAMETER :: icmp_co = 6 INTEGER(iwp), PARAMETER :: icmp_no3 = 7 INTEGER(iwp), PARAMETER :: icmp_hno3 = 8 INTEGER(iwp), PARAMETER :: icmp_n2o5 = 9 INTEGER(iwp), PARAMETER :: icmp_h2o2 = 10 ! Define component number SELECT CASE ( trim(compnam) ) CASE ( 'O3', 'o3' ) icmp = icmp_o3 CASE ( 'SO2', 'so2' ) icmp = icmp_so2 CASE ( 'NO2', 'no2' ) icmp = icmp_no2 CASE ( 'NO', 'no' ) icmp = icmp_no CASE ( 'NH3', 'nh3' ) icmp = icmp_nh3 CASE ( 'CO', 'co' ) icmp = icmp_co CASE ( 'NO3', 'no3' ) icmp = icmp_no3 CASE ( 'HNO3', 'hno3' ) icmp = icmp_hno3 CASE ( 'N2O5', 'n2o5' ) icmp = icmp_n2o5 CASE ( 'H2O2', 'h2o2' ) icmp = icmp_h2o2 CASE default !Component not part of DEPAC --> not deposited RETURN END SELECT ! inititalize gw = 0. gstom = 0. gsoil_eff = 0. gc_tot = 0. cw = 0. cstom = 0. csoil = 0. ! Check whether vegetation is present (in that CASE the leaf or surface area index > 0): LAI_present = (lai .gt. 0.0) SAI_present = (sai .gt. 0.0) ! Set Rc (i.e. rc_tot) in special cases: CALL rc_special(icmp,compnam,lu,t, nwet,rc_tot,ready,ccomp_tot) ! set effective resistance equal to total resistance, this will be changed for compensation points !IF ( present(rc_eff) ) then ! rc_eff = rc_tot !END IF ! IF Rc is not set: IF (.not. ready) then ! External conductance: CALL rc_gw(compnam,iratns,t,rh,nwet,SAI_present,sai,gw) ! Stomatal conductance: ! CALL rc_gstom(icmp,compnam,lu,LAI_present,lai,glrad,sinphi,t,rh,diffc,gstom,p,& ! smi=smi,calc_stom_o3flux=calc_stom_o3flux) CALL rc_gstom(icmp,compnam,lu,LAI_present,lai,glrad,sinphi,t,rh,diffc,gstom,p) ! Effective soil conductance: CALL rc_gsoil_eff(icmp,lu,sai,ust,nwet,t,gsoil_eff) ! Total canopy conductance (gc_tot) and resistance Rc (rc_tot): CALL rc_rctot(gstom,gsoil_eff,gw,gc_tot,rc_tot) ! Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3): ! CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,& ! LAI_present, SAI_present, & ! ccomp_tot, & ! catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, & ! c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, & ! tsea=tsea,cw=cw,cstom=cstom,csoil=csoil ) ! Effective Rc based on compensation points: ! IF ( present(rc_eff) ) then ! check on required arguments: !IF ( (.not. present(catm)) .or. (.not. present(ra)) .or. (.not. present(rb)) ) then ! stop 'output argument rc_eff requires input arguments catm, ra and rb' !END IF ! compute rc_eff : !CALL rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff) !ENDIF ENDIF END SUBROUTINE drydepos_gas_depac !------------------------------------------------------------------- ! rc_special: compute total canopy resistance in special CASEs !------------------------------------------------------------------- SUBROUTINE rc_special(icmp,compnam,lu,t,nwet,rc_tot,ready,ccomp_tot) INTEGER(iwp) , INTENT(in) :: icmp ! component index CHARACTER(len=*) , INTENT(in) :: compnam ! component name INTEGER(iwp) , INTENT(in) :: lu ! land use type, lu = 1,...,nlu REAL(kind=wp) , INTENT(in) :: t ! temperature (C) ! = 1 -> constant Rc ! = 2 -> temperature dependent Rc INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow REAL(kind=wp) , INTENT(out) :: rc_tot ! total canopy resistance Rc (s/m) LOGICAL , INTENT(out) :: ready ! Rc has been set REAL(kind=wp) , INTENT(out) :: ccomp_tot ! total compensation point (ug/m3) ! rc_tot is not yet set: ready = .false. ! Default compensation point in special CASEs = 0: ccomp_tot = 0.0 SELECT CASE(trim(compnam)) CASE('HNO3','N2O5','NO3','H2O2') ! No separate resistances for HNO3; just one total canopy resistance: IF (t .lt. -5.0 .and. nwet .eq. 9) then ! T < 5 C and snow: rc_tot = 50. ELSE ! all other circumstances: rc_tot = 10.0 ENDIF ready = .true. CASE('NO','CO') IF (lu .eq. ilu_water_sea .or. lu .eq. ilu_water_inland) then ! water rc_tot = 2000. ready = .true. ELSEIF (nwet .eq. 1) then ! wet rc_tot = 2000. ready = .true. ENDIF CASE('NO2','O3','SO2','NH3') ! snow surface: IF (nwet.eq.9) then !CALL rc_snow(ipar_snow(icmp),t,rc_tot) ready = .true. ENDIF CASE default message_string = 'Component '// trim(compnam) // ' not supported' CALL message( 'rc_special', 'CM0457', 1, 2, 0, 6, 0 ) END SELECT END SUBROUTINE rc_special !------------------------------------------------------------------- ! rc_gw: compute external conductance !------------------------------------------------------------------- SUBROUTINE rc_gw( compnam, iratns,t,rh,nwet, SAI_present, sai, gw ) ! Input/output variables: CHARACTER(len=*) , INTENT(in) :: compnam ! component name INTEGER(iwp) , INTENT(in) :: iratns ! index for NH3/SO2 ratio; ! iratns = 1: low NH3/SO2 ! iratns = 2: high NH3/SO2 ! iratns = 3: very low NH3/SO2 REAL(kind=wp) , INTENT(in) :: t ! temperature (C) REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow LOGICAL , INTENT(in) :: SAI_present REAL(kind=wp) , INTENT(in) :: sai ! one-sided leaf area index (-) REAL(kind=wp) , INTENT(out) :: gw ! external leaf conductance (m/s) SELECT CASE(trim(compnam)) !CASE('HNO3','N2O5','NO3','H2O2') this routine is not CALLed for HNO3 CASE('NO2') CALL rw_constant(2000.0_wp,SAI_present,gw) CASE('NO','CO') CALL rw_constant(-9999.0_wp,SAI_present,gw) ! see Erisman et al, 1994 section 3.2.3 CASE('O3') CALL rw_constant(2500.0_wp,SAI_present,gw) CASE('SO2') CALL rw_so2( t, nwet, rh, iratns, SAI_present, gw ) CASE('NH3') CALL rw_nh3_sutton(t,rh,SAI_present,gw) ! conversion from leaf resistance to canopy resistance by multiplying with SAI: Gw = sai*gw CASE default message_string = 'Component '// trim(compnam) // ' not supported' CALL message( 'rc_gw', 'CM0458', 1, 2, 0, 6, 0 ) END SELECT END SUBROUTINE rc_gw !------------------------------------------------------------------- ! rw_so2: compute external leaf conductance for SO2 !------------------------------------------------------------------- SUBROUTINE rw_so2( t, nwet, rh, iratns, SAI_present, gw ) ! Input/output variables: REAL(kind=wp) , INTENT(in) :: t ! temperature (C) INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) INTEGER(iwp) , INTENT(in) :: iratns ! index for NH3/SO2 ratio; ! iratns = 1: low NH3/SO2 ! iratns = 2: high NH3/SO2 ! iratns = 3: very low NH3/SO2 LOGICAL, INTENT(in) :: SAI_present REAL(kind=wp) , INTENT(out) :: gw ! external leaf conductance (m/s) ! Variables from module: ! SAI_present: vegetation is present ! Local variables: REAL(kind=wp) :: rw ! external leaf resistance (s/m) ! Check IF vegetation present: IF (SAI_present) then IF (nwet .eq. 0) then !-------------------------- ! dry surface !-------------------------- ! T > -1 C IF (t .gt. -1.0) then IF (rh .lt. 81.3) then rw = 25000*exp(-0.0693*rh) ELSE rw = 0.58e12*exp(-0.278*rh) + 10. ENDIF ELSE ! -5 C < T <= -1 C IF (t .gt. -5.0) then rw=200 ELSE ! T <= -5 C rw=500 ENDIF ENDIF ELSE !-------------------------- ! wet surface !-------------------------- rw = 10. !see Table 5, Erisman et al, 1994 Atm. Environment, 0 is impl. as 10 ENDIF ! very low NH3/SO2 ratio: IF (iratns == iratns_very_low ) rw = rw+50. ! Conductance: gw = 1./rw ELSE ! no vegetation: gw = 0.0 ENDIF END SUBROUTINE rw_so2 !------------------------------------------------------------------- ! rw_nh3_sutton: compute external leaf conductance for NH3, ! following Sutton & Fowler, 1993 !------------------------------------------------------------------- SUBROUTINE rw_nh3_sutton(tsurf,rh,SAI_present,gw) ! Input/output variables: REAL(kind=wp) , INTENT(in) :: tsurf ! surface temperature (C) REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) LOGICAL, INTENT(in) :: SAI_present REAL(kind=wp) , INTENT(out) :: gw ! external leaf conductance (m/s) ! Variables from module: ! SAI_present: vegetation is present ! Local variables: REAL(kind=wp) :: rw ! external leaf resistance (s/m) REAL(kind=wp) :: sai_grass_haarweg ! surface area index at experimental site Haarweg ! Fix SAI_grass at value valid for Haarweg data for which gamma_w parametrization is derived sai_grass_haarweg = 3.5 ! 100 - rh ! rw = 2.0 * exp(----------) ! 12 IF (SAI_present) then ! External resistance according to Sutton & Fowler, 1993 rw = 2.0 * exp((100.0 - rh)/12.0) rw = sai_grass_haarweg * rw ! Frozen soil (from Depac v1): IF (tsurf .lt. 0) rw = 200 ! Conductance: gw = 1./rw ELSE ! no vegetation: gw = 0.0 ENDIF END SUBROUTINE rw_nh3_sutton !------------------------------------------------------------------- ! rw_constant: compute constant external leaf conductance !------------------------------------------------------------------- SUBROUTINE rw_constant(rw_val,SAI_present,gw) ! Input/output variables: REAL(kind=wp) , INTENT(in) :: rw_val ! constant value of Rw LOGICAL , INTENT(in) :: SAI_present REAL(kind=wp) , INTENT(out) :: gw ! wernal leaf conductance (m/s) ! Variables from module: ! SAI_present: vegetation is present ! Compute conductance: IF (SAI_present .and. .not.missing(rw_val)) then gw = 1./rw_val ELSE gw = 0. ENDIF END SUBROUTINE rw_constant !------------------------------------------------------------------- ! rc_gstom: compute stomatal conductance !------------------------------------------------------------------- SUBROUTINE rc_gstom( icmp, compnam, lu, LAI_present, lai, glrad, sinphi, t, rh, diffc, & gstom, & p) ! input/output variables: INTEGER(iwp), INTENT(in) :: icmp ! component index CHARACTER(len=*), INTENT(in) :: compnam ! component name INTEGER(iwp), INTENT(in) :: lu ! land use type , lu = 1,...,nlu LOGICAL, INTENT(in) :: LAI_present REAL(kind=wp), INTENT(in) :: lai ! one-sided leaf area index REAL(kind=wp), INTENT(in) :: glrad ! global radiation (W/m2) REAL(kind=wp), INTENT(in) :: sinphi ! sin of solar elevation angle REAL(kind=wp), INTENT(in) :: t ! temperature (C) REAL(kind=wp), INTENT(in) :: rh ! relative humidity (%) REAL(kind=wp), INTENT(in) :: diffc ! diffusion coefficient of the gas involved REAL(kind=wp), INTENT(out):: gstom ! stomatal conductance (m/s) REAL(kind=wp), OPTIONAL,INTENT(in) :: p ! pressure (Pa) ! Local variables REAL(kind=wp) :: vpd ! vapour pressure deficit (kPa) REAL(kind=wp), PARAMETER :: dO3 = 0.13e-4 ! diffusion coefficient of ozon (m2/s) SELECT CASE(trim(compnam)) !CASE('HNO3','N2O5','NO3','H2O2') this routine is not CALLed for HNO3 CASE('NO','CO') ! for NO stomatal uptake is neglected: gstom = 0.0 CASE('NO2','O3','SO2','NH3') ! IF vegetation present: IF (LAI_present) then IF (glrad .gt. 0.0) then CALL rc_get_vpd(t,rh,vpd) CALL rc_gstom_emb( lu, glrad, t, vpd, LAI_present, lai, sinphi, gstom, p ) gstom = gstom*diffc/dO3 ! Gstom of Emberson is derived for ozone ELSE gstom = 0.0 ENDIF ELSE ! no vegetation; zero conductance (infinite resistance): gstom = 0.0 ENDIF CASE default message_string = 'Component '// trim(compnam) // ' not supported' CALL message( 'rc_gstom', 'CM0459', 1, 2, 0, 6, 0 ) END SELECT END SUBROUTINE rc_gstom !------------------------------------------------------------------- ! rc_gstom_emb: stomatal conductance according to Emberson !------------------------------------------------------------------- SUBROUTINE rc_gstom_emb( lu, glrad, T, vpd, LAI_present, lai, sinp, & Gsto, p ) ! History ! Original code from Lotos-Euros, TNO, M. Schaap ! 2009-08, M.C. van Zanten, Rivm ! Updated and extended. ! 2009-09, Arjo Segers, TNO ! Limitted temperature influence to range to avoid ! floating point exceptions. ! ! Method ! ! Code based on Emberson et al, 2000, Env. Poll., 403-413 ! Notation conform Unified EMEP Model Description Part 1, ch 8 ! ! In the calculation of F_light the modification of L. Zhang 2001, AE to the PARshade and PARsun ! parametrizations of Norman 1982 are applied ! ! F_phen and F_SWP are set to 1 ! ! Land use types DEPAC versus Emberson (Table 5.1, EMEP model description) ! DEPAC Emberson ! 1 = grass GR = grassland ! 2 = arable land TC = temperate crops ( LAI according to RC = rootcrops) ! 3 = permanent crops TC = temperate crops ( LAI according to RC = rootcrops) ! 4 = coniferous forest CF = tempareate/boREAL(kind=wp) coniferous forest ! 5 = deciduous forest DF = temperate/boREAL(kind=wp) deciduous forest ! 6 = water W = water ! 7 = urban U = urban ! 8 = other GR = grassland ! 9 = desert DE = desert ! ! Emberson specific declarations ! Input/output variables: INTEGER(iwp), INTENT(in) :: lu ! land use type, lu = 1,...,nlu REAL(kind=wp), INTENT(in) :: glrad ! global radiation (W/m2) REAL(kind=wp), INTENT(in) :: T ! temperature (C) REAL(kind=wp), INTENT(in) :: vpd ! vapour pressure deficit (kPa) LOGICAL, INTENT(in) :: LAI_present REAL(kind=wp), INTENT(in) :: lai ! one-sided leaf area index REAL(kind=wp), INTENT(in) :: sinp ! sin of solar elevation angle REAL(kind=wp), INTENT(out):: Gsto ! stomatal conductance (m/s) REAL(kind=wp), OPTIONAL, INTENT(in) :: p ! pressure (Pa) ! Local variables: REAL(kind=wp) :: F_light, F_phen, F_temp, F_vpd, F_swp, bt, F_env REAL(kind=wp) :: pardir, pardiff, parshade, parsun, LAIsun, LAIshade, sinphi REAL(kind=wp) :: pres REAL(kind=wp), PARAMETER :: p_sealevel = 1.01325e5 ! Pa ! Check whether vegetation is present: IF (LAI_present) then ! calculation of correction factors for stomatal conductance IF (sinp .le. 0.0) then sinphi = 0.0001 ELSE sinphi = sinp END IF ! ratio between actual and sea-level pressure is used ! to correct for height in the computation of par ; ! should not exceed sea-level pressure therefore ... IF ( present(p) ) then pres = min( p, p_sealevel ) ELSE pres = p_sealevel ENDIF ! Direct and diffuse par, Photoactive (=visible) radiation: CALL par_dir_diff( glrad, sinphi, pres, p_sealevel, pardir, pardiff ) ! par for shaded leaves (canopy averaged): parshade = pardiff*exp(-0.5*LAI**0.7)+0.07*pardir*(1.1-0.1*LAI)*exp(-sinphi) ! Norman, 1982 IF (glrad .gt. 200 .and. LAI .gt. 2.5) then parshade = pardiff*exp(-0.5*LAI**0.8)+0.07*pardir*(1.1-0.1*LAI)*exp(-sinphi) ! Zhang et al., 2001 END IF ! par for sunlit leaves (canopy averaged): ! alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5 parsun = pardir*0.5/sinphi + parshade ! Norman, 1982 IF (glrad .gt. 200 .and. LAI .gt. 2.5) then parsun = pardir**0.8*0.5/sinphi + parshade ! Zhang et al., 2001 END IF ! leaf area index for sunlit and shaded leaves: IF (sinphi .gt. 0) then LAIsun = 2*sinphi*(1-exp(-0.5*LAI/sinphi )) LAIshade = LAI -LAIsun ELSE LAIsun = 0 LAIshade = LAI END IF F_light = (LAIsun*(1 - exp(-1.*alpha(lu)*parsun)) + LAIshade*(1 - exp(-1.*alpha(lu)*parshade)))/LAI F_light = max(F_light,F_min(lu)) ! temperature influence; only non-zero within range [Tmin,Tmax]: IF ( (Tmin(lu) < T) .and. (T < Tmax(lu)) ) then BT = (Tmax(lu)-Topt(lu))/(Topt(lu)-Tmin(lu)) F_temp = ((T-Tmin(lu))/(Topt(lu)-Tmin(lu))) * ((Tmax(lu)-T)/(Tmax(lu)-Topt(lu)))**BT ELSE F_temp = 0.0 END IF F_temp = max(F_temp,F_min(lu)) ! vapour pressure deficit influence F_vpd = MIN( 1.0_wp, ((1.0_wp-F_min(lu))*(vpd_min(lu)-vpd)/(vpd_min(lu)-vpd_max(lu)) + F_min(lu) ) ) F_vpd = max(F_vpd,F_min(lu)) F_swp = 1. ! influence of phenology on stom. conductance ! ignored for now in DEPAC since influence of F_phen on lu classes in use is negligible. ! When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore F_phen = 1. ! evaluate total stomatal conductance F_env = F_temp*F_vpd*F_swp F_env = max(F_env,F_min(lu)) gsto = G_max(lu) * F_light * F_phen * F_env ! gstom expressed per m2 leafarea; ! this is converted with LAI to m2 surface. Gsto = lai*gsto ! in m/s ELSE GSto = 0.0 ENDIF END SUBROUTINE rc_gstom_emb !------------------------------------------------------------------- ! par_dir_diff !------------------------------------------------------------------- SUBROUTINE par_dir_diff(glrad,sinphi,pres,pres_0,par_dir,par_diff) ! ! Weiss, A., Norman, J.M. (1985) Partitioning solar radiation into direct and ! diffuse, visible and near-infrared components. Agric. Forest Meteorol. ! 34, 205-213. ! ! From a SUBROUTINE obtained from Leiming Zhang (via Willem Asman), ! MeteoroLOGICAL Service of Canada ! e-mail: leiming.zhang@ec.gc.ca ! ! Leiming uses solar irradiance. This should be equal to global radiation and ! Willem Asman set it to global radiation REAL(kind=wp), INTENT(in) :: glrad ! global radiation (W m-2) REAL(kind=wp), INTENT(in) :: sinphi ! sine of the solar elevation REAL(kind=wp), INTENT(in) :: pres ! actual pressure (to correct for height) (Pa) REAL(kind=wp), INTENT(in) :: pres_0 ! pressure at sea level (Pa) REAL(kind=wp), INTENT(out) :: par_dir ! PAR direct : visible (photoactive) direct beam radiation (W m-2) REAL(kind=wp), INTENT(out) :: par_diff ! PAR diffuse: visible (photoactive) diffuse radiation (W m-2) ! fn = near-infrared direct beam fraction (DIMENSIONless) ! fv = PAR direct beam fraction (DIMENSIONless) ! ratio = ratio measured to potential solar radiation (DIMENSIONless) ! rdm = potential direct beam near-infrared radiation (W m-2); "potential" means clear-sky ! rdn = potential diffuse near-infrared radiation (W m-2) ! rdu = visible (PAR) direct beam radiation (W m-2) ! rdv = potential visible (PAR) diffuse radiation (W m-2) ! rn = near-infrared radiation (W m-2) ! rv = visible radiation (W m-2) ! ww = water absorption in the near infrared for 10 mm of precipitable water REAL(kind=wp) :: rdu,rdv,ww,rdm,rdn,rv,rn,ratio,sv,fv ! Calculate visible (PAR) direct beam radiation ! 600 W m-2 represents average amount of PAR (400-700 nm wavelength) ! at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2) rdu=600.*exp(-0.185*(pres/pres_0)/sinphi)*sinphi ! Calculate potential visible diffuse radiation rdv=0.4*(600.- rdu)*sinphi ! Calculate the water absorption in the-near infrared ww=1320*10**( -1.195+0.4459*log10(1./sinphi)-0.0345*(log10(1./sinphi))**2 ) ! Calculate potential direct beam near-infrared radiation rdm=(720.*exp(-0.06*(pres/pres_0)/sinphi)-ww)*sinphi !720 = solar constant - 600 ! Calculate potential diffuse near-infrared radiation rdn=0.6*(720-rdm-ww)*sinphi ! Compute visible and near-infrared radiation rv = MAX( 0.1_wp, rdu+rdv ) rn = MAX( 0.01_wp, rdm+rdn ) ! Compute ratio between input global radiation and total radiation computed here ratio = MIN( 0.9_wp, glrad/(rv+rn) ) ! Calculate total visible radiation sv=ratio*rv ! Calculate fraction of PAR in the direct beam fv = MIN( 0.99_wp, (0.9_wp-ratio)/0.7_wp ) ! help variable fv = MAX( 0.01_wp, rdu/rv*(1.0_wp-fv**0.6667_wp) ) ! fraction of PAR in the direct beam ! Compute direct and diffuse parts of PAR par_dir=fv*sv par_diff=sv-par_dir END SUBROUTINE par_dir_diff !------------------------------------------------------------------- ! rc_get_vpd: get vapour pressure deficit (kPa) !------------------------------------------------------------------- SUBROUTINE rc_get_vpd(temp, relh,vpd) IMPLICIT NONE ! Input/output variables: REAL(kind=wp) , INTENT(in) :: temp ! temperature (C) REAL(kind=wp) , INTENT(in) :: relh ! relative humidity (%) REAL(kind=wp) , INTENT(out) :: vpd ! vapour pressure deficit (kPa) ! Local variables: REAL(kind=wp) :: esat ! fit PARAMETERs: REAL(kind=wp), PARAMETER :: a1 = 6.113718e-1 REAL(kind=wp), PARAMETER :: a2 = 4.43839e-2 REAL(kind=wp), PARAMETER :: a3 = 1.39817e-3 REAL(kind=wp), PARAMETER :: a4 = 2.9295e-5 REAL(kind=wp), PARAMETER :: a5 = 2.16e-7 REAL(kind=wp), PARAMETER :: a6 = 3.0e-9 ! esat is saturation vapour pressure (kPa) at temp(C) following monteith(1973) esat = a1 + a2*temp + a3*temp**2 + a4*temp**3 + a5*temp**4 + a6*temp**5 vpd = esat * (1-relh/100) END SUBROUTINE rc_get_vpd !------------------------------------------------------------------- ! rc_gsoil_eff: compute effective soil conductance !------------------------------------------------------------------- SUBROUTINE rc_gsoil_eff(icmp,lu,sai,ust,nwet,t,gsoil_eff) ! Input/output variables: INTEGER(iwp) , INTENT(in) :: icmp ! component index INTEGER(iwp) , INTENT(in) :: lu ! land use type, lu = 1,...,nlu REAL(kind=wp) , INTENT(in) :: sai ! surface area index REAL(kind=wp) , INTENT(in) :: ust ! friction velocity (m/s) INTEGER(iwp) , INTENT(in) :: nwet ! index for wetness ! nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow ! N.B. this routine cannot be CALLed with nwet = 9, ! nwet = 9 should be handled outside this routine. REAL(kind=wp) , INTENT(in) :: t ! temperature (C) REAL(kind=wp) , INTENT(out) :: gsoil_eff ! effective soil conductance (m/s) ! local variables: REAL(kind=wp) :: rinc ! in canopy resistance (s/m) REAL(kind=wp) :: rsoil_eff ! effective soil resistance (s/m) ! Soil resistance (numbers matched with lu_classes and component numbers) REAL(kind=wp), PARAMETER :: rsoil(nlu_dep,ncmp) = reshape( (/ & ! grs ara crp cnf dec wat urb oth des ice sav trf wai med sem 1000., 200., 200., 200., 200., 2000., 400., 1000., 2000., 2000., 1000., 200., 2000., 200., 400., & ! O3 1000., 1000., 1000., 1000., 1000., 10., 1000., 1000., 1000., 500., 1000., 1000., 10., 1000., 1000., & ! SO2 1000., 1000., 1000., 1000., 1000., 2000., 1000., 1000., 1000., 2000., 1000., 1000., 2000., 1000., 1000., & ! NO2 -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., & ! NO 100., 100., 100., 100., 100., 10., 100., 100., 100., 1000., 100., 100., 10., 100., 100., & ! NH3 -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., & ! CO -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! NO3 -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! HNO3 -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! N2O5 -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999. /),& ! H2O2 (/nlu_dep,ncmp/) ) ! o3 so2 no2 no nh3 co no3 hno3 n2o5 h2o2 REAL(kind=wp), PARAMETER :: rsoil_wet(ncmp) = (/2000., 10., 2000., -999., 10., -999., -999., -999., -999., -999./) REAL(kind=wp), PARAMETER :: rsoil_frozen(ncmp) = (/2000., 500., 2000., -999., 1000., -999., -999., -999., -999., -999./) ! Compute in canopy (in crop) resistance: CALL rc_rinc(lu,sai,ust,rinc) ! Check for missing deposition path: IF (missing(rinc)) then rsoil_eff = -9999. ELSE ! Frozen soil (temperature below 0): IF (t .lt. 0.0) then IF (missing(rsoil_frozen(icmp))) then rsoil_eff = -9999. ELSE rsoil_eff = rsoil_frozen(icmp) + rinc ENDIF ELSE ! Non-frozen soil; dry: IF (nwet .eq. 0) then IF (missing(rsoil(lu,icmp))) then rsoil_eff = -9999. ELSE rsoil_eff = rsoil(lu,icmp) + rinc ENDIF ! Non-frozen soil; wet: ELSEIF (nwet .eq. 1) then IF (missing(rsoil_wet(icmp))) then rsoil_eff = -9999. ELSE rsoil_eff = rsoil_wet(icmp) + rinc ENDIF ELSE message_string = 'nwet can only be 0 or 1' CALL message( 'rc_gsoil_eff', 'CM0460', 1, 2, 0, 6, 0 ) ENDIF ENDIF ENDIF ! Compute conductance: IF (rsoil_eff .gt. 0.0) then gsoil_eff = 1./rsoil_eff ELSE gsoil_eff = 0.0 ENDIF END SUBROUTINE rc_gsoil_eff !------------------------------------------------------------------- ! rc_rinc: compute in canopy (or in crop) resistance ! van Pul and Jacobs, 1993, BLM !------------------------------------------------------------------- SUBROUTINE rc_rinc(lu,sai,ust,rinc) ! Input/output variables: INTEGER(iwp) , INTENT(in) :: lu ! land use class, lu = 1, ..., nlu REAL(kind=wp) , INTENT(in) :: sai ! surface area index REAL(kind=wp) , INTENT(in) :: ust ! friction velocity (m/s) REAL(kind=wp) , INTENT(out) :: rinc ! in canopy resistance (s/m) ! b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 IF not applicable) ! h = vegetation height (m) grass arabl crops conIF decid water urba othr desr ice sav trf wai med semi REAL(kind=wp), DIMENSION(nlu_dep), PARAMETER :: b = (/ -999, 14, 14, 14, 14, -999, -999, -999, -999, -999, -999, 14, -999, 14, 14 /) REAL(kind=wp), DIMENSION(nlu_dep), PARAMETER :: h = (/ -999, 1, 1, 20, 20, -999, -999, -999, -999, -999, -999, 20, -999, 1 , 1 /) ! Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0: IF (b(lu) .gt. 0.0) then ! Check for u* > 0 (otherwise denominator = 0): IF (ust .gt. 0.0) then rinc = b(lu)*h(lu)*sai/ust ELSE rinc = 1000.0 ENDIF ELSE IF (lu .eq. ilu_grass .or. lu .eq. ilu_other ) then rinc = -999. ! no deposition path for grass, other, and semi-natural ELSE rinc = 0.0 ! no in-canopy resistance ENDIF ENDIF END SUBROUTINE rc_rinc !------------------------------------------------------------------- ! rc_rctot: compute total canopy (or surface) resistance Rc !------------------------------------------------------------------- SUBROUTINE rc_rctot(gstom,gsoil_eff,gw,gc_tot,rc_tot) ! Input/output variables: REAL(kind=wp), INTENT(in) :: gstom ! stomatal conductance (s/m) REAL(kind=wp), INTENT(in) :: gsoil_eff ! effective soil conductance (s/m) REAL(kind=wp), INTENT(in) :: gw ! external leaf conductance (s/m) REAL(kind=wp), INTENT(out) :: gc_tot ! total canopy conductance (m/s) REAL(kind=wp), INTENT(out) :: rc_tot ! total canopy resistance Rc (s/m) ! Total conductance: gc_tot = gstom + gsoil_eff + gw ! Total resistance (note: gw can be negative, but no total emission allowed here): IF (gc_tot .le. 0.0 .or. gw .lt. 0.0) then rc_tot = -9999. ELSE rc_tot = 1./gc_tot ENDIF END SUBROUTINE rc_rctot !------------------------------------------------------------------- ! rc_comp_point_rc_eff: calculate the effective resistance Rc ! based on one or more compensation points !------------------------------------------------------------------- ! ! NH3rc (see depac v3.6 is based on Avero workshop Marc Sutton. p. 173. ! Sutton 1998 AE 473-480) ! ! Documentation by Ferd Sauter, 2008; see also documentation block in header of depac subroutine. ! FS 2009-01-29: variable names made consistent with DEPAC ! FS 2009-03-04: use total compensation point ! ! C: with total compensation point ! D: approximation of C ! ! with classical approach ! zr --------- Catm ! zr --------- Catm ! | ! | ! Ra ! Ra ! | ! | ! Rb ! Rb ! | ! | ! z0 --------- Cc ! z0 --------- Cc ! | ! | ! Rc ! Rc_eff ! | ! | ! --------- Ccomp_tot ! --------- C=0 ! ! ! The effective Rc is defined such that instead of using ! ! F = -vd*[Catm - Ccomp_tot] (1) ! ! we can use the 'normal' flux formula ! ! F = -vd'*Catm, (2) ! ! with vd' = 1/(Ra + Rb + Rc') (3) ! ! and Rc' the effective Rc (rc_eff). ! (Catm - Ccomp_tot) ! vd'*Catm = vd*(Catm - Ccomp_tot) <=> vd' = vd* ------------------ ! Catm ! ! (Catm - Ccomp_tot) ! 1/(Ra + Rb + Rc') = (1/Ra + Rb + Rc) * ------------------ ! Catm ! ! Catm ! (Ra + Rb + Rc') = (Ra + Rb + Rc) * ------------------ ! (Catm - Ccomp_tot) ! ! Catm ! Rc' = (Ra + Rb + Rc) * ------------------ - Ra - Rb ! (Catm - Ccomp_tot) ! ! Catm Catm ! Rc' = (Ra + Rb) [------------------ - 1 ] + Rc * ------------------ ! (Catm - Ccomp_tot) (Catm - Ccomp_tot) ! ! Rc' = [(Ra + Rb)*Ccomp_tot + Rc*Catm ] / (Catm - Ccomp_tot) ! ! This is not the most efficient way to DO this; ! in the current LE version, (1) is used directly in the CALLing routine ! and this routine is not used anymore. ! ! ------------------------------------------------------------------------------------------- ! In fact it is the question IF this correct; note the difference in differential equations ! ! dCatm ! dCatm ! H ----- = F = -vd'*Catm, ! H ----- = F = -vd*[Catm - Ccomp_tot], ! dt ! dt ! ! where in the left colum vd' is a function of Catm and not a constant anymore. ! ! Another problem is that IF Catm -> 0, we end up with an infinite value of the exchange ! velocity vd. ! ------------------------------------------------------------------------------------------- SUBROUTINE rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff) IMPLICIT NONE ! Input/output variables: REAL(kind=wp), INTENT(in) :: ccomp_tot ! total compensation point (weighed average of separate compensation points) (ug/m3) REAL(kind=wp), INTENT(in) :: catm ! atmospheric concentration (ug/m3) REAL(kind=wp), INTENT(in) :: ra ! aerodynamic resistance (s/m) REAL(kind=wp), INTENT(in) :: rb ! boundary layer resistance (s/m) REAL(kind=wp), INTENT(in) :: rc_tot ! total canopy resistance (s/m) REAL(kind=wp), INTENT(out) :: rc_eff ! effective total canopy resistance (s/m) ! Compute effective resistance: IF ( ccomp_tot == 0.0 ) then ! trace with no compensiation point ( or compensation point equal to zero) rc_eff = rc_tot ELSE IF ( ccomp_tot > 0 .and. ( abs( catm - ccomp_tot ) < 1.e-8 ) ) then ! surface concentration (almoast) equal to atmospheric concentration ! no exchange between surface and atmosphere, infinite RC --> vd=0 rc_eff = 9999999999. ELSE IF ( ccomp_tot > 0 ) then ! compensation point available, calculate effective restisance rc_eff = ((ra + rb)*ccomp_tot + rc_tot*catm)/(catm-ccomp_tot) ELSE rc_eff = -999. message_string = 'This should not be possible, check ccomp_tot' CALL message( 'rc_comp_point_rc_eff', 'CM0461', 1, 2, 0, 6, 0 ) ENDIF return END SUBROUTINE rc_comp_point_rc_eff !------------------------------------------------------------------- ! missing: check for data that correspond with a missing deposition path ! this data is represented by -999 !------------------------------------------------------------------- LOGICAL function missing(x) REAL(kind=wp), INTENT(in) :: x ! bandwidth for checking (in)equalities of floats REAL(kind=wp), PARAMETER :: EPS = 1.0e-5 missing = (abs(x + 999.) .le. EPS) END function missing ELEMENTAL FUNCTION sedimentation_velocity(rhopart, partsize, slipcor, visc) RESULT (vs) IMPLICIT NONE ! --- in/out --------------------------------- REAL(kind=wp), INTENT(in) :: rhopart ! kg/m3 REAL(kind=wp), INTENT(in) :: partsize ! m REAL(kind=wp), INTENT(in) :: slipcor ! m REAL(kind=wp), INTENT(in) :: visc ! viscosity REAL(kind=wp) :: vs ! acceleration of gravity: REAL(kind=wp), PARAMETER :: grav = 9.80665 ! m/s2 ! --- begin ---------------------------------- !viscosity & sedimentation velocity vs = rhopart * (partsize**2) * grav * slipcor / (18*visc) END FUNCTION sedimentation_velocity !------------------------------------------------------------------------ ! Boundary-layer deposition resistance following Zhang (2001) !------------------------------------------------------------------------ SUBROUTINE drydepo_aero_zhang_vd( vd, Rs, vs1, partsize, slipcor, & nwet, tsurf, dens1, viscos1, & luc, ftop_lu, ustar_lu) IMPLICIT NONE ! --- in/out --------------------------------- REAL(kind=wp), INTENT(out) :: vd ! deposition velocity (m/s) REAL(kind=wp), INTENT(out) :: Rs ! Sedimentaion resistance (s/m) REAL(kind=wp), INTENT(in) :: vs1 ! sedimentation velocity in lowest layer REAL(kind=wp), INTENT(in) :: partsize ! particle diameter (m) REAL(kind=wp), INTENT(in) :: slipcor ! slip correction factor REAL(kind=wp), INTENT(in) :: tsurf ! surface temperature (K) REAL(kind=wp), INTENT(in) :: dens1 ! air density (kg/m3) in lowest layer REAL(kind=wp), INTENT(in) :: viscos1 ! air viscosity in lowest layer REAL(kind=wp), INTENT(in) :: ftop_lu ! atmospheric resistnace Ra REAL(kind=wp), INTENT(in) :: ustar_lu ! friction velocity u* INTEGER(iwp), INTENT(in) :: nwet ! 1=rain, 9=snowcover INTEGER(iwp), INTENT(in) :: luc ! DEPAC LU ! --- const ---------------------------------- ! acceleration of gravity: REAL(kind=wp), PARAMETER :: grav = 9.80665 ! m/s2 REAL(kind=wp), PARAMETER :: beta = 2.0 REAL(kind=wp), PARAMETER :: epsilon0 = 3.0 REAL(kind=wp), PARAMETER :: kb = 1.38066e-23 REAL(kind=wp), PARAMETER :: pi = 3.141592654_wp !< PI REAL, PARAMETER :: alfa_lu(nlu_dep) = & (/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/) ! grass arabl crops conif decid water urba othr desr ice sav trf wai med sem REAL, PARAMETER :: gamma_lu(nlu_dep) = & (/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/) ! grass arabl crops conif decid water urba othr desr ice sav trf wai med sem REAL, PARAMETER ::A_lu(nlu_dep) = & (/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./) ! grass arabl crops conif decid water urba othr desr ice sav trf wai med sem ! --- local ---------------------------------- REAL(kind=wp) :: kinvisc, diff_part REAL(kind=wp) :: schmidt,stokes, Ebrown, Eimpac, Einterc, Reffic ! --- begin ---------------------------------- ! kinetic viscosity & diffusivity kinvisc = viscos1 / dens1 ! only needed at surface diff_part = kb * tsurf * slipcor / (3*pi*viscos1*partsize) ! Schmidt number schmidt = kinvisc / diff_part ! calculate collection efficiencie E Ebrown = Schmidt**(-gamma_lu(luc)) ! Brownian diffusion ! determine Stokes number, interception efficiency ! and sticking efficiency R (1 = no rebound) IF ( luc == ilu_ice .or. nwet.eq.9 .or. luc == ilu_water_sea .or. luc == ilu_water_inland ) THEN stokes=vs1*ustar_lu**2/(grav*kinvisc) Einterc=0.0 Reffic=1.0 ELSE IF ( luc == ilu_other .or. luc == ilu_desert ) THEN !tundra of desert stokes=vs1*ustar_lu**2/(grav*kinvisc) Einterc=0.0 Reffic=exp(-Stokes**0.5) ELSE stokes=vs1*ustar_lu/(grav*A_lu(luc)*1.e-3) Einterc=0.5*(partsize/(A_lu(luc)*1e-3))**2 Reffic=exp(-Stokes**0.5) END IF ! when surface is wet all particles DO not rebound: IF(nwet.eq.1) Reffic = 1.0 ! determine impaction efficiency: Eimpac = ( stokes / (alfa_lu(luc)+stokes) )**beta ! sedimentation resistance: Rs = 1.0 / ( epsilon0 * ustar_lu * (Ebrown+Eimpac+Einterc) * Reffic ) ! deposition velocity according to Seinfeld and Pandis (= SP, 2006; eq 19.7): ! ! 1 ! vd = ------------------ + vs ! Ra + Rs + Ra*Rs*vs ! ! where: Rs = Rb (in SP, 2006) ! ! vd = 1.0 / ( ftop_lu + Rs + ftop_lu*Rs*vs1) + vs1 END SUBROUTINE drydepo_aero_zhang_vd SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffc, Rb ) ! Compute quasi-laminar boundary layer resistance as a function of landuse and tracer. ! ! Original EMEP formulation by (Simpson et al, 2003) is used. ! IMPLICIT NONE ! --- in/out --------------------------------- LOGICAL , INTENT(in) :: is_water REAL(kind=wp), INTENT(in) :: z0h ! roughness length for heat REAL(kind=wp), INTENT(in) :: ustar ! friction velocity REAL(kind=wp), INTENT(in) :: diffc ! coefficient of diffusivity REAL(kind=wp), INTENT(out) :: Rb ! boundary layer resistance ! --- const ---------------------------------- ! thermal diffusivity of dry air 20 C REAL(kind=wp), PARAMETER :: thk = 0.19e-4 ! http://en.wikipedia.org/wiki/Thermal_diffusivity (@ T=300K) REAL(kind=wp), PARAMETER :: kappa_stab = 0.35 ! von Karman constant ! --- begin --------------------------------- ! Use Simpson et al. (2003) !IF ( is_water ) THEN !Rb = 1.0_wp / (kappa_stab*max(0.01_wp,ustar)) * log(z0h/diffc*kappa_stab*max(0.01_wp,ustar)) !TODO: Check Rb over water calculation!!! ! Rb = 1.0_wp / (kappa_stab*max(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*max(0.1_wp,ustar)) !ELSE Rb = 5.0_wp / max(0.01_wp,ustar) * (thk/diffc)**0.67_wp !END IF END SUBROUTINE get_rb_cell !----------------------------------------------------------------- ! Compute water vapor partial pressure (e_w) ! given specific humidity Q [(kg water)/(kg air)]. ! ! Use that gas law for volume V with temperature T ! holds for the total mixture as well as the water part: ! ! R T / V = p_air / n_air = p_water / n_water ! ! thus: ! ! p_water = p_air n_water / n_air ! ! Use: ! n_air = m_air / xm_air ! [kg air] / [(kg air)/(mole air)] ! and: ! n_water = m_air * Q / xm_water ! [kg water] / [(kg water)/(mole water)] ! thus: ! p_water = p_air Q / (xm_water/xm_air) ! ELEMENTAL FUNCTION watervaporpartialpressure( Q, p ) result ( p_w ) ! --- in/out --------------------------------- REAL(kind=wp), INTENT(in) :: Q ! specific humidity [(kg water)/(kg air)] REAL(kind=wp), INTENT(in) :: p ! air pressure [Pa] REAL(kind=wp) :: p_w ! water vapor partial pressure [Pa] ! --- const ---------------------------------- ! mole mass ratio: REAL(kind=wp), PARAMETER :: eps = xm_H2O / xm_air ! ~ 0.622 ! --- begin ---------------------------------- ! partial pressure of water vapor: p_w = p * Q / eps END function watervaporpartialpressure !------------------------------------------------------------------ ! Saturation vapor pressure. ! From (Stull 1988, eq. 7.5.2d): ! ! e_sat = p0 exp( 17.67 * (T-273.16) / (T-29.66) ) [Pa] ! ! where: ! p0 = 611.2 [Pa] : reference pressure ! ! Arguments: ! T [K] : air temperature ! Result: ! e_sat_w [Pa] : saturation vapor pressure ! ! References: ! Roland B. Stull, 1988 ! An introduction to boundary layer meteorology. ! ELEMENTAL FUNCTION saturationvaporpressure( T ) result( e_sat_w ) ! --- in/out --------------------------------- REAL(kind=wp), INTENT(in) :: T ! temperature [K] REAL(kind=wp) :: e_sat_w ! saturation vapor pressure [Pa] ! --- const ---------------------------------- ! base pressure: REAL(kind=wp), PARAMETER :: p0 = 611.2 ! [Pa] ! --- begin ---------------------------------- ! saturation vapor pressure: e_sat_w = p0 * exp( 17.67 * (T-273.16) / (T-29.66) ) ! [Pa] END FUNCTION saturationvaporpressure !------------------------------------------------------------------------ ! Relative humidity RH [%] is by definition: ! ! e_w # water vapor partial pressure ! Rh = -------- * 100 ! e_sat_w # saturation vapor pressure ! ! Compute here from: ! Q specific humidity [(kg water)/(kg air)] ! T temperature [K] ! P pressure [Pa] ! ELEMENTAL FUNCTION relativehumidity_from_specifichumidity( Q, T, p ) result( Rh ) ! --- in/out --------------------------------- REAL(kind=wp), INTENT(in) :: Q ! specific humidity [(kg water)/(kg air)] REAL(kind=wp), INTENT(in) :: T ! temperature [K] REAL(kind=wp), INTENT(in) :: p ! air pressure [Pa] REAL(kind=wp) :: Rh ! relative humidity [%] ! --- begin ---------------------------------- ! relative humidity: Rh = watervaporpartialpressure( Q, p ) / saturationvaporpressure( T ) * 100.0 END FUNCTION relativehumidity_from_specifichumidity END MODULE chemistry_model_mod