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

Last change on this file since 2894 was 2894, checked in by Giersch, 4 years ago

Reading/Writing? data in case of restart runs revised

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