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

Last change on this file since 3209 was 3209, checked in by suehring, 3 years ago

Additional namelist parameter to switch on/off the nesting of chemical species; Enable the input of soil data from dynamic input files independent on atmosphere in order to initialize soil properties in nested child domains from dynamic input; Revise error message number for static/dynamic input; Revise and add checks for static/dynamic input; Bugfix, add netcdf4_parallel directive for collective read operation

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