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

Last change on this file since 3198 was 3183, checked in by suehring, 6 years ago

last commit documented

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