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

Last change on this file since 2932 was 2932, checked in by maronga, 3 years ago

renamed all Fortran NAMELISTS

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