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

Last change on this file since 3004 was 3004, checked in by Giersch, 3 years ago

precipitation_rate removed, further allocation checks for data output of averaged quantities implemented, double CALL of flow_statistics at the beginning of time_integration removed, further minor bugfixes, comments added

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