[2657] | 1 | !> @file chemistry_model_mod.f90 |
---|
| 2 | !------------------------------------------------------------------------------! |
---|
[2696] | 3 | ! This file is part of the PALM model system. |
---|
[2657] | 4 | ! |
---|
| 5 | ! PALM is free software: you can redistribute it and/or modify it under the |
---|
| 6 | ! terms of the GNU General Public License as published by the Free Software |
---|
| 7 | ! Foundation, either version 3 of the License, or (at your option) any later |
---|
| 8 | ! version. |
---|
| 9 | ! |
---|
| 10 | ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY |
---|
| 11 | ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR |
---|
| 12 | ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. |
---|
| 13 | ! |
---|
| 14 | ! You should have received a copy of the GNU General Public License along with |
---|
| 15 | ! PALM. If not, see <http://www.gnu.org/licenses/>. |
---|
| 16 | ! |
---|
[2718] | 17 | ! Copyright 1997-2018 Leibniz Universitaet Hannover |
---|
| 18 | ! Copyright 2017-2018 Karlsruhe Institute of Technology |
---|
| 19 | ! Copyright 2017-2018 Freie Universitaet Berlin |
---|
[2657] | 20 | !------------------------------------------------------------------------------! |
---|
| 21 | ! |
---|
| 22 | ! Current revisions: |
---|
| 23 | ! ----------------- |
---|
| 24 | ! |
---|
[2696] | 25 | ! |
---|
[2657] | 26 | ! Former revisions: |
---|
| 27 | ! ----------------- |
---|
[2696] | 28 | ! $Id: chemistry_model_mod.f90 2768 2018-01-24 15:38:29Z raasch $ |
---|
[2768] | 29 | ! Shorten lines to maximum length of 132 characters |
---|
| 30 | ! |
---|
| 31 | ! 2766 2018-01-22 17:17:47Z kanani |
---|
[2766] | 32 | ! Removed preprocessor directive __chem |
---|
| 33 | ! |
---|
| 34 | ! 2756 2018-01-16 18:11:14Z suehring |
---|
[2756] | 35 | ! Fill values in 3D output introduced. |
---|
| 36 | ! |
---|
| 37 | ! 2718 2018-01-02 08:49:38Z maronga |
---|
[2696] | 38 | ! Initial revision |
---|
| 39 | ! |
---|
[2657] | 40 | ! |
---|
[2696] | 41 | ! |
---|
| 42 | ! |
---|
| 43 | ! Authors: |
---|
| 44 | ! -------- |
---|
| 45 | ! @author Renate Forkel |
---|
| 46 | ! @author Farah Kanani-Suehring |
---|
| 47 | ! @author Klaus Ketelsen |
---|
| 48 | ! @author Basit Khan |
---|
| 49 | ! |
---|
| 50 | ! |
---|
[2657] | 51 | !------------------------------------------------------------------------------! |
---|
| 52 | ! Description: |
---|
| 53 | ! ------------ |
---|
| 54 | !> Chemistry model for PALM-4U |
---|
[2668] | 55 | !> @todo Update/clean-up todo list! (FK) |
---|
| 56 | !> @todo Chemistry for prognostic_equations_vector! (FK) |
---|
| 57 | !> @todo Set proper fill values (/= 0) for chem output arrays! (FK) |
---|
| 58 | !> @todo Add routine chem_check_parameters, add checks for inconsistent or |
---|
| 59 | !> unallowed parameter settings (e.g. "only cache optimized version |
---|
| 60 | !> available"). CALL of chem_check_parameters from check_parameters. (FK) |
---|
| 61 | !> @todo Make routine chem_header available, CALL from header.f90 |
---|
| 62 | !> (see e.g. how it is done in routine lsm_header in |
---|
| 63 | !> land_surface_model_mod.f90). chem_header should include all setup |
---|
| 64 | !> info about chemistry parameter settings. (FK) |
---|
[2657] | 65 | !> @todo Implement turbulent inflow of chem spcs in inflow_turbulence. |
---|
| 66 | !> @todo Separate boundary conditions for each chem spcs to be implemented |
---|
| 67 | !> @todo Currently only total concentration are calculated. Resolved, parameterized |
---|
| 68 | !> and chemistry fluxes although partially and some completely coded but |
---|
| 69 | !> are not operational/activated in this version. bK. |
---|
| 70 | !> @todo slight differences in passive scalar and chem spcs when chem reactions |
---|
| 71 | !> turned off. Need to be fixed. bK |
---|
| 72 | !> @todo introduce nesting for chem spcs. bK |
---|
| 73 | !> @todo subroutine set_const_initial_values to be taken out from chemistry_model_mod !bK. |
---|
| 74 | !> @todo chemistry error messages |
---|
| 75 | !> @todo Format this module according to PALM coding standard (see e.g. module |
---|
| 76 | !> template under http://palm.muk.uni-hannover.de/mosaik/downloads/8 or |
---|
| 77 | !> D3_coding_standard.pdf under https://palm.muk.uni-hannover.de/trac/downloads/16) |
---|
| 78 | ! |
---|
| 79 | !------------------------------------------------------------------------------! |
---|
| 80 | |
---|
| 81 | MODULE chemistry_model_mod |
---|
[2766] | 82 | |
---|
[2657] | 83 | USE kinds, ONLY: wp, iwp |
---|
| 84 | USE indices, ONLY: nz, nzb,nzt,nysg,nyng,nxlg,nxrg,nys,nyn |
---|
| 85 | USE pegrid, ONLY: myid, threads_per_task |
---|
| 86 | USE control_parameters, ONLY: dt_3d, ws_scheme_sca, initializing_actions, message_string, & |
---|
| 87 | omega, tsc, intermediate_timestep_count, & |
---|
| 88 | intermediate_timestep_count_max, & |
---|
| 89 | timestep_scheme, use_prescribed_profile_data |
---|
| 90 | USE arrays_3d, ONLY: hyp, pt, rdf_sc, tend, zu |
---|
| 91 | USE chem_gasphase_mod, ONLY: nspec, spc_names, nkppctrl, nmaxfixsteps, & |
---|
| 92 | t_steps, fill_temp, chem_gasphase_integrate, & |
---|
| 93 | nvar, atol, rtol, nphot, phot_names |
---|
| 94 | USE cpulog, ONLY: cpu_log, log_point |
---|
| 95 | |
---|
| 96 | USE chem_modules |
---|
| 97 | |
---|
| 98 | |
---|
| 99 | IMPLICIT NONE |
---|
| 100 | PRIVATE |
---|
| 101 | SAVE |
---|
| 102 | |
---|
| 103 | !- Define chemical variables |
---|
| 104 | |
---|
| 105 | TYPE species_def |
---|
| 106 | CHARACTER(LEN=8) :: name |
---|
| 107 | CHARACTER(LEN=16) :: unit |
---|
| 108 | REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc |
---|
| 109 | REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc_av |
---|
| 110 | REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc_p |
---|
| 111 | REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: tconc_m |
---|
| 112 | REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: cssws_av |
---|
| 113 | REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: flux_s_cs |
---|
| 114 | REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: diss_s_cs |
---|
| 115 | REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:) :: flux_l_cs |
---|
| 116 | REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:) :: diss_l_cs |
---|
| 117 | REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: conc_pr_init |
---|
| 118 | END TYPE species_def |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | TYPE photols_def |
---|
| 122 | CHARACTER(LEN=8) :: name |
---|
| 123 | CHARACTER(LEN=16) :: unit |
---|
| 124 | REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: freq |
---|
| 125 | END TYPE photols_def |
---|
| 126 | |
---|
| 127 | |
---|
| 128 | PUBLIC species_def |
---|
| 129 | PUBLIC photols_def |
---|
| 130 | |
---|
| 131 | |
---|
| 132 | TYPE(species_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC :: chem_species |
---|
| 133 | TYPE(photols_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC :: phot_frequen |
---|
| 134 | |
---|
| 135 | REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_1 |
---|
| 136 | REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_2 |
---|
| 137 | REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_3 |
---|
| 138 | REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_av |
---|
| 139 | REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: freq_1 |
---|
| 140 | |
---|
| 141 | INTEGER,DIMENSION(nkppctrl) :: icntrl ! Fine tuning kpp |
---|
| 142 | REAL(kind=wp),DIMENSION(nkppctrl) :: rcntrl ! Fine tuning kpp |
---|
| 143 | |
---|
| 144 | CHARACTER(10), PUBLIC :: photolysis_scheme |
---|
| 145 | ! 'constant', |
---|
| 146 | ! 'simple' (Simple parameterisation from MCM, Saunders et al., 2003, Atmos. Chem. Phys., 3, 161-180 |
---|
| 147 | ! 'fastj' (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED |
---|
| 148 | |
---|
| 149 | PUBLIC nspec |
---|
| 150 | PUBLIC nvar |
---|
| 151 | PUBLIC spc_names |
---|
| 152 | |
---|
| 153 | !- Interface section |
---|
| 154 | INTERFACE chem_boundary_conds |
---|
| 155 | MODULE PROCEDURE chem_boundary_conds |
---|
| 156 | END INTERFACE chem_boundary_conds |
---|
| 157 | |
---|
| 158 | INTERFACE chem_init |
---|
| 159 | MODULE PROCEDURE chem_init |
---|
| 160 | END INTERFACE chem_init |
---|
| 161 | |
---|
| 162 | INTERFACE chem_init_profiles |
---|
| 163 | MODULE PROCEDURE chem_init_profiles |
---|
[2682] | 164 | END INTERFACE chem_init_profiles |
---|
[2657] | 165 | |
---|
| 166 | INTERFACE chem_parin |
---|
| 167 | MODULE PROCEDURE chem_parin |
---|
| 168 | END INTERFACE chem_parin |
---|
| 169 | |
---|
| 170 | INTERFACE chem_integrate |
---|
| 171 | MODULE PROCEDURE chem_integrate_ij |
---|
| 172 | END INTERFACE chem_integrate |
---|
| 173 | |
---|
| 174 | INTERFACE chem_swap_timelevel |
---|
| 175 | MODULE PROCEDURE chem_swap_timelevel |
---|
| 176 | END INTERFACE chem_swap_timelevel |
---|
| 177 | |
---|
| 178 | INTERFACE chem_define_netcdf_grid |
---|
| 179 | MODULE PROCEDURE chem_define_netcdf_grid |
---|
| 180 | END INTERFACE chem_define_netcdf_grid |
---|
| 181 | |
---|
| 182 | INTERFACE chem_data_output_3d |
---|
| 183 | MODULE PROCEDURE chem_data_output_3d |
---|
| 184 | END INTERFACE chem_data_output_3d |
---|
| 185 | |
---|
| 186 | INTERFACE chem_check_data_output |
---|
| 187 | MODULE PROCEDURE chem_check_data_output |
---|
| 188 | END INTERFACE chem_check_data_output |
---|
| 189 | |
---|
| 190 | INTERFACE chem_check_data_output_pr |
---|
| 191 | MODULE PROCEDURE chem_check_data_output_pr |
---|
| 192 | END INTERFACE chem_check_data_output_pr |
---|
| 193 | |
---|
| 194 | INTERFACE chem_3d_data_averaging |
---|
| 195 | MODULE PROCEDURE chem_3d_data_averaging |
---|
| 196 | END INTERFACE chem_3d_data_averaging |
---|
| 197 | |
---|
| 198 | INTERFACE chem_last_actions |
---|
| 199 | MODULE PROCEDURE chem_last_actions |
---|
| 200 | END INTERFACE chem_last_actions |
---|
| 201 | |
---|
| 202 | INTERFACE chem_read_restart_data |
---|
| 203 | MODULE PROCEDURE chem_read_restart_data |
---|
| 204 | END INTERFACE chem_read_restart_data |
---|
| 205 | |
---|
| 206 | INTERFACE chem_tendency |
---|
| 207 | MODULE PROCEDURE chem_tendency |
---|
| 208 | END INTERFACE chem_tendency |
---|
| 209 | |
---|
| 210 | INTERFACE chem_header |
---|
| 211 | MODULE PROCEDURE chem_header |
---|
| 212 | END INTERFACE chem_header |
---|
| 213 | |
---|
[2668] | 214 | INTERFACE chem_emissions |
---|
| 215 | MODULE PROCEDURE chem_emissions |
---|
| 216 | END INTERFACE chem_emissions |
---|
[2657] | 217 | |
---|
[2668] | 218 | ! INTERFACE chem_write_var_list |
---|
| 219 | ! MODULE PROCEDURE chem_write_var_list |
---|
| 220 | ! END INTERFACE chem_write_var_list |
---|
| 221 | ! |
---|
| 222 | ! INTERFACE chem_read_var_list |
---|
| 223 | ! MODULE PROCEDURE chem_read_var_list |
---|
| 224 | ! END INTERFACE chem_read_var_list |
---|
| 225 | ! |
---|
| 226 | ! INTERFACE chem_skip_var_list |
---|
| 227 | ! MODULE PROCEDURE chem_skip_var_list |
---|
| 228 | ! END INTERFACE chem_skip_var_list |
---|
[2657] | 229 | |
---|
| 230 | |
---|
| 231 | PUBLIC chem_boundary_conds, chem_check_data_output, & |
---|
| 232 | chem_check_data_output_pr, chem_data_output_3d, & |
---|
| 233 | chem_define_netcdf_grid, chem_header, chem_init, & |
---|
| 234 | chem_init_profiles, chem_integrate, chem_parin, & |
---|
| 235 | chem_swap_timelevel, chem_last_actions, chem_read_restart_data, & |
---|
[2668] | 236 | chem_tendency, chem_3d_data_averaging, chem_emissions!, chem_write_var_list, & |
---|
| 237 | !chem_read_var_list, chem_skip_var_list |
---|
[2657] | 238 | |
---|
| 239 | |
---|
| 240 | CONTAINS |
---|
| 241 | |
---|
| 242 | !------------------------------------------------------------------------------! |
---|
| 243 | ! |
---|
| 244 | ! Description: |
---|
| 245 | ! ------------ |
---|
| 246 | !> Subroutine to initialize and set all boundary conditions for chemical species |
---|
| 247 | !------------------------------------------------------------------------------! |
---|
| 248 | |
---|
| 249 | SUBROUTINE chem_boundary_conds( mode ) |
---|
| 250 | |
---|
| 251 | USE control_parameters, & |
---|
| 252 | ONLY: air_chemistry, outflow_l, outflow_n, outflow_r, outflow_s |
---|
| 253 | USE indices, & |
---|
| 254 | ONLY: nxl, nxr, nxlg, nxrg, nyng, nysg, nzt |
---|
| 255 | |
---|
| 256 | ! USE prognostic_equations_mod, & |
---|
| 257 | |
---|
| 258 | USE arrays_3d, & |
---|
| 259 | ONLY: dzu |
---|
| 260 | USE surface_mod, & |
---|
| 261 | ONLY: bc_h |
---|
| 262 | |
---|
| 263 | CHARACTER (len=*), INTENT(IN) :: mode |
---|
| 264 | INTEGER(iwp) :: i !< grid index x direction. |
---|
| 265 | INTEGER(iwp) :: j !< grid index y direction. |
---|
| 266 | INTEGER(iwp) :: k !< grid index z direction. |
---|
| 267 | INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing. |
---|
| 268 | INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls. |
---|
| 269 | INTEGER(iwp) :: m !< running index surface elements. |
---|
| 270 | INTEGER(iwp) :: lsp !< running index for chem spcs. |
---|
| 271 | INTEGER(iwp) :: lph !< running index for photolysis frequencies |
---|
| 272 | |
---|
| 273 | |
---|
| 274 | SELECT CASE ( TRIM( mode ) ) |
---|
| 275 | CASE ( 'init' ) |
---|
| 276 | DO lsp = 1, nspec |
---|
| 277 | IF ( surface_csflux(lsp) == 9999999.9_wp ) THEN |
---|
| 278 | constant_csflux(lsp) = .FALSE. |
---|
| 279 | ENDIF |
---|
| 280 | ENDDO |
---|
| 281 | |
---|
| 282 | IF ( bc_cs_b == 'dirichlet' ) THEN |
---|
| 283 | ibc_cs_b = 0 |
---|
| 284 | ELSEIF ( bc_cs_b == 'neumann' ) THEN |
---|
| 285 | ibc_cs_b = 1 |
---|
| 286 | ELSE |
---|
| 287 | ! message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' ! bK commented |
---|
| 288 | CALL message( 'chem_boundary_conds', 'CHEM001', 1, 2, 0, 6, 0 ) !< chemistry_model_mod should have special error numbers --> "CHEM###", |
---|
| 289 | ENDIF |
---|
| 290 | ! |
---|
| 291 | !-- Set Integer flags and check for possible erroneous settings for top |
---|
| 292 | !-- boundary condition. bK added *_cs_* here. |
---|
| 293 | IF ( bc_cs_t == 'dirichlet' ) THEN |
---|
| 294 | ibc_cs_t = 0 |
---|
| 295 | ELSEIF ( bc_cs_t == 'neumann' ) THEN |
---|
| 296 | ibc_cs_t = 1 |
---|
| 297 | ELSEIF ( bc_cs_t == 'initial_gradient' ) THEN |
---|
| 298 | ibc_cs_t = 2 |
---|
| 299 | ELSEIF ( bc_cs_t == 'nested' ) THEN |
---|
| 300 | ibc_cs_t = 3 |
---|
| 301 | ELSE |
---|
| 302 | ! message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"' |
---|
| 303 | CALL message( 'check_parameters', 'CHEM002', 1, 2, 0, 6, 0 ) |
---|
| 304 | ENDIF |
---|
| 305 | |
---|
| 306 | |
---|
| 307 | CASE ( 'set_bc_bottomtop' ) |
---|
| 308 | !-- Bottom boundary condtions for chemical species |
---|
| 309 | DO lsp = 1, nspec |
---|
| 310 | IF ( ibc_cs_b == 0 ) THEN |
---|
| 311 | DO l = 0, 1 |
---|
| 312 | !-- Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e. |
---|
| 313 | !-- the chem_species(nsp)%conc_p value at the topography top (k-1) |
---|
| 314 | !-- is set; for downward-facing surfaces (l=1), kb=1, i.e. the |
---|
| 315 | !-- value at the topography bottom (k+1) is set. |
---|
| 316 | |
---|
| 317 | kb = MERGE( -1, 1, l == 0 ) |
---|
| 318 | !$OMP PARALLEL DO PRIVATE( i, j, k ) |
---|
| 319 | DO m = 1, bc_h(l)%ns |
---|
| 320 | i = bc_h(l)%i(m) |
---|
| 321 | j = bc_h(l)%j(m) |
---|
| 322 | k = bc_h(l)%k(m) |
---|
| 323 | chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i) |
---|
| 324 | ENDDO |
---|
| 325 | ENDDO |
---|
| 326 | |
---|
| 327 | ELSEIF ( ibc_cs_b == 1 ) THEN |
---|
| 328 | ! in boundary_conds there is som extra loop over m here for passive tracer |
---|
| 329 | DO l = 0, 1 |
---|
| 330 | kb = MERGE( -1, 1, l == 0 ) |
---|
| 331 | !$OMP PARALLEL DO PRIVATE( i, j, k ) |
---|
| 332 | DO m = 1, bc_h(l)%ns |
---|
| 333 | i = bc_h(l)%i(m) |
---|
| 334 | j = bc_h(l)%j(m) |
---|
| 335 | k = bc_h(l)%k(m) |
---|
| 336 | chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc_p(k,j,i) |
---|
| 337 | |
---|
| 338 | !< @todo: chem_species(nsp)%conc_p(k+kb,j,i) = chem_species(nsp)%conc(k+kb,j,i), & |
---|
| 339 | !< pls loop over nsp=1, NSPEC. |
---|
| 340 | !< @todo: We should also think about the possibility to have & |
---|
| 341 | !< individual boundary conditions for each species? This means, bc_cs_b, & |
---|
| 342 | !< bc_cs_t, ibc_cs_b, ibc_cs_t would need to be added to TYPE chem_species(nsp)%, & |
---|
| 343 | !< and the loop over nsp would have to be put outside of this IF-clause. |
---|
| 344 | !< i think its better we keep the same bonundary cond i.e. dirichlet or neumann |
---|
| 345 | !< for all chem spcs. ... !bK |
---|
| 346 | |
---|
| 347 | ENDDO |
---|
| 348 | ENDDO |
---|
| 349 | ENDIF |
---|
| 350 | ENDDO ! end lsp loop |
---|
| 351 | |
---|
| 352 | !-- Top boundary conditions for chemical species - Should this not be done for all species? |
---|
| 353 | IF ( ibc_cs_t == 0 ) THEN |
---|
| 354 | DO lsp = 1, nspec |
---|
| 355 | chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:) |
---|
| 356 | ENDDO |
---|
| 357 | ELSEIF ( ibc_cs_t == 1 ) THEN |
---|
| 358 | DO lsp = 1, nspec |
---|
[2678] | 359 | chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) |
---|
[2657] | 360 | ENDDO |
---|
| 361 | ELSEIF ( ibc_cs_t == 2 ) THEN |
---|
| 362 | DO lsp = 1, nspec |
---|
[2678] | 363 | chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1) |
---|
[2657] | 364 | ENDDO |
---|
| 365 | !@todo: bc_cs_t_val needs to be calculated, |
---|
| 366 | !see bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1) |
---|
| 367 | !(in time_integration). pt_init would be the counterpart to |
---|
| 368 | !chem_species(i)%conc_pr_init (see kchem_driver_FKa1408.f90 of my |
---|
| 369 | !"Hints: prescribing initial concentration. |
---|
| 370 | ENDIF |
---|
| 371 | ! |
---|
| 372 | CASE ( 'set_bc_lateral' ) ! bK commented it |
---|
| 373 | !-- Lateral boundary conditions for chem species at inflow boundary |
---|
| 374 | !-- are automatically set when chem_species concentration is |
---|
| 375 | !-- initialized. The initially set value at the inflow boundary is not |
---|
| 376 | !-- touched during time integration, hence, this boundary value remains |
---|
| 377 | !-- at a constant value, which is the concentration that flows into the |
---|
| 378 | !-- domain. |
---|
| 379 | !-- Lateral boundary conditions for chem species at outflow boundary |
---|
| 380 | |
---|
| 381 | IF ( outflow_s ) THEN |
---|
| 382 | DO lsp = 1, nspec |
---|
| 383 | chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:) |
---|
| 384 | ENDDO |
---|
| 385 | ELSEIF ( outflow_n ) THEN |
---|
| 386 | DO lsp = 1, nspec |
---|
| 387 | chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:) |
---|
| 388 | ENDDO |
---|
| 389 | ELSEIF ( outflow_l ) THEN |
---|
| 390 | DO lsp = 1, nspec |
---|
| 391 | chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl) |
---|
| 392 | ENDDO |
---|
| 393 | ELSEIF ( outflow_r ) THEN |
---|
| 394 | DO lsp = 1, nspec |
---|
| 395 | chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr) |
---|
| 396 | ENDDO |
---|
| 397 | ENDIF |
---|
| 398 | |
---|
| 399 | END SELECT |
---|
| 400 | |
---|
| 401 | END SUBROUTINE chem_boundary_conds |
---|
| 402 | ! |
---|
| 403 | !------------------------------------------------------------------------------! |
---|
| 404 | ! |
---|
| 405 | ! Description: |
---|
| 406 | ! ------------ |
---|
| 407 | !> Subroutine defining initial vertical profiles of chemical species (given by |
---|
| 408 | !> namelist parameters chem_profiles and chem_heights) --> which should work |
---|
| 409 | !> analogue to parameters u_profile, v_profile and uv_heights) |
---|
| 410 | !------------------------------------------------------------------------------! |
---|
| 411 | SUBROUTINE chem_init_profiles !< SUBROUTINE is called from chem_init in case of |
---|
| 412 | !< TRIM( initializing_actions ) /= 'read_restart_data' |
---|
| 413 | !< We still need to see what has to be done in case of restart run |
---|
| 414 | USE chem_modules |
---|
| 415 | |
---|
| 416 | IMPLICIT NONE |
---|
| 417 | |
---|
| 418 | !-- Local variables |
---|
| 419 | INTEGER :: lsp !< running index for number of species in derived data type species_def |
---|
| 420 | INTEGER :: lsp_pr !< running index for number of species in cs_names, cs_profiles etc |
---|
| 421 | INTEGER :: lpr_lev !< running index for profile level for each chem spcs. |
---|
| 422 | INTEGER :: npr_lev !< the next available profile lev |
---|
| 423 | |
---|
| 424 | !----------------- |
---|
| 425 | !-- To prescribe/initialize a vertically constant 'cs_profile', another parameter |
---|
| 426 | !-- "cs_surface" is introduced. If "cs_profile" and "cs_heights" are prescribed, |
---|
| 427 | !-- their values will override the constant profile given by "cs_surface". |
---|
| 428 | |
---|
[2668] | 429 | ! IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN |
---|
[2657] | 430 | lsp_pr = 1 |
---|
| 431 | DO WHILE ( TRIM( cs_name( lsp_pr ) ) /= 'novalue' ) !'novalue' is the default |
---|
| 432 | DO lsp = 1, nspec ! |
---|
| 433 | !-- Create initial profile (conc_pr_init) for each chemical species |
---|
| 434 | IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_pr) ) ) THEN ! |
---|
[2682] | 435 | ! IF ( cs_profile(1,1) == 9999999.9_wp ) THEN |
---|
[2657] | 436 | !-- Set a vertically constant profile based on the surface conc (cs_surface(lsp_pr)) of each species |
---|
| 437 | DO lpr_lev = 0, nzt+1 |
---|
| 438 | chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_pr) |
---|
| 439 | ENDDO |
---|
| 440 | |
---|
[2682] | 441 | ! ELSE |
---|
| 442 | ! IF ( cs_heights(lsp,1) /= 0.0_wp ) THEN |
---|
| 443 | ! message_string = 'cs_heights(1,1) must be 0.0' |
---|
| 444 | ! CALL message( 'chem_check_parameters', 'CHEM0004', 1, 2, 0, 6, 0 ) |
---|
| 445 | ! ENDIF |
---|
| 446 | ! |
---|
| 447 | ! IF ( omega /= 0.0_wp ) THEN |
---|
| 448 | ! message_string = 'Coriolis force must be switched off (by setting omega=0.0)' // & |
---|
| 449 | ! ' when prescribing the forcing by u_profile and v_profile' |
---|
| 450 | ! CALL message( 'check_parameters', 'PA0347', 1, 2, 0, 6, 0 ) |
---|
| 451 | ! ENDIF |
---|
| 452 | ! |
---|
| 453 | ! use_prescribed_profile_data = .TRUE. |
---|
| 454 | ! |
---|
| 455 | ! npr_lev = 1 |
---|
| 456 | ! ! chem_sddpecies(lsp)%conc_pr_init(0) = 0.0_wp |
---|
| 457 | ! DO lpr_lev = 1, nz+1 |
---|
| 458 | ! IF ( npr_lev < 100 ) THEN |
---|
| 459 | ! DO WHILE ( cs_heights(lsp, npr_lev+1) <= zu(lpr_lev) ) |
---|
| 460 | ! npr_lev = npr_lev + 1 |
---|
| 461 | ! IF ( npr_lev == 100 ) EXIT |
---|
| 462 | ! ENDDO |
---|
| 463 | ! ENDIF |
---|
| 464 | ! |
---|
| 465 | ! IF ( npr_lev < 100 .AND. cs_heights(lsp, npr_lev + 1) /= 9999999.9_wp ) THEN |
---|
| 466 | ! chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp, npr_lev) + & |
---|
| 467 | ! ( zu(lpr_lev) - cs_heights(lsp, npr_lev) ) / & |
---|
| 468 | ! ( cs_heights(lsp, npr_lev + 1) - cs_heights(lsp, npr_lev ) ) * & |
---|
| 469 | ! ( cs_profile(lsp, npr_lev + 1) - cs_profile(lsp, npr_lev ) ) |
---|
| 470 | ! ELSE |
---|
| 471 | ! chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp, npr_lev) |
---|
| 472 | ! ENDIF |
---|
| 473 | ! ENDDO |
---|
| 474 | ! ENDIF |
---|
[2657] | 475 | !-- If a profile is prescribed explicity using cs_profiles and cs_heights, we then have to fill the |
---|
| 476 | !-- chem_species(lsp)%conc_pr_init for the specific "lsp" based on the cs_profiles(lsp_pr,:) |
---|
| 477 | !-- and cs_heights(lsp_pr,:). |
---|
| 478 | ENDIF |
---|
| 479 | ENDDO |
---|
| 480 | lsp_pr = lsp_pr + 1 |
---|
| 481 | ENDDO |
---|
[2668] | 482 | ! ELSE |
---|
[2657] | 483 | ! IF (chem_debug1 ) print*,'code to be added for initializing_actions == read_restart_data' !bK |
---|
[2668] | 484 | ! ENDIF |
---|
[2657] | 485 | |
---|
| 486 | !-- Now, go back to chem_init and use the contents of chem_species(lsp)%conc_pr_init to |
---|
| 487 | !-- initialize the 3D conc arrays, as it is currently taken care of in chem_set_constant_values. |
---|
| 488 | !-- After initializing the 3D arrays, these can be used to set the boundary conditions in the |
---|
| 489 | !-- subroutine kchem_boundary_conds, which should be called from boundary_conds.f90. |
---|
| 490 | |
---|
| 491 | |
---|
| 492 | END SUBROUTINE chem_init_profiles |
---|
| 493 | ! |
---|
| 494 | !------------------------------------------------------------------------------! |
---|
| 495 | ! |
---|
| 496 | ! Description: |
---|
| 497 | ! ------------ |
---|
| 498 | !> Subroutine initializating chemistry_model_mod |
---|
| 499 | !------------------------------------------------------------------------------! |
---|
| 500 | SUBROUTINE chem_init |
---|
| 501 | |
---|
| 502 | |
---|
| 503 | USE control_parameters, & |
---|
| 504 | ONLY: message_string, io_blocks, io_group, turbulent_inflow |
---|
| 505 | |
---|
| 506 | USE arrays_3d, & |
---|
| 507 | ONLY: mean_inflow_profiles |
---|
| 508 | |
---|
| 509 | USE pegrid |
---|
| 510 | |
---|
| 511 | IMPLICIT none |
---|
| 512 | !-- local variables |
---|
| 513 | INTEGER :: i,j !< running index for for horiz numerical grid points |
---|
| 514 | INTEGER :: lsp !< running index for chem spcs |
---|
| 515 | INTEGER :: lpr_lev !< running index for chem spcs profile level |
---|
| 516 | ! |
---|
| 517 | !-- NOPOINTER version not implemented yet |
---|
[2682] | 518 | ! #if defined( __nopointer ) |
---|
| 519 | ! message_string = 'The chemistry module only runs with POINTER version' |
---|
| 520 | ! CALL message( 'chemistry_model_mod', 'CM0001', 1, 2, 0, 6, 0 ) |
---|
| 521 | ! #endif |
---|
| 522 | ! |
---|
| 523 | !-- Allocate memory for chemical species |
---|
| 524 | ALLOCATE( chem_species(nspec) ) |
---|
| 525 | ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) |
---|
| 526 | ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) |
---|
| 527 | ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) |
---|
| 528 | ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) |
---|
| 529 | ALLOCATE( phot_frequen(nphot) ) |
---|
| 530 | ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) ) |
---|
| 531 | ALLOCATE( bc_cs_t_val(nspec) ) |
---|
| 532 | ! |
---|
| 533 | !-- Initialize arrays |
---|
| 534 | spec_conc_1 (:,:,:,:) = 0.0_wp |
---|
| 535 | spec_conc_2 (:,:,:,:) = 0.0_wp |
---|
| 536 | spec_conc_3 (:,:,:,:) = 0.0_wp |
---|
| 537 | spec_conc_av(:,:,:,:) = 0.0_wp |
---|
[2657] | 538 | |
---|
[2682] | 539 | |
---|
[2657] | 540 | DO lsp = 1, nspec |
---|
| 541 | chem_species(lsp)%name = spc_names(lsp) |
---|
| 542 | |
---|
| 543 | chem_species(lsp)%conc (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1 (:,:,:,lsp) |
---|
| 544 | chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2 (:,:,:,lsp) |
---|
| 545 | chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_3 (:,:,:,lsp) |
---|
| 546 | chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_av(:,:,:,lsp) |
---|
| 547 | |
---|
| 548 | ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg)) |
---|
[2682] | 549 | chem_species(lsp)%cssws_av = 0.0_wp |
---|
| 550 | ! |
---|
| 551 | !-- (todo (FK): This needs to be revised. This block must go somewhere else) |
---|
[2657] | 552 | ! IF ( ws_scheme_sca ) THEN |
---|
| 553 | ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1)) |
---|
| 554 | ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1)) |
---|
| 555 | ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1)) |
---|
| 556 | ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1)) |
---|
| 557 | chem_species(lsp)%flux_s_cs = 0.0_wp |
---|
| 558 | chem_species(lsp)%flux_l_cs = 0.0_wp |
---|
| 559 | chem_species(lsp)%diss_s_cs = 0.0_wp |
---|
| 560 | chem_species(lsp)%diss_l_cs = 0.0_wp |
---|
| 561 | ! ENDIF |
---|
[2682] | 562 | ! |
---|
| 563 | !-- Allocate memory for initial concentration profiles |
---|
| 564 | !-- (concentration values come from namelist) |
---|
| 565 | !-- (todo (FK): Because of this, chem_init is called in palm before |
---|
| 566 | !-- check_parameters, since conc_pr_init is used there. |
---|
| 567 | !-- We have to find another solution since chem_init should |
---|
| 568 | !-- eventually be called from init_3d_model!!) |
---|
| 569 | ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) ) |
---|
| 570 | chem_species(lsp)%conc_pr_init(:) = 0.0_wp |
---|
[2657] | 571 | |
---|
| 572 | ENDDO |
---|
| 573 | |
---|
| 574 | ! |
---|
| 575 | !-- Set initial concentration of profiles prescribed by parameters cs_profile |
---|
| 576 | !-- and cs_heights in the namelist &chemistry_par |
---|
[2682] | 577 | !-- (todo (FK): chem_init_profiles not ready yet, has some bugs) |
---|
| 578 | ! CALL chem_init_profiles |
---|
[2657] | 579 | ! |
---|
| 580 | !-- Initialize model variables |
---|
| 581 | |
---|
| 582 | |
---|
| 583 | IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & |
---|
| 584 | TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN |
---|
| 585 | |
---|
| 586 | |
---|
| 587 | !-- First model run of a possible job queue. |
---|
| 588 | !-- Initial profiles of the variables must be computes. |
---|
| 589 | IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN |
---|
[2668] | 590 | ! CALL location_message( 'initializing with 1D model profiles', .FALSE. ) |
---|
[2657] | 591 | ! |
---|
| 592 | ! CALL init_1d_model ...... decide to call or not later !bK |
---|
| 593 | |
---|
| 594 | !-- Transfer initial profiles to the arrays of the 3D model |
---|
| 595 | DO lsp = 1, nspec |
---|
| 596 | DO i = nxlg, nxrg |
---|
| 597 | DO j = nysg, nyng |
---|
| 598 | DO lpr_lev = 1, nz + 1 |
---|
| 599 | chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev) |
---|
| 600 | ENDDO |
---|
| 601 | ENDDO |
---|
| 602 | ENDDO |
---|
| 603 | ENDDO |
---|
| 604 | |
---|
| 605 | ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) & |
---|
| 606 | THEN |
---|
[2668] | 607 | ! CALL location_message( 'initializing with constant profiles', .FALSE. ) |
---|
[2657] | 608 | |
---|
| 609 | |
---|
| 610 | |
---|
| 611 | !-- Set initial horizontal velocities at the lowest computational grid |
---|
| 612 | !-- levels to zero in order to avoid too small time steps caused by the |
---|
| 613 | !-- diffusion limit in the initial phase of a run (at k=1, dz/2 occurs |
---|
| 614 | !-- in the limiting formula!). |
---|
| 615 | |
---|
| 616 | DO lsp = 1, nspec |
---|
| 617 | DO i = nxlg, nxrg |
---|
| 618 | DO j = nysg, nyng |
---|
| 619 | chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init !ITS THERE bK |
---|
| 620 | ENDDO |
---|
| 621 | ENDDO |
---|
| 622 | ENDDO |
---|
| 623 | |
---|
[2668] | 624 | ! ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) & |
---|
| 625 | ! THEN |
---|
| 626 | ! CALL location_message( 'initializing by user', .FALSE. ) |
---|
[2657] | 627 | ! |
---|
| 628 | !-- Initialization will completely be done by the user |
---|
[2668] | 629 | !-- (FK: This should be called only once, in init_3d_model, i.e. remove it here) |
---|
| 630 | ! CALL user_init_3d_model |
---|
| 631 | ! CALL location_message( 'finished', .TRUE. ) |
---|
[2657] | 632 | |
---|
| 633 | ENDIF |
---|
| 634 | |
---|
| 635 | !-- Store initial chem spcs profile |
---|
| 636 | ! DO lsp = 1, nvar |
---|
| 637 | ! hom_cs(:,1,115,:) = SPREAD( chem_species(lsp)%conc(:,nys,nxl), 2, statistic_regions+1 ) |
---|
| 638 | ! ENDDO |
---|
| 639 | ! |
---|
| 640 | !-- If required, change the surface chem spcs at the start of the 3D run |
---|
| 641 | IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN |
---|
| 642 | DO lsp = 1, nspec |
---|
| 643 | chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) + & |
---|
| 644 | cs_surface_initial_change(lsp) |
---|
| 645 | ENDDO |
---|
| 646 | ENDIF |
---|
| 647 | ! |
---|
| 648 | !-- Initiale old and new time levels. |
---|
| 649 | DO lsp = 1, nvar |
---|
| 650 | chem_species(lsp)%tconc_m = 0.0_wp |
---|
[2682] | 651 | chem_species(lsp)%conc_p = chem_species(lsp)%conc |
---|
[2657] | 652 | ENDDO |
---|
[2682] | 653 | ! CALL location_message( 'finished', .TRUE. ) |
---|
[2657] | 654 | |
---|
[2682] | 655 | ! |
---|
| 656 | !-- (todo (FK): Restarts not available yet. It is not correct to call |
---|
| 657 | !-- read_3d_binary here) |
---|
[2657] | 658 | ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. & |
---|
| 659 | TRIM( initializing_actions ) == 'cyclic_fill' ) & |
---|
| 660 | THEN |
---|
| 661 | |
---|
[2682] | 662 | message_string = 'so far, restarts are not possible with chemistry' |
---|
| 663 | CALL message( 'chem_init', 'CM0002', 1, 2, 0, 6, 0 ) |
---|
| 664 | ! CALL location_message( 'initializing in case of restart / cyclic_fill', & |
---|
| 665 | ! .FALSE. ) |
---|
| 666 | ! ! |
---|
| 667 | ! !-- When reading data for cyclic fill of 3D prerun data files, read |
---|
| 668 | ! !-- some of the global variables from the restart file which are required |
---|
| 669 | ! !-- for initializing the inflow |
---|
| 670 | ! IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN |
---|
| 671 | ! DO i = 0, io_blocks-1 |
---|
| 672 | ! IF ( i == io_group ) THEN |
---|
| 673 | ! CALL read_parts_of_var_list |
---|
| 674 | ! CALL close_file( 13 ) |
---|
| 675 | ! ENDIF |
---|
| 676 | ! #if defined( __parallel ) |
---|
| 677 | ! CALL MPI_BARRIER( comm2d, ierr ) |
---|
| 678 | ! #endif |
---|
| 679 | ! ENDDO |
---|
| 680 | ! ENDIF |
---|
| 681 | ! ! |
---|
| 682 | ! !-- Read binary data from restart file |
---|
| 683 | ! DO i = 0, io_blocks-1 |
---|
| 684 | ! IF ( i == io_group ) THEN |
---|
| 685 | ! CALL read_3d_binary |
---|
| 686 | ! ENDIF |
---|
| 687 | ! #if defined( __parallel ) |
---|
| 688 | ! CALL MPI_BARRIER( comm2d, ierr ) |
---|
| 689 | ! #endif |
---|
| 690 | ! ENDDO |
---|
| 691 | ! ! |
---|
| 692 | ! !-- Initialization of the turbulence recycling method |
---|
| 693 | ! IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & |
---|
| 694 | ! turbulent_inflow ) THEN |
---|
| 695 | ! ! |
---|
| 696 | ! !-- First store the profiles to be used at the inflow. |
---|
| 697 | ! !-- These profiles are the (temporally) and horizontally averaged vertical |
---|
| 698 | ! !-- profiles from the prerun. Alternatively, prescribed profiles |
---|
| 699 | ! !-- for u,v-components can be used. |
---|
| 700 | ! ALLOCATE( mean_inflow_profiles(nzb:nzt+1,8) ) |
---|
| 701 | ! ! mean_inflow_profiles(:,8) = hom_sum(:,115,0) !cs |
---|
| 702 | ! |
---|
| 703 | ! ! IF ( inflow_l) THEN |
---|
| 704 | ! ! DO j = nysg, nyng |
---|
| 705 | ! ! DO k = nzb, nzt+1 |
---|
| 706 | ! ! DO lsp = 1, nvar |
---|
| 707 | ! ! chem_species(lsp)%conc(k,j,nxlg:-1) = mean_inflow_profiles(k,8) ??? |
---|
| 708 | ! ! ENDDO |
---|
| 709 | ! ! ENDDO |
---|
| 710 | ! ! ENDDO |
---|
| 711 | ! ! ENDIF |
---|
| 712 | ! ENDIF |
---|
| 713 | ! |
---|
| 714 | ! !-- Inside buildings set velocities and TKE back to zero |
---|
| 715 | ! ! IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & |
---|
| 716 | ! ! topography /= 'flat' ) THEN |
---|
| 717 | ! ! ...... |
---|
| 718 | ! ! ENDIF |
---|
| 719 | ! ! |
---|
| 720 | ! !-- Initialize new time levels (only done in order to set boundary values |
---|
| 721 | ! !-- including ghost points) |
---|
| 722 | ! DO lsp = 1, nvar |
---|
| 723 | ! chem_species(lsp)%conc_p = chem_species(lsp)%conc |
---|
| 724 | ! ENDDO |
---|
| 725 | ! ! |
---|
| 726 | ! !-- Allthough tendency arrays are set in prognostic_equations, they have |
---|
| 727 | ! !-- have to be predefined here because they are used (but multiplied with 0) |
---|
| 728 | ! !-- there before they are set. |
---|
| 729 | ! |
---|
| 730 | ! DO lsp = 1, nvar |
---|
| 731 | ! chem_species(lsp)%tconc_m = 0.0_wp |
---|
| 732 | ! ENDDO |
---|
[2657] | 733 | |
---|
| 734 | ELSE |
---|
| 735 | ! |
---|
| 736 | !-- Actually this part of the programm should not be reached |
---|
| 737 | message_string = 'unknown initializing problem' |
---|
[2682] | 738 | CALL message( 'chem_init', 'CM0003', 1, 2, 0, 6, 0 ) |
---|
[2657] | 739 | ENDIF |
---|
| 740 | |
---|
[2682] | 741 | |
---|
| 742 | |
---|
[2657] | 743 | !--- new code add above this line |
---|
| 744 | DO lsp = 1, nphot |
---|
| 745 | phot_frequen(lsp)%name = phot_names(lsp) |
---|
| 746 | ! IF( myid == 0 ) THEN |
---|
| 747 | ! WRITE(6,'(a,i4,3x,a)') 'Photolysis: ',lsp,trim(phot_names(lsp)) |
---|
| 748 | ! ENDIF |
---|
| 749 | phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => freq_1(:,:,:,lsp) |
---|
| 750 | ENDDO |
---|
| 751 | |
---|
| 752 | !-- Set initial values |
---|
| 753 | ! Not required any more, this can now be done with the namelist by setting cs_surface |
---|
| 754 | ! and cs_name without specifying cs_profile (Nevertheless, we still want to keep it for a while) |
---|
| 755 | ! IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN |
---|
| 756 | ! CALL set_const_initial_values |
---|
| 757 | ! ENDIF |
---|
[2682] | 758 | |
---|
[2657] | 759 | RETURN |
---|
| 760 | |
---|
| 761 | CONTAINS |
---|
| 762 | !------------------------------------------------------------------------------! |
---|
| 763 | ! |
---|
| 764 | ! Description: |
---|
| 765 | ! ------------ |
---|
| 766 | !> Subroutine setting constant initial values of chemical species |
---|
| 767 | !------------------------------------------------------------------------------! |
---|
| 768 | SUBROUTINE set_const_initial_values |
---|
| 769 | ! Not required any more, this can now be done with the namelist by setting cs_surface |
---|
| 770 | ! and cs_name without specifying cs_profile (Nevertheless, we still want to keep it for a while) |
---|
| 771 | IMPLICIT none |
---|
| 772 | |
---|
| 773 | !-- local variables |
---|
| 774 | INTEGER :: lsp |
---|
| 775 | |
---|
| 776 | IF(myid == 0) THEN |
---|
| 777 | write(6,'(/,a,/)') ' chemics >>>> Set constant Initial Values: ' |
---|
| 778 | ENDIF |
---|
| 779 | |
---|
| 780 | ! Default values are taken from smog.def from supplied kpp example |
---|
| 781 | DO lsp = 1, nspec |
---|
| 782 | IF(trim(chem_species(lsp)%name) == 'NO') THEN |
---|
| 783 | ! chem_species(lsp)%conc = 8.725*1.0E+08 |
---|
| 784 | ! chem_species(lsp)%conc = 0.05_wp !added by bK |
---|
| 785 | chem_species(lsp)%conc = 0.01_wp !added by RFo |
---|
| 786 | ELSEIF (trim(chem_species(lsp)%name) == 'NO2') THEN |
---|
| 787 | ! chem_species(lsp)%conc = 2.240*1.0E+08 |
---|
| 788 | ! chem_species(lsp)%conc = 0.01_wp !added by bK |
---|
| 789 | chem_species(lsp)%conc = 0.05_wp !added by RFo |
---|
| 790 | ELSEIF( trim( chem_species(lsp)%name ) == 'O3' ) THEN |
---|
| 791 | chem_species(lsp)%conc = 0.05_wp !added by bK |
---|
| 792 | ELSEIF(trim(chem_species(lsp)%name) == 'H2O') THEN |
---|
| 793 | ! chem_species(lsp)%conc = 5.326*1.0E+11 |
---|
| 794 | chem_species(lsp)%conc = 1.30*1.0E+4_wp !added by bK |
---|
| 795 | ELSEIF(trim(chem_species(lsp)%name) == 'O2') THEN |
---|
| 796 | chem_species(lsp)%conc = 2.0*1.0E+5_wp !added by bK |
---|
| 797 | ELSEIF(trim(chem_species(lsp)%name) == 'RH') THEN |
---|
| 798 | chem_species(lsp)%conc = 0.001_wp !added by RFo |
---|
| 799 | ELSEIF(trim(chem_species(lsp)%name) == 'CO') THEN |
---|
| 800 | chem_species(lsp)%conc = 0.5_wp !added by RFo |
---|
| 801 | ELSEIF(trim(chem_species(lsp)%name) == 'RCHO') THEN |
---|
| 802 | ! chem_species(lsp)%conc = 2.0_wp !added by bK |
---|
| 803 | chem_species(lsp)%conc = 0.01_wp !added by RFo |
---|
| 804 | ! ELSEIF(trim(chem_species(lsp)%name) == 'OH') THEN |
---|
| 805 | ! chem_species(lsp)%conc = 1.0*1.0E-07_wp !added by bK |
---|
| 806 | ! ELSEIF(trim(chem_species(lsp)%name) == 'HO2') THEN |
---|
| 807 | ! chem_species(lsp)%conc = 1*1.0E-7_wp !added by bK |
---|
| 808 | ! ELSEIF(trim(chem_species(lsp)%name) == 'RCOO2') THEN ! corrected RFo |
---|
| 809 | ! chem_species(lsp)%conc = 1.0*1.0E-7_wp !added by bK |
---|
| 810 | ! ELSEIF(trim(chem_species(lsp)%name) == 'RCOO2NO2') THEN |
---|
| 811 | ! chem_species(lsp)%conc = 1.0*1.0E-7_wp !added by bK |
---|
| 812 | ELSE |
---|
| 813 | ! H2O = 2.0e+04; |
---|
| 814 | chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) = 0.0_wp |
---|
| 815 | ENDIF |
---|
| 816 | |
---|
| 817 | IF(myid == 0) THEN |
---|
[2768] | 818 | WRITE(6,'(a,3x,a,3x,a,e12.4)') ' Species: ',chem_species(lsp)%name(1:7),' & |
---|
| 819 | Initial Value = ',chem_species(lsp)%conc(nzb,nysg,nxlg) |
---|
[2657] | 820 | ENDIF |
---|
| 821 | ENDDO |
---|
| 822 | |
---|
| 823 | ! #if defined( __nopointer ) |
---|
| 824 | ! !kk Hier mit message abbrechen |
---|
| 825 | ! if(myid == 0) then |
---|
| 826 | ! write(6,*) ' KPP does only run with POINTER Version' |
---|
| 827 | ! end if |
---|
| 828 | ! stop 'error' |
---|
| 829 | ! #endif |
---|
| 830 | |
---|
| 831 | RETURN |
---|
| 832 | END SUBROUTINE set_const_initial_values |
---|
| 833 | |
---|
| 834 | |
---|
| 835 | END SUBROUTINE chem_init |
---|
| 836 | ! |
---|
| 837 | !------------------------------------------------------------------------------! |
---|
| 838 | ! |
---|
| 839 | ! Description: |
---|
| 840 | ! ------------ |
---|
| 841 | !> Subroutine defining parin for &chemistry_par for chemistry model |
---|
| 842 | !------------------------------------------------------------------------------! |
---|
| 843 | SUBROUTINE chem_parin |
---|
[2682] | 844 | |
---|
[2657] | 845 | USE control_parameters, & |
---|
| 846 | ONLY: air_chemistry |
---|
[2682] | 847 | |
---|
[2657] | 848 | USE chem_modules |
---|
[2682] | 849 | |
---|
| 850 | USE kinds |
---|
[2657] | 851 | |
---|
| 852 | IMPLICIT none |
---|
| 853 | |
---|
[2682] | 854 | CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file |
---|
| 855 | |
---|
| 856 | REAL(wp), DIMENSION(nmaxfixsteps) :: my_steps !< List of fixed timesteps my_step(1) = 0.0 automatic stepping |
---|
[2657] | 857 | |
---|
[2668] | 858 | NAMELIST /chemistry_par/ bc_cs_b, & |
---|
| 859 | bc_cs_t, & |
---|
| 860 | call_chem_at_all_substeps, & |
---|
| 861 | chem_debug0, & |
---|
| 862 | chem_debug1, & |
---|
| 863 | chem_debug2, & |
---|
| 864 | chem_gasphase_on, & |
---|
| 865 | cs_heights, & |
---|
| 866 | cs_name, & |
---|
| 867 | cs_profile, & |
---|
| 868 | cs_profile_name, & |
---|
| 869 | cs_surface, & |
---|
| 870 | emiss_factor_main, & |
---|
| 871 | emiss_factor_side, & |
---|
| 872 | icntrl, & |
---|
| 873 | main_street_id, & |
---|
[2678] | 874 | max_street_id, & |
---|
[2668] | 875 | my_steps, & |
---|
| 876 | rcntrl, & |
---|
| 877 | side_street_id, & |
---|
| 878 | photolysis_scheme, & |
---|
| 879 | wall_csflux, & |
---|
| 880 | cs_vertical_gradient, & |
---|
| 881 | top_csflux, & |
---|
| 882 | surface_csflux, & |
---|
[2678] | 883 | surface_csflux_name, & |
---|
[2668] | 884 | cs_surface_initial_change, & |
---|
| 885 | cs_vertical_gradient_level |
---|
[2657] | 886 | |
---|
| 887 | !-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj) |
---|
| 888 | !-- so this way we could prescribe a specific flux value for each species |
---|
| 889 | !> chemistry_par for initial profiles |
---|
| 890 | !> cs_names = 'O3', 'NO2', 'NO', ... to set initial profiles) |
---|
| 891 | !> cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3) |
---|
| 892 | !> cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.) |
---|
| 893 | !> cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, ..... (chem spcs conc at height lvls chem_heights(1,:)) etc. |
---|
| 894 | !> If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)" |
---|
| 895 | !> then write these cs_surface values to chem_species(lsp)%conc_pr_init(:) |
---|
| 896 | |
---|
[2668] | 897 | ! |
---|
| 898 | !-- Read chem namelist |
---|
| 899 | !-- (todo: initialize these parameters in declaration part, do this for |
---|
| 900 | !-- all chemistry_par namelist parameters) |
---|
[2657] | 901 | icntrl = 0 |
---|
| 902 | rcntrl = 0.0_wp |
---|
| 903 | my_steps = 0.0_wp |
---|
| 904 | icntrl(2) = 1 |
---|
| 905 | photolysis_scheme = 'simple' |
---|
| 906 | atol = 1.0_wp |
---|
| 907 | rtol = 0.01_wp |
---|
| 908 | ! |
---|
| 909 | !-- Try to find chemistry package |
---|
| 910 | REWIND ( 11 ) |
---|
| 911 | line = ' ' |
---|
| 912 | DO WHILE ( INDEX( line, '&chemistry_par' ) == 0 ) |
---|
[2668] | 913 | READ ( 11, '(A)', END=10 ) line |
---|
[2657] | 914 | ENDDO |
---|
| 915 | BACKSPACE ( 11 ) |
---|
| 916 | ! |
---|
[2668] | 917 | !-- Read chemistry namelist |
---|
[2657] | 918 | READ ( 11, chemistry_par ) |
---|
[2668] | 919 | ! |
---|
| 920 | !-- Enable chemistry model |
---|
| 921 | air_chemistry = .TRUE. |
---|
[2657] | 922 | |
---|
[2668] | 923 | ! |
---|
| 924 | !-- FK: Not in use/needed (yet), since chemistry_par namelist is always read in |
---|
| 925 | ! IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN |
---|
| 926 | ! CALL chem_read_var_list |
---|
| 927 | ! ENDIF |
---|
| 928 | |
---|
[2657] | 929 | 10 CONTINUE |
---|
| 930 | |
---|
[2668] | 931 | t_steps = my_steps !(todo: Why not directly make t_steps a |
---|
| 932 | ! namelist parameter?) |
---|
[2657] | 933 | |
---|
| 934 | END SUBROUTINE chem_parin |
---|
| 935 | |
---|
| 936 | ! |
---|
| 937 | !------------------------------------------------------------------------------! |
---|
| 938 | ! |
---|
| 939 | ! Description: |
---|
| 940 | ! ------------ |
---|
| 941 | !> Subroutine to integrate chemical species in the given chemical mechanism |
---|
| 942 | !------------------------------------------------------------------------------! |
---|
| 943 | |
---|
| 944 | SUBROUTINE chem_integrate_ij (i, j) |
---|
| 945 | |
---|
| 946 | USE cpulog, & |
---|
| 947 | ONLY: cpu_log, log_point |
---|
| 948 | USE statistics, & ! ## RFo |
---|
| 949 | ONLY: weight_pres |
---|
| 950 | USE control_parameters, & ! ## RFo |
---|
| 951 | ONLY: dt_3d, intermediate_timestep_count |
---|
| 952 | |
---|
| 953 | IMPLICIT none |
---|
| 954 | INTEGER,INTENT(IN) :: i,j |
---|
| 955 | |
---|
| 956 | !-- local variables |
---|
| 957 | INTEGER(iwp) :: lsp !< running index for chem spcs. |
---|
| 958 | INTEGER(iwp) :: lph !< running index for photolysis frequencies |
---|
| 959 | INTEGER :: k,m,istatf |
---|
| 960 | INTEGER,dimension(20) :: istatus |
---|
| 961 | REAL(kind=wp),dimension(nzb+1:nzt,nspec) :: tmp_conc |
---|
| 962 | REAL(kind=wp),dimension(nzb+1:nzt) :: tmp_temp |
---|
| 963 | REAL(kind=wp),dimension(nzb+1:nzt,nphot) :: tmp_phot |
---|
| 964 | REAL(kind=wp),dimension(nzb+1:nzt) :: tmp_fact |
---|
| 965 | REAL(kind=wp),dimension(nzb+1:nzt) :: tmp_fact_i !< conversion factor between molecules cm^{-3} and ppm |
---|
| 966 | REAL(wp) :: conv !< conversion factor |
---|
| 967 | REAL(wp), PARAMETER :: ppm2fr = 1.0e-6_wp !< Conversion factor ppm to fraction |
---|
| 968 | REAL(wp), PARAMETER :: pref_i = 1._wp / 100000.0_wp !< inverse reference pressure (1/Pa) |
---|
| 969 | REAL(wp), PARAMETER :: t_std = 273.15_wp !< standard pressure (Pa) |
---|
| 970 | REAL(wp), PARAMETER :: p_std = 101325.0_wp !< standard pressure (Pa) |
---|
| 971 | REAL(wp), PARAMETER :: r_cp = 0.286_wp !< R / cp (exponent for potential temperature) |
---|
| 972 | REAL(wp), PARAMETER :: vmolcm = 22.414e3_wp !< Mole volume (22.414 l) in cm^{-3} |
---|
| 973 | REAL(wp), PARAMETER :: xna = 6.022e23_wp !< Avogadro number (molecules/mol) |
---|
| 974 | |
---|
| 975 | |
---|
| 976 | REAL(kind=wp) :: dt_chem |
---|
| 977 | |
---|
| 978 | CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'start' ) |
---|
| 979 | !< Set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry |
---|
| 980 | IF (chem_gasphase_on) THEN |
---|
| 981 | |
---|
| 982 | tmp_temp(:) = pt(:,j,i) * ( hyp(nzb+1:nzt) / 100000.0_wp )**0.286_wp |
---|
| 983 | ! ppm to molecules/cm**3 |
---|
| 984 | ! 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 ) |
---|
| 985 | conv = ppm2fr * xna / vmolcm |
---|
| 986 | tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std) |
---|
| 987 | tmp_fact_i = 1.0_wp/tmp_fact |
---|
| 988 | |
---|
| 989 | CALL fill_temp (istatf, tmp_temp) !< Load constant temperature into kpp context |
---|
| 990 | ! CALL fill_temp (istatf, pt(nzb+1:nzt,j,i)) !< Load temperature into kpp context |
---|
| 991 | |
---|
| 992 | DO lsp = 1,nspec |
---|
| 993 | tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:) |
---|
| 994 | ENDDO |
---|
| 995 | |
---|
| 996 | DO lph = 1,nphot |
---|
| 997 | tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i) |
---|
| 998 | ENDDO |
---|
| 999 | |
---|
| 1000 | IF(myid == 0 .AND. chem_debug0 ) THEN |
---|
| 1001 | IF (i == 10 .and. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d |
---|
| 1002 | ENDIF |
---|
| 1003 | |
---|
| 1004 | !-- Compute length of time step # RFo |
---|
| 1005 | IF ( call_chem_at_all_substeps ) THEN |
---|
| 1006 | dt_chem = dt_3d * weight_pres(intermediate_timestep_count) |
---|
| 1007 | ELSE |
---|
| 1008 | dt_chem = dt_3d |
---|
| 1009 | ENDIF |
---|
| 1010 | |
---|
| 1011 | CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'start' ) |
---|
| 1012 | |
---|
| 1013 | |
---|
| 1014 | CALL chem_gasphase_integrate (dt_chem, tmp_conc, tmp_temp, tmp_phot, istatus=istatus) |
---|
| 1015 | |
---|
| 1016 | |
---|
| 1017 | CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'stop' ) |
---|
| 1018 | |
---|
| 1019 | DO lsp = 1,nspec |
---|
| 1020 | chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:) ! RFo |
---|
| 1021 | ENDDO |
---|
| 1022 | |
---|
| 1023 | ! IF (myid == 0 .AND. chem_debug2 ) THEN |
---|
| 1024 | ! IF (i == 10 .and. j == 10) WRITE(6,'(a,8i7)') ' KPP Status ',istatus(1:8) |
---|
| 1025 | ! write(6,'(a,8i7)') ' KPP Status ',istatus(1:8) |
---|
| 1026 | ! ENDIF |
---|
| 1027 | |
---|
| 1028 | ENDIF |
---|
| 1029 | CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'stop' ) |
---|
| 1030 | |
---|
| 1031 | RETURN |
---|
| 1032 | END SUBROUTINE chem_integrate_ij |
---|
| 1033 | ! |
---|
| 1034 | !------------------------------------------------------------------------------! |
---|
| 1035 | ! |
---|
| 1036 | ! Description: |
---|
| 1037 | ! ------------ |
---|
| 1038 | !> Subroutine for swapping of timelevels for chemical species |
---|
| 1039 | !> called out from subroutine swap_timelevel |
---|
| 1040 | !------------------------------------------------------------------------------! |
---|
| 1041 | |
---|
| 1042 | SUBROUTINE chem_swap_timelevel (level) |
---|
| 1043 | IMPLICIT none |
---|
| 1044 | |
---|
| 1045 | INTEGER,INTENT(IN) :: level |
---|
| 1046 | |
---|
| 1047 | !-- local variables |
---|
| 1048 | |
---|
| 1049 | INTEGER :: lsp |
---|
| 1050 | |
---|
| 1051 | ! print*,' *** entering chem_swap_timelevel ***) ' |
---|
| 1052 | if(level == 0) then |
---|
| 1053 | do lsp=1, nvar |
---|
| 1054 | chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp) |
---|
| 1055 | chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp) |
---|
| 1056 | end do |
---|
| 1057 | else |
---|
| 1058 | do lsp=1, nvar |
---|
| 1059 | chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp) |
---|
| 1060 | chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp) |
---|
| 1061 | end do |
---|
| 1062 | end if |
---|
| 1063 | |
---|
| 1064 | RETURN |
---|
| 1065 | END SUBROUTINE chem_swap_timelevel |
---|
| 1066 | |
---|
| 1067 | !------------------------------------------------------------------------------! |
---|
| 1068 | ! |
---|
| 1069 | ! Description: |
---|
| 1070 | ! ------------ |
---|
| 1071 | !> Subroutine defining appropriate grid for netcdf variables. |
---|
| 1072 | !> It is called out from subroutine netcdf. |
---|
| 1073 | !------------------------------------------------------------------------------! |
---|
| 1074 | SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) |
---|
| 1075 | |
---|
| 1076 | IMPLICIT NONE |
---|
| 1077 | |
---|
| 1078 | CHARACTER (LEN=*), INTENT(IN) :: var !< |
---|
| 1079 | LOGICAL, INTENT(OUT) :: found !< |
---|
| 1080 | CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< |
---|
| 1081 | CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< |
---|
| 1082 | CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< |
---|
| 1083 | |
---|
| 1084 | found = .TRUE. |
---|
| 1085 | |
---|
| 1086 | if(var(1:3) == 'kc_') then !< always the same grid for chemistry variables |
---|
| 1087 | grid_x = 'x' |
---|
| 1088 | grid_y = 'y' |
---|
| 1089 | grid_z = 'zu' !< kk Use same z axis as u variables. Has to be checked if OK |
---|
| 1090 | else |
---|
| 1091 | found = .FALSE. |
---|
| 1092 | grid_x = 'none' |
---|
| 1093 | grid_y = 'none' |
---|
| 1094 | grid_z = 'none' |
---|
| 1095 | end if |
---|
| 1096 | |
---|
| 1097 | ! write(6,*) 'chem_define_netcdf_grid ',TRIM(var),' ',trim(grid_x),' ',found |
---|
| 1098 | |
---|
| 1099 | END SUBROUTINE chem_define_netcdf_grid |
---|
| 1100 | ! |
---|
| 1101 | !------------------------------------------------------------------------------! |
---|
| 1102 | ! |
---|
| 1103 | ! Description: |
---|
| 1104 | ! ------------ |
---|
| 1105 | !> Subroutine for checking data output for chemical species |
---|
| 1106 | !------------------------------------------------------------------------------! |
---|
| 1107 | |
---|
| 1108 | SUBROUTINE chem_check_data_output( var, unit, i, ilen, k ) |
---|
| 1109 | |
---|
| 1110 | |
---|
| 1111 | USE control_parameters, & |
---|
| 1112 | ONLY: data_output, message_string |
---|
| 1113 | |
---|
| 1114 | IMPLICIT NONE |
---|
| 1115 | |
---|
| 1116 | CHARACTER (LEN=*) :: unit !< |
---|
| 1117 | CHARACTER (LEN=*) :: var !< |
---|
| 1118 | |
---|
| 1119 | INTEGER(iwp) :: i, lsp |
---|
| 1120 | INTEGER(iwp) :: ilen |
---|
| 1121 | INTEGER(iwp) :: k |
---|
| 1122 | |
---|
| 1123 | CHARACTER(len=16) :: spec_name |
---|
| 1124 | |
---|
| 1125 | unit = 'illegal' |
---|
| 1126 | |
---|
| 1127 | spec_name = TRIM(var(4:)) !< var 1:3 is 'kc_'. |
---|
| 1128 | |
---|
| 1129 | DO lsp=1,nspec |
---|
| 1130 | IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN |
---|
| 1131 | unit = 'ppm' |
---|
| 1132 | ENDIF |
---|
| 1133 | ! It is possible to plant PM10 and PM25 into the gasphase chemistry code |
---|
| 1134 | ! ! as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): |
---|
| 1135 | ! ! set unit to micrograms per m**3 for PM10 and PM25 (PM2.5) |
---|
| 1136 | IF (spec_name(1:2) == 'PM') THEN |
---|
| 1137 | unit = 'ug m-3' |
---|
| 1138 | ENDIF |
---|
| 1139 | ENDDO |
---|
| 1140 | |
---|
| 1141 | DO lsp=1,nphot |
---|
| 1142 | IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name)) THEN |
---|
| 1143 | unit = 'sec-1' |
---|
| 1144 | ENDIF |
---|
| 1145 | ENDDO |
---|
| 1146 | |
---|
| 1147 | |
---|
| 1148 | RETURN |
---|
| 1149 | END SUBROUTINE chem_check_data_output |
---|
| 1150 | ! |
---|
| 1151 | !------------------------------------------------------------------------------! |
---|
| 1152 | ! |
---|
| 1153 | ! Description: |
---|
| 1154 | ! ------------ |
---|
| 1155 | !> Subroutine for checking data output of profiles for chemistry model |
---|
| 1156 | !------------------------------------------------------------------------------! |
---|
| 1157 | |
---|
| 1158 | SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit ) |
---|
| 1159 | |
---|
| 1160 | |
---|
| 1161 | USE arrays_3d, & |
---|
| 1162 | ONLY: zu |
---|
| 1163 | |
---|
| 1164 | USE control_parameters, & |
---|
| 1165 | ONLY: data_output_pr, message_string, air_chemistry |
---|
| 1166 | |
---|
| 1167 | USE indices |
---|
| 1168 | |
---|
| 1169 | USE profil_parameter |
---|
| 1170 | |
---|
| 1171 | USE statistics |
---|
| 1172 | |
---|
| 1173 | |
---|
| 1174 | IMPLICIT NONE |
---|
| 1175 | |
---|
| 1176 | CHARACTER (LEN=*) :: unit !< |
---|
| 1177 | CHARACTER (LEN=*) :: variable !< |
---|
| 1178 | CHARACTER (LEN=*) :: dopr_unit |
---|
| 1179 | CHARACTER(len=16) :: spec_name |
---|
| 1180 | |
---|
| 1181 | INTEGER(iwp) :: var_count, lsp !< |
---|
| 1182 | |
---|
| 1183 | |
---|
| 1184 | spec_name = TRIM(variable(4:)) |
---|
| 1185 | ! write(9,*) 'fm #32 .. air_chemistry ', air_chemistry |
---|
| 1186 | |
---|
| 1187 | IF ( .NOT. air_chemistry ) THEN |
---|
| 1188 | message_string = 'data_output_pr = ' // & |
---|
| 1189 | TRIM( data_output_pr(var_count) ) // ' is not imp' // & |
---|
| 1190 | 'lemented for air_chemistry = .FALSE.' |
---|
| 1191 | ! CALL message( 'check_parameters', 'PA0185', 1, 2, 0, 6, 0 ) |
---|
| 1192 | ELSE |
---|
| 1193 | DO lsp = 1, nspec |
---|
| 1194 | IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN |
---|
| 1195 | dopr_index(var_count) = 900 |
---|
| 1196 | dopr_unit = 'ppm' |
---|
| 1197 | hom(:,2,900,:) = SPREAD( zu, 2, statistic_regions+1 ) |
---|
| 1198 | ENDIF |
---|
| 1199 | ENDDO |
---|
| 1200 | ENDIF |
---|
| 1201 | |
---|
| 1202 | END SUBROUTINE chem_check_data_output_pr |
---|
| 1203 | ! |
---|
| 1204 | !------------------------------------------------------------------------------! |
---|
| 1205 | ! |
---|
| 1206 | ! Description: |
---|
| 1207 | ! ------------ |
---|
| 1208 | !> Subroutine defining 3D output variables for chemical species |
---|
| 1209 | !------------------------------------------------------------------------------! |
---|
| 1210 | |
---|
| 1211 | |
---|
[2756] | 1212 | SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value ) |
---|
[2657] | 1213 | |
---|
| 1214 | |
---|
| 1215 | USE indices |
---|
| 1216 | |
---|
| 1217 | USE kinds |
---|
| 1218 | |
---|
| 1219 | |
---|
| 1220 | IMPLICIT NONE |
---|
| 1221 | |
---|
| 1222 | CHARACTER (LEN=*) :: variable !< |
---|
| 1223 | LOGICAL :: found !< |
---|
| 1224 | INTEGER(iwp) :: av !< |
---|
| 1225 | |
---|
[2756] | 1226 | REAL(wp) :: fill_value !< |
---|
| 1227 | REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) :: local_pf |
---|
| 1228 | |
---|
| 1229 | |
---|
[2657] | 1230 | !-- local variables |
---|
| 1231 | |
---|
| 1232 | INTEGER :: i, j, k, lsp |
---|
| 1233 | CHARACTER(len=16) :: spec_name |
---|
| 1234 | |
---|
| 1235 | |
---|
| 1236 | found = .FALSE. |
---|
| 1237 | |
---|
| 1238 | spec_name = TRIM(variable(4:)) |
---|
| 1239 | !av == 0 |
---|
| 1240 | |
---|
| 1241 | DO lsp=1,nspec |
---|
| 1242 | IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN |
---|
[2768] | 1243 | IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ',TRIM(variable),' & |
---|
| 1244 | ',TRIM(chem_species(lsp)%name) |
---|
[2657] | 1245 | |
---|
| 1246 | IF (av == 0) THEN |
---|
| 1247 | DO i = nxl, nxr |
---|
| 1248 | DO j = nys, nyn |
---|
| 1249 | DO k = nzb, nzt+1 |
---|
[2756] | 1250 | local_pf(i,j,k) = MERGE( & |
---|
| 1251 | chem_species(lsp)%conc(k,j,i), & |
---|
| 1252 | REAL( fill_value, KIND = wp ), & |
---|
| 1253 | BTEST( wall_flags_0(k,j,i), 0 ) ) |
---|
[2657] | 1254 | ENDDO |
---|
| 1255 | ENDDO |
---|
| 1256 | ENDDO |
---|
| 1257 | |
---|
| 1258 | ELSE |
---|
| 1259 | DO i = nxl, nxr |
---|
| 1260 | DO j = nys, nyn |
---|
| 1261 | DO k = nzb, nzt+1 |
---|
[2756] | 1262 | local_pf(i,j,k) = MERGE( & |
---|
| 1263 | chem_species(lsp)%conc_av(k,j,i),& |
---|
| 1264 | REAL( fill_value, KIND = wp ), & |
---|
| 1265 | BTEST( wall_flags_0(k,j,i), 0 ) ) |
---|
[2657] | 1266 | ENDDO |
---|
| 1267 | ENDDO |
---|
| 1268 | ENDDO |
---|
| 1269 | ENDIF |
---|
| 1270 | |
---|
| 1271 | found = .TRUE. |
---|
| 1272 | ENDIF |
---|
| 1273 | ENDDO |
---|
| 1274 | |
---|
| 1275 | RETURN |
---|
| 1276 | END SUBROUTINE chem_data_output_3d |
---|
| 1277 | ! |
---|
| 1278 | !------------------------------------------------------------------------------! |
---|
| 1279 | ! |
---|
| 1280 | ! Description: |
---|
| 1281 | ! ------------ |
---|
| 1282 | !> Subroutine for averaging 3D data of chemical species |
---|
| 1283 | !------------------------------------------------------------------------------! |
---|
| 1284 | |
---|
| 1285 | SUBROUTINE chem_3d_data_averaging ( mode, variable ) |
---|
| 1286 | |
---|
| 1287 | USE control_parameters |
---|
| 1288 | |
---|
| 1289 | USE indices |
---|
| 1290 | |
---|
| 1291 | USE kinds |
---|
| 1292 | |
---|
| 1293 | USE surface_mod, & |
---|
| 1294 | ONLY: surf_def_h, surf_lsm_h, surf_usm_h |
---|
| 1295 | |
---|
| 1296 | IMPLICIT NONE |
---|
| 1297 | |
---|
| 1298 | CHARACTER (LEN=*) :: mode !< |
---|
| 1299 | CHARACTER (LEN=*) :: variable !< |
---|
| 1300 | |
---|
| 1301 | |
---|
| 1302 | INTEGER(iwp) :: i !< grid index x direction |
---|
| 1303 | INTEGER(iwp) :: j !< grid index y direction |
---|
| 1304 | INTEGER(iwp) :: k !< grid index z direction |
---|
| 1305 | INTEGER(iwp) :: m !< running index surface type |
---|
| 1306 | INTEGER(iwp) :: lsp !< running index for chem spcs |
---|
| 1307 | INTEGER(iwp) :: lsp_2 !< it looks like redundent .. will be delted ..bK |
---|
| 1308 | |
---|
| 1309 | IF ( mode == 'allocate' ) THEN |
---|
| 1310 | DO lsp = 1, nspec |
---|
| 1311 | IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN |
---|
| 1312 | ! lsp_2 = lsp |
---|
| 1313 | chem_species(lsp)%conc_av = 0.0_wp |
---|
| 1314 | |
---|
| 1315 | ENDIF |
---|
| 1316 | ENDDO |
---|
| 1317 | |
---|
| 1318 | ELSEIF ( mode == 'sum' ) THEN |
---|
| 1319 | |
---|
| 1320 | DO lsp = 1, nspec |
---|
| 1321 | IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN |
---|
| 1322 | ! lsp_2 = lsp |
---|
| 1323 | DO i = nxlg, nxrg |
---|
| 1324 | DO j = nysg, nyng |
---|
| 1325 | DO k = nzb, nzt+1 |
---|
| 1326 | chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) + & |
---|
| 1327 | chem_species(lsp)%conc(k,j,i) |
---|
| 1328 | ENDDO |
---|
| 1329 | ENDDO |
---|
| 1330 | ENDDO |
---|
| 1331 | ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') ) THEN |
---|
| 1332 | DO m = 1, surf_def_h(0)%ns |
---|
| 1333 | i = surf_def_h(0)%i(m) |
---|
| 1334 | j = surf_def_h(0)%j(m) |
---|
| 1335 | chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) + surf_def_h(0)%cssws(lsp,m) |
---|
| 1336 | ENDDO |
---|
| 1337 | DO m = 1, surf_lsm_h%ns |
---|
| 1338 | i = surf_lsm_h%i(m) |
---|
| 1339 | j = surf_lsm_h%j(m) |
---|
| 1340 | chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) + surf_lsm_h%cssws(lsp,m) |
---|
| 1341 | ENDDO |
---|
| 1342 | DO m = 1, surf_usm_h%ns |
---|
| 1343 | i = surf_usm_h%i(m) |
---|
| 1344 | j = surf_usm_h%j(m) |
---|
| 1345 | chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) + surf_usm_h%cssws(lsp,m) |
---|
| 1346 | ENDDO |
---|
| 1347 | |
---|
| 1348 | ENDIF |
---|
| 1349 | ENDDO |
---|
| 1350 | |
---|
| 1351 | ELSEIF ( mode == 'average' ) THEN |
---|
| 1352 | DO lsp = 1, nspec |
---|
| 1353 | IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN |
---|
| 1354 | ! lsp_2 = lsp |
---|
| 1355 | DO i = nxlg, nxrg |
---|
| 1356 | DO j = nysg, nyng |
---|
| 1357 | DO k = nzb, nzt+1 |
---|
| 1358 | chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) / REAL( average_count_3d, KIND=wp ) |
---|
| 1359 | ENDDO |
---|
| 1360 | ENDDO |
---|
| 1361 | ENDDO |
---|
| 1362 | |
---|
| 1363 | ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') ) THEN |
---|
| 1364 | DO i = nxlg, nxrg |
---|
| 1365 | DO j = nysg, nyng |
---|
| 1366 | chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp ) |
---|
| 1367 | ENDDO |
---|
| 1368 | ENDDO |
---|
| 1369 | CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp ) |
---|
| 1370 | ENDIF |
---|
| 1371 | ENDDO |
---|
| 1372 | |
---|
| 1373 | ENDIF |
---|
| 1374 | |
---|
| 1375 | END SUBROUTINE chem_3d_data_averaging |
---|
| 1376 | |
---|
| 1377 | !------------------------------------------------------------------------------! |
---|
| 1378 | ! |
---|
| 1379 | ! Description: |
---|
| 1380 | ! ------------ |
---|
| 1381 | !> Subroutine to write restart data for chemistry model |
---|
| 1382 | !------------------------------------------------------------------------------! |
---|
[2668] | 1383 | SUBROUTINE chem_last_actions |
---|
[2657] | 1384 | |
---|
| 1385 | |
---|
| 1386 | USE control_parameters |
---|
| 1387 | |
---|
| 1388 | USE kinds |
---|
| 1389 | |
---|
| 1390 | IMPLICIT NONE |
---|
| 1391 | |
---|
| 1392 | INTEGER(iwp) :: lsp !< |
---|
[2668] | 1393 | CHARACTER(LEN=20) :: cspcs_name |
---|
| 1394 | CHARACTER(LEN=20) :: cspcs_name_av |
---|
[2657] | 1395 | ! REAL(kind=wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: chems_conc |
---|
| 1396 | |
---|
| 1397 | IF ( write_binary ) THEN |
---|
[2668] | 1398 | DO lsp = 1, nspec |
---|
| 1399 | cspcs_name = ' ' |
---|
| 1400 | cspcs_name_av = ' ' |
---|
| 1401 | cspcs_name = TRIM(chem_species(lsp)%name) |
---|
| 1402 | cspcs_name_av = TRIM(chem_species(lsp)%name)//'_av' |
---|
| 1403 | |
---|
[2657] | 1404 | WRITE(14) cspcs_name ; WRITE(14) chem_species(lsp)%conc |
---|
| 1405 | WRITE(14) cspcs_name_av; WRITE(14) chem_species(lsp)%conc_av |
---|
| 1406 | ENDDO |
---|
| 1407 | |
---|
[2668] | 1408 | WRITE ( 14 ) '*** end chem *** ' |
---|
[2657] | 1409 | |
---|
| 1410 | ENDIF |
---|
| 1411 | |
---|
[2668] | 1412 | END SUBROUTINE chem_last_actions |
---|
[2657] | 1413 | |
---|
| 1414 | !------------------------------------------------------------------------------! |
---|
| 1415 | ! |
---|
| 1416 | ! Description: |
---|
| 1417 | ! ------------ |
---|
| 1418 | !> Subroutine to read restart data of chemical species |
---|
| 1419 | !------------------------------------------------------------------------------! |
---|
| 1420 | |
---|
[2668] | 1421 | SUBROUTINE chem_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file, & |
---|
| 1422 | nynfa, nyn_on_file, nysfa, nys_on_file, & |
---|
| 1423 | offset_xa, offset_ya, overlap_count, tmp_2d, tmp_3d) |
---|
[2657] | 1424 | |
---|
| 1425 | USE control_parameters |
---|
| 1426 | |
---|
| 1427 | USE indices |
---|
| 1428 | |
---|
| 1429 | USE kinds |
---|
| 1430 | |
---|
| 1431 | USE pegrid |
---|
| 1432 | |
---|
| 1433 | IMPLICIT NONE |
---|
| 1434 | |
---|
| 1435 | CHARACTER (LEN=20) :: field_char !< |
---|
| 1436 | CHARACTER (LEN=20) :: spc_name_av !< |
---|
| 1437 | |
---|
| 1438 | INTEGER(iwp) :: i, lsp !< |
---|
| 1439 | INTEGER(iwp) :: k !< |
---|
| 1440 | INTEGER(iwp) :: nxlc !< |
---|
| 1441 | INTEGER(iwp) :: nxlf !< |
---|
| 1442 | INTEGER(iwp) :: nxl_on_file !< |
---|
| 1443 | INTEGER(iwp) :: nxrc !< |
---|
| 1444 | INTEGER(iwp) :: nxrf !< |
---|
| 1445 | INTEGER(iwp) :: nxr_on_file !< |
---|
| 1446 | INTEGER(iwp) :: nync !< |
---|
| 1447 | INTEGER(iwp) :: nynf !< |
---|
| 1448 | INTEGER(iwp) :: nyn_on_file !< |
---|
| 1449 | INTEGER(iwp) :: nysc !< |
---|
| 1450 | INTEGER(iwp) :: nysf !< |
---|
| 1451 | INTEGER(iwp) :: nys_on_file !< |
---|
| 1452 | INTEGER(iwp) :: overlap_count !< |
---|
| 1453 | |
---|
| 1454 | INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxlfa !< |
---|
| 1455 | INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxrfa !< |
---|
| 1456 | INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nynfa !< |
---|
| 1457 | INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nysfa !< |
---|
| 1458 | INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_xa !< |
---|
| 1459 | INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_ya !< |
---|
| 1460 | !! |
---|
| 1461 | REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp, & !< 2D array to temp store data |
---|
| 1462 | nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d !< |
---|
| 1463 | REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp, & !< |
---|
| 1464 | nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !< 3D array to temp store data |
---|
| 1465 | |
---|
| 1466 | |
---|
| 1467 | IF ( initializing_actions == 'read_restart_data' ) THEN |
---|
| 1468 | READ ( 13 ) field_char |
---|
| 1469 | |
---|
| 1470 | DO WHILE ( TRIM( field_char ) /= '*** end chem ***') |
---|
| 1471 | DO k = 1, overlap_count |
---|
| 1472 | nxlf = nxlfa(i,k) |
---|
| 1473 | nxlc = nxlfa(i,k) + offset_xa(i,k) |
---|
| 1474 | nxrf = nxrfa(i,k) |
---|
| 1475 | nxrc = nxrfa(i,k) + offset_xa(i,k) |
---|
| 1476 | nysf = nysfa(i,k) |
---|
| 1477 | nysc = nysfa(i,k) + offset_ya(i,k) |
---|
| 1478 | nynf = nynfa(i,k) |
---|
| 1479 | nync = nynfa(i,k) + offset_ya(i,k) |
---|
| 1480 | |
---|
[2668] | 1481 | DO lsp = 1, nspec |
---|
[2657] | 1482 | spc_name_av = TRIM(chem_species(lsp)%name)//'_av' |
---|
[2668] | 1483 | IF (TRIM( field_char ) == TRIM(chem_species(lsp)%name) ) THEN |
---|
[2657] | 1484 | IF ( k == 1 ) READ ( 13 ) tmp_3d !< read data into tmp_3d |
---|
| 1485 | chem_species(lsp)%conc = tmp_3d !< fill ..%conc in the restart run |
---|
| 1486 | ELSEIF (TRIM( field_char ) == spc_name_av ) THEN |
---|
| 1487 | IF( k == 1 ) READ ( 13 ) tmp_3d |
---|
| 1488 | chem_species(lsp)%conc_av = tmp_3d |
---|
| 1489 | ENDIF |
---|
| 1490 | |
---|
| 1491 | ENDDO |
---|
| 1492 | |
---|
| 1493 | ENDDO |
---|
| 1494 | |
---|
| 1495 | READ ( 13 ) field_char |
---|
| 1496 | |
---|
| 1497 | ENDDO |
---|
| 1498 | |
---|
[2668] | 1499 | ! IF (TRIM( field_char ) == TRIM('surface_csflux') ) THEN |
---|
| 1500 | ! READ ( 13 ) surface_csflux |
---|
| 1501 | ! ENDIF |
---|
[2657] | 1502 | |
---|
| 1503 | ENDIF |
---|
| 1504 | |
---|
[2668] | 1505 | END SUBROUTINE chem_read_restart_data |
---|
[2657] | 1506 | ! |
---|
| 1507 | !------------------------------------------------------------------------------! |
---|
| 1508 | ! |
---|
| 1509 | ! Description: |
---|
| 1510 | ! ------------ |
---|
| 1511 | !> Subroutine calculating tendencies for chemical species |
---|
| 1512 | !------------------------------------------------------------------------------! |
---|
| 1513 | |
---|
[2768] | 1514 | SUBROUTINE chem_tendency ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs, & |
---|
| 1515 | i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs, & |
---|
| 1516 | flux_l_cs, diss_l_cs ) |
---|
[2657] | 1517 | USE pegrid |
---|
| 1518 | USE advec_ws, ONLY: advec_s_ws |
---|
| 1519 | USE advec_s_pw_mod, ONLY: advec_s_pw |
---|
| 1520 | USE advec_s_up_mod, ONLY: advec_s_up |
---|
| 1521 | USE diffusion_s_mod, ONLY: diffusion_s |
---|
| 1522 | USE indices, ONLY: wall_flags_0 |
---|
[2768] | 1523 | USE surface_mod, ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & |
---|
| 1524 | surf_usm_v |
---|
[2657] | 1525 | !> Only one chem_spcs recieved prog_eqns at time. |
---|
| 1526 | !> cem_tendency is called in prog_eqns over loop => nvar |
---|
| 1527 | |
---|
| 1528 | IMPLICIT NONE |
---|
| 1529 | |
---|
[2768] | 1530 | REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar_p, cs_scalar, tcs_scalar_m |
---|
[2657] | 1531 | |
---|
| 1532 | INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn, ilsp |
---|
| 1533 | REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: flux_s_cs !< |
---|
| 1534 | REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: diss_s_cs !< |
---|
| 1535 | REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: flux_l_cs !< |
---|
| 1536 | REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: diss_l_cs !< |
---|
| 1537 | REAL(wp), DIMENSION(0:nz+1) :: pr_init_cs !< |
---|
| 1538 | |
---|
| 1539 | !-- local variables |
---|
| 1540 | |
---|
| 1541 | INTEGER :: k |
---|
| 1542 | ! |
---|
| 1543 | !-- Tendency-terms for chem spcs. |
---|
| 1544 | tend(:,j,i) = 0.0_wp |
---|
| 1545 | ! |
---|
| 1546 | !-- Advection terms |
---|
| 1547 | IF ( timestep_scheme(1:5) == 'runge' ) THEN |
---|
| 1548 | IF ( ws_scheme_sca ) THEN |
---|
| 1549 | CALL advec_s_ws( i, j, cs_scalar, 'kc', flux_s_cs, diss_s_cs, & |
---|
| 1550 | flux_l_cs, diss_l_cs, i_omp_start, tn ) |
---|
| 1551 | ELSE |
---|
| 1552 | CALL advec_s_pw( i, j, cs_scalar ) |
---|
| 1553 | ENDIF |
---|
| 1554 | ELSE |
---|
| 1555 | CALL advec_s_up( i, j, cs_scalar ) |
---|
| 1556 | ENDIF |
---|
| 1557 | |
---|
| 1558 | ! |
---|
| 1559 | |
---|
| 1560 | !-- Diffusion terms ( ie hinteren 3 sind 0 ) |
---|
| 1561 | |
---|
| 1562 | CALL diffusion_s( i, j, cs_scalar, & |
---|
| 1563 | surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:), & |
---|
| 1564 | surf_def_h(2)%cssws(ilsp,:), & |
---|
| 1565 | surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:), & |
---|
| 1566 | surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:), & |
---|
| 1567 | surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:), & |
---|
| 1568 | surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:), & |
---|
| 1569 | surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:), & |
---|
| 1570 | surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:), & |
---|
| 1571 | surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) ) |
---|
| 1572 | |
---|
| 1573 | ! |
---|
| 1574 | !-- Prognostic equation for chem spcs |
---|
| 1575 | DO k = nzb+1, nzt |
---|
| 1576 | cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d * & |
---|
| 1577 | ( tsc(2) * tend(k,j,i) + & |
---|
| 1578 | tsc(3) * tcs_scalar_m(k,j,i) ) & |
---|
| 1579 | - tsc(5) * rdf_sc(k) & |
---|
| 1580 | * ( cs_scalar(k,j,i) - pr_init_cs(k) ) & |
---|
| 1581 | ) & |
---|
| 1582 | * MERGE( 1.0_wp, 0.0_wp, & |
---|
| 1583 | BTEST( wall_flags_0(k,j,i), 0 ) & |
---|
| 1584 | ) |
---|
| 1585 | |
---|
| 1586 | 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 |
---|
| 1587 | ENDDO |
---|
| 1588 | |
---|
| 1589 | ! |
---|
| 1590 | !-- Calculate tendencies for the next Runge-Kutta step |
---|
| 1591 | IF ( timestep_scheme(1:5) == 'runge' ) THEN |
---|
| 1592 | IF ( intermediate_timestep_count == 1 ) THEN |
---|
| 1593 | DO k = nzb+1, nzt |
---|
| 1594 | tcs_scalar_m(k,j,i) = tend(k,j,i) |
---|
| 1595 | ENDDO |
---|
| 1596 | ELSEIF ( intermediate_timestep_count < & |
---|
| 1597 | intermediate_timestep_count_max ) THEN |
---|
| 1598 | DO k = nzb+1, nzt |
---|
| 1599 | tcs_scalar_m(k,j,i) = -9.5625_wp * tend(k,j,i) + & |
---|
| 1600 | 5.3125_wp * tcs_scalar_m(k,j,i) |
---|
| 1601 | ENDDO |
---|
| 1602 | ENDIF |
---|
| 1603 | ENDIF |
---|
| 1604 | ! |
---|
| 1605 | END SUBROUTINE chem_tendency |
---|
| 1606 | |
---|
| 1607 | !------------------------------------------------------------------------------! |
---|
| 1608 | ! |
---|
| 1609 | ! Description: |
---|
| 1610 | ! ------------ |
---|
| 1611 | !> Subroutine defining header output for chemistry model |
---|
| 1612 | !------------------------------------------------------------------------------! |
---|
| 1613 | SUBROUTINE chem_header ( io ) |
---|
| 1614 | |
---|
| 1615 | IMPLICIT NONE |
---|
| 1616 | |
---|
| 1617 | INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file |
---|
| 1618 | |
---|
| 1619 | ! print*,'the header subroutine is still not operational' |
---|
| 1620 | !! |
---|
| 1621 | !!-- Write chemistry model header |
---|
| 1622 | ! WRITE( io, 3 ) |
---|
| 1623 | ! |
---|
| 1624 | ! IF ( radiation_scheme == "constant" ) THEN |
---|
| 1625 | ! WRITE( io, 4 ) net_radiation |
---|
| 1626 | ! ELSEIF ( radiation_scheme == "clear-sky" ) THEN |
---|
| 1627 | ! WRITE( io, 5 ) |
---|
| 1628 | ! ELSEIF ( radiation_scheme == "rrtmg" ) THEN |
---|
| 1629 | ! WRITE( io, 6 ) |
---|
| 1630 | ! IF ( .NOT. lw_radiation ) WRITE( io, 10 ) |
---|
| 1631 | ! IF ( .NOT. sw_radiation ) WRITE( io, 11 ) |
---|
| 1632 | ! ENDIF |
---|
| 1633 | ! |
---|
| 1634 | ! IF ( albedo_type == 0 ) THEN |
---|
| 1635 | ! WRITE( io, 7 ) albedo |
---|
| 1636 | ! ELSE |
---|
| 1637 | ! WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) ) |
---|
| 1638 | ! ENDIF |
---|
| 1639 | ! IF ( constant_albedo ) THEN |
---|
| 1640 | ! WRITE( io, 9 ) |
---|
| 1641 | ! ENDIF |
---|
| 1642 | ! |
---|
| 1643 | ! IF ( radiation .AND. radiation_scheme /= 'constant' ) THEN |
---|
| 1644 | ! WRITE ( io, 1 ) lambda |
---|
| 1645 | ! WRITE ( io, 2 ) day_init, time_utc_init |
---|
| 1646 | ! ENDIF |
---|
| 1647 | ! |
---|
| 1648 | ! WRITE( io, 12 ) dt_radiation |
---|
| 1649 | ! |
---|
| 1650 | ! 1 FORMAT (' Geograph. longitude : lambda = ',F4.1,' degr') |
---|
| 1651 | ! 2 FORMAT (' Day of the year at model start : day_init = ',I3 & |
---|
| 1652 | ! /' UTC time at model start : time_utc_init = ',F7.1' s') |
---|
| 1653 | ! 3 FORMAT (//' Radiation model information:'/ & |
---|
| 1654 | ! ' ----------------------------'/) |
---|
| 1655 | ! 4 FORMAT (' --> Using constant net radiation: net_radiation = ', F6.2, & |
---|
| 1656 | ! // 'W/m**2') |
---|
| 1657 | ! 5 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,', & |
---|
| 1658 | ! ' default)') |
---|
| 1659 | ! 6 FORMAT (' --> RRTMG scheme is used') |
---|
| 1660 | ! 7 FORMAT (/' User-specific surface albedo: albedo =', F6.3) |
---|
| 1661 | ! 8 FORMAT (/' Albedo is set for land surface type: ', A) |
---|
| 1662 | ! 9 FORMAT (/' --> Albedo is fixed during the run') |
---|
| 1663 | !10 FORMAT (/' --> Longwave radiation is disabled') |
---|
| 1664 | !11 FORMAT (/' --> Shortwave radiation is disabled.') |
---|
| 1665 | !12 FORMAT (' Timestep: dt_radiation = ', F6.2, ' s') |
---|
| 1666 | ! |
---|
| 1667 | ! |
---|
| 1668 | END SUBROUTINE chem_header |
---|
[2668] | 1669 | |
---|
| 1670 | !------------------------------------------------------------------------------ |
---|
| 1671 | ! Description: |
---|
| 1672 | ! ------------ |
---|
| 1673 | !> Subroutine reading restart data for chemistry model input parameters |
---|
| 1674 | ! (FK: To make restarts work, I had to comment this routine. We actually |
---|
| 1675 | ! don't need it, since the namelist parameters are always read in, |
---|
| 1676 | ! also in case of a restart run) |
---|
| 1677 | !------------------------------------------------------------------------------ |
---|
| 1678 | ! SUBROUTINE chem_read_var_list |
---|
| 1679 | ! |
---|
| 1680 | ! ! USE kinds |
---|
| 1681 | ! ! USE chem_modules |
---|
| 1682 | ! |
---|
| 1683 | ! |
---|
| 1684 | ! IMPLICIT NONE |
---|
| 1685 | ! |
---|
| 1686 | ! CHARACTER (LEN=30) :: variable_chr !< dummy variable to read string |
---|
| 1687 | ! |
---|
| 1688 | ! |
---|
| 1689 | ! READ ( 13 ) variable_chr |
---|
| 1690 | ! DO WHILE ( TRIM( variable_chr ) /= '*** end chemistry ***' ) |
---|
| 1691 | ! |
---|
| 1692 | ! SELECT CASE ( TRIM( variable_chr ) ) |
---|
| 1693 | ! |
---|
| 1694 | ! CASE ( 'bc_cs_b' ) |
---|
| 1695 | ! READ ( 13 ) bc_cs_b |
---|
| 1696 | ! |
---|
| 1697 | ! |
---|
| 1698 | ! END SELECT |
---|
| 1699 | ! |
---|
| 1700 | ! READ ( 13 ) variable_chr |
---|
| 1701 | ! |
---|
| 1702 | ! ENDDO |
---|
| 1703 | ! |
---|
| 1704 | ! END SUBROUTINE chem_read_var_list |
---|
| 1705 | |
---|
| 1706 | !------------------------------------------------------------------------------ |
---|
| 1707 | ! Description: |
---|
| 1708 | ! ------------ |
---|
| 1709 | !> Skipping the chemistry-module parameters from restart-file (binary format). |
---|
| 1710 | ! (FK: To make restarts work, I had to comment this routine. We actually |
---|
| 1711 | ! don't need it, since the namelist parameters are always read in, |
---|
| 1712 | ! also in case of a restart run) |
---|
| 1713 | !------------------------------------------------------------------------------ |
---|
| 1714 | ! SUBROUTINE chem_skip_var_list |
---|
| 1715 | ! |
---|
| 1716 | ! IMPLICIT NONE |
---|
| 1717 | ! |
---|
| 1718 | ! CHARACTER (LEN=1) :: cdum |
---|
| 1719 | ! CHARACTER (LEN=30) :: variable_chr |
---|
| 1720 | ! |
---|
| 1721 | ! READ ( 13 ) variable_chr |
---|
| 1722 | ! |
---|
| 1723 | ! DO WHILE ( TRIM( variable_chr ) /= '*** end chemistry ***' ) |
---|
| 1724 | ! |
---|
| 1725 | ! READ ( 13 ) cdum |
---|
| 1726 | ! READ ( 13 ) variable_chr |
---|
| 1727 | ! |
---|
| 1728 | ! ENDDO |
---|
| 1729 | ! |
---|
| 1730 | ! END SUBROUTINE chem_skip_var_list |
---|
| 1731 | |
---|
| 1732 | |
---|
| 1733 | !------------------------------------------------------------------------------! |
---|
[2657] | 1734 | ! |
---|
[2668] | 1735 | ! Description: |
---|
| 1736 | ! ------------ |
---|
| 1737 | !> Subroutine writing restart data for chemistry model input parameters |
---|
| 1738 | ! (FK: To make restarts work, I had to comment this routine. We actually |
---|
| 1739 | ! don't need it, since the namelist parameters are always read in, |
---|
| 1740 | ! also in case of a restart run) |
---|
[2657] | 1741 | !------------------------------------------------------------------------------! |
---|
[2668] | 1742 | ! SUBROUTINE chem_write_var_list |
---|
| 1743 | ! |
---|
| 1744 | ! USE chem_modules |
---|
| 1745 | ! |
---|
| 1746 | ! USE kinds |
---|
| 1747 | ! |
---|
| 1748 | ! |
---|
| 1749 | ! IMPLICIT NONE |
---|
| 1750 | ! |
---|
| 1751 | ! INTEGER(iwp) :: lsp !< running index for chem spcs |
---|
| 1752 | ! |
---|
| 1753 | ! ! |
---|
| 1754 | ! !-- Writing out input parameters that are not part of chemistry_par namelist |
---|
| 1755 | ! !-- (namelist parameters are anyway read in again in case of restart) |
---|
| 1756 | ! DO lsp = 1, nvar |
---|
| 1757 | ! WRITE ( 14 ) 'conc_pr_init_'//chem_species(lsp)%name |
---|
| 1758 | ! WRITE ( 14 ) chem_species(lsp)%conc_pr_init |
---|
| 1759 | ! ENDDO |
---|
| 1760 | ! |
---|
| 1761 | ! WRITE ( 14 ) '*** end chemistry *** ' |
---|
| 1762 | ! |
---|
| 1763 | ! END SUBROUTINE chem_write_var_list |
---|
| 1764 | |
---|
| 1765 | |
---|
| 1766 | !------------------------------------------------------------------------------! |
---|
[2657] | 1767 | ! |
---|
| 1768 | ! Description: |
---|
| 1769 | ! ------------ |
---|
[2668] | 1770 | !> Subroutine for emission |
---|
[2657] | 1771 | !------------------------------------------------------------------------------! |
---|
[2668] | 1772 | SUBROUTINE chem_emissions |
---|
[2657] | 1773 | |
---|
[2668] | 1774 | USE chem_modules |
---|
| 1775 | |
---|
| 1776 | USE netcdf_data_input_mod, & |
---|
| 1777 | ONLY: street_type_f |
---|
| 1778 | |
---|
| 1779 | USE surface_mod, & |
---|
| 1780 | ONLY: surf_lsm_h |
---|
[2657] | 1781 | |
---|
| 1782 | |
---|
[2668] | 1783 | IMPLICIT NONE |
---|
[2657] | 1784 | |
---|
[2678] | 1785 | INTEGER(iwp) :: i !< running index for grid in x-direction |
---|
| 1786 | INTEGER(iwp) :: j !< running index for grid in y-direction |
---|
| 1787 | INTEGER(iwp) :: m !< running index for horizontal surfaces |
---|
[2668] | 1788 | INTEGER(iwp) :: lsp !< running index for chem spcs |
---|
[2657] | 1789 | |
---|
[2668] | 1790 | ! |
---|
| 1791 | !-- Comment??? (todo) |
---|
[2678] | 1792 | IF ( street_type_f%from_file ) THEN |
---|
| 1793 | ! |
---|
| 1794 | !-- Streets are lsm surfaces, hence, no usm surface treatment required |
---|
| 1795 | DO m = 1, surf_lsm_h%ns |
---|
| 1796 | i = surf_lsm_h%i(m) |
---|
| 1797 | j = surf_lsm_h%j(m) |
---|
| 1798 | |
---|
| 1799 | IF ( street_type_f%var(j,i) >= main_street_id .AND. & |
---|
| 1800 | street_type_f%var(j,i) < max_street_id ) THEN |
---|
| 1801 | DO lsp = 1, nvar |
---|
| 1802 | surf_lsm_h%cssws(lsp,m) = emiss_factor_main * surface_csflux(lsp) |
---|
| 1803 | ENDDO |
---|
| 1804 | ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. & |
---|
| 1805 | street_type_f%var(j,i) < main_street_id ) THEN |
---|
| 1806 | DO lsp = 1, nvar |
---|
| 1807 | surf_lsm_h%cssws(lsp,m) = emiss_factor_side * surface_csflux(lsp) |
---|
| 1808 | ENDDO |
---|
[2682] | 1809 | ELSE |
---|
| 1810 | surf_lsm_h%cssws(:,m) = 0.0_wp |
---|
[2678] | 1811 | ENDIF |
---|
| 1812 | ENDDO |
---|
| 1813 | |
---|
| 1814 | ENDIF |
---|
[2657] | 1815 | |
---|
[2668] | 1816 | END SUBROUTINE chem_emissions |
---|
[2657] | 1817 | |
---|
| 1818 | |
---|
| 1819 | END MODULE chemistry_model_mod |
---|
| 1820 | |
---|