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

Last change on this file since 3129 was 3045, checked in by Giersch, 6 years ago

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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