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

Last change on this file since 3285 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

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