source: palm/trunk/SOURCE/chemistry_model_mod.f90 @ 3246

Last change on this file since 3246 was 3246, checked in by sward, 3 years ago

Added error handling for wrong input parameters

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