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

Last change on this file since 2756 was 2756, checked in by suehring, 4 years ago

Fill values for 3D data output of chemical species introduced

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