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

Last change on this file since 2823 was 2815, checked in by kanani, 6 years ago

Enable restarts and vector optimization for chemistry model

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