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

Last change on this file since 2805 was 2773, checked in by suehring, 7 years ago

Nesting for chemical species implemented; Bugfix passive scalar boundary condition after anterpolation; Timeseries output of surface temperature; Enable initialization of 3D topography (was commented out so far)

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