Index: /palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 2814)
+++ /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 2815)
@@ -27,4 +27,10 @@
 ! -----------------
 ! $Id$
+! Bugfix in restart mechanism,
+! rename chem_tendency to chem_prognostic_equations,
+! implement vector-optimized version of chem_prognostic_equations,
+! some clean up (incl. todo list)
+! 
+! 2773 2018-01-30 14:12:54Z suehring
 ! Declare variables required for nesting as public
 ! 
@@ -60,9 +66,8 @@
 !> Chemistry model for PALM-4U
 !> @todo Update/clean-up todo list! (FK)
-!> @todo Chemistry for prognostic_equations_vector! (FK)
 !> @todo Set proper fill values (/= 0) for chem output arrays! (FK)
 !> @todo Add routine chem_check_parameters, add checks for inconsistent or
-!>       unallowed parameter settings (e.g. "only cache optimized version
-!>       available"). CALL of chem_check_parameters from check_parameters. (FK)
+!>       unallowed parameter settings. 
+!>       CALL of chem_check_parameters from check_parameters. (FK)
 !> @todo Make routine chem_header available, CALL from header.f90
 !>       (see e.g. how it is done in routine lsm_header in 
@@ -76,5 +81,5 @@
 !> @todo slight differences in passive scalar and chem spcs when chem reactions
 !>       turned off. Need to be fixed. bK
-!> @todo introduce nesting for chem spcs. bK
+!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
 !> @todo subroutine set_const_initial_values to be taken out from chemistry_model_mod !bK.
 !> @todo chemistry error messages
@@ -211,7 +216,8 @@
    END INTERFACE chem_read_restart_data
 
-   INTERFACE chem_tendency
-      MODULE PROCEDURE chem_tendency
-   END INTERFACE chem_tendency
+   INTERFACE chem_prognostic_equations
+      MODULE PROCEDURE chem_prognostic_equations
+      MODULE PROCEDURE chem_prognostic_equations_ij
+   END INTERFACE chem_prognostic_equations
 
    INTERFACE chem_header
@@ -236,11 +242,11 @@
 
 
-   PUBLIC chem_boundary_conds, chem_check_data_output,                         &
+   PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, &
           chem_check_data_output_pr, chem_data_output_3d,                      &
-          chem_define_netcdf_grid, chem_header, chem_init,                     &
-          chem_init_profiles, chem_integrate, chem_parin,                      &
-          chem_swap_timelevel, chem_last_actions, chem_read_restart_data,      &
-          chem_tendency, chem_3d_data_averaging, chem_emissions!, chem_write_var_list,          &
-          !chem_read_var_list, chem_skip_var_list
+          chem_define_netcdf_grid, chem_emissions, chem_header, chem_init,     &
+          chem_init_profiles, chem_integrate, chem_last_actions,               &
+          chem_parin, chem_prognostic_equations,                               &
+          chem_read_restart_data, chem_swap_timelevel
+          
 
 
@@ -293,5 +299,5 @@
           ELSE
 !             message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"'  ! bK commented
-             CALL message( 'chem_boundary_conds', 'CHEM001', 1, 2, 0, 6, 0 )     !< chemistry_model_mod should have special error numbers --> "CHEM###",
+             CALL message( 'chem_boundary_conds', 'CM0010', 1, 2, 0, 6, 0 )     !< chemistry_model_mod should have special error numbers --> "CHEM###",
           ENDIF                                                                  
 !
@@ -308,5 +314,5 @@
           ELSE
 !            message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'  
-             CALL message( 'check_parameters', 'CHEM002', 1, 2, 0, 6, 0 )
+             CALL message( 'check_parameters', 'CM0011', 1, 2, 0, 6, 0 )
           ENDIF
 
@@ -449,5 +455,5 @@
 !                   IF ( cs_heights(lsp,1) /= 0.0_wp )  THEN
 !                      message_string = 'cs_heights(1,1) must be 0.0'
-!                      CALL message( 'chem_check_parameters', 'CHEM0004', 1, 2, 0, 6, 0 )
+!                      CALL message( 'chem_check_parameters', 'CM0012', 1, 2, 0, 6, 0 )
 !                   ENDIF
 !  
@@ -658,90 +664,5 @@
           chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
        ENDDO
-!        CALL location_message( 'finished', .TRUE. )
-
-!
-!-- (todo (FK): Restarts not available yet. It is not correct to call 
-!--             read_3d_binary here)
-    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
-             TRIM( initializing_actions ) == 'cyclic_fill' )                   &
-    THEN
-
-       message_string = 'so far, restarts are not possible with chemistry'
-       CALL message( 'chem_init', 'CM0002', 1, 2, 0, 6, 0 )
-!        CALL location_message( 'initializing in case of restart / cyclic_fill', &
-!                               .FALSE. )
-! !
-! !--    When reading data for cyclic fill of 3D prerun data files, read
-! !--    some of the global variables from the restart file which are required
-! !--    for initializing the inflow
-!        IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
-!           DO  i = 0, io_blocks-1
-!              IF ( i == io_group )  THEN
-!                 CALL read_parts_of_var_list
-!                 CALL close_file( 13 )
-!              ENDIF
-! #if defined( __parallel )
-!              CALL MPI_BARRIER( comm2d, ierr )
-! #endif
-!           ENDDO
-!        ENDIF
-! !
-! !--    Read binary data from restart file
-!        DO  i = 0, io_blocks-1
-!           IF ( i == io_group )  THEN
-!              CALL read_3d_binary
-!           ENDIF
-! #if defined( __parallel )
-!           CALL MPI_BARRIER( comm2d, ierr )
-! #endif
-!        ENDDO
-! !
-! !--    Initialization of the turbulence recycling method
-!        IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.               &
-!             turbulent_inflow )  THEN
-! !
-! !--       First store the profiles to be used at the inflow.
-! !--       These profiles are the (temporally) and horizontally averaged vertical
-! !--       profiles from the prerun. Alternatively, prescribed profiles
-! !--       for u,v-components can be used.
-!           ALLOCATE( mean_inflow_profiles(nzb:nzt+1,8) )
-! !         mean_inflow_profiles(:,8) = hom_sum(:,115,0)   !cs           
-! 
-! !         IF ( inflow_l) THEN
-! !            DO  j = nysg, nyng
-! !               DO  k = nzb, nzt+1
-! !                     DO lsp = 1, nvar 
-! !                       chem_species(lsp)%conc(k,j,nxlg:-1)  = mean_inflow_profiles(k,8)  ???
-! !                     ENDDO
-! !               ENDDO
-! !            ENDDO
-! !         ENDIF
-!        ENDIF
-! 
-! !--    Inside buildings set velocities and TKE back to zero         
-! !       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.                &
-! !            topography /= 'flat' )  THEN
-! !           ......
-! !       ENDIF
-! ! 
-! !--    Initialize new time levels (only done in order to set boundary values
-! !--    including ghost points)
-!        DO lsp = 1, nvar
-!           chem_species(lsp)%conc_p  = chem_species(lsp)%conc
-!        ENDDO
-! !
-! !--    Allthough tendency arrays are set in prognostic_equations, they have
-! !--    have to be predefined here because they are used (but multiplied with 0)
-! !--    there before they are set. 
-! 
-!        DO lsp = 1, nvar
-!           chem_species(lsp)%tconc_m  = 0.0_wp
-!        ENDDO
-
-    ELSE
-!
-!--    Actually this part of the programm should not be reached
-       message_string = 'unknown initializing problem'
-       CALL message( 'chem_init', 'CM0003', 1, 2, 0, 6, 0 )
+
     ENDIF
 
@@ -1397,9 +1318,10 @@
     IMPLICIT NONE 
 
-    INTEGER(iwp) :: lsp !< 
+    INTEGER(iwp) ::  lsp !< 
     CHARACTER(LEN=20) ::  cspcs_name
     CHARACTER(LEN=20) ::  cspcs_name_av
 !      REAL(kind=wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   :: chems_conc
-  
+
+
     IF ( write_binary )  THEN
        DO  lsp = 1, nspec
@@ -1428,5 +1350,6 @@
  SUBROUTINE chem_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file, &
                                     nynfa, nyn_on_file, nysfa, nys_on_file,    &
-                                    offset_xa, offset_ya, overlap_count, tmp_2d, tmp_3d)    
+                                    offset_xa, offset_ya, overlap_count,       &
+                                    tmp_2d, tmp_3d )    
                                      
     USE control_parameters
@@ -1464,17 +1387,23 @@
     INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nysfa       !< 
     INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_xa   !<   
-    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !<   
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !<  
+    
+    LOGICAL ::  chem_found 
 !!
-    REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,                    &   !< 2D array to temp store data 
-                        nxl_on_file-nbgp:nxr_on_file+nbgp)     ::    tmp_2d       !<
-    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
+    REAL(wp),                                                                  &
+       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
+          tmp_2d   !< 2D array to temp store data
+
+    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
 
 
     IF ( initializing_actions == 'read_restart_data' )  THEN
        READ ( 13 )  field_char
-
-       DO WHILE ( TRIM( field_char ) /= '*** end chem ***')
-          DO k = 1, overlap_count
+       DO  WHILE ( TRIM( field_char ) /= '*** end chem ***    ' )
+
+          DO  k = 1, overlap_count
+
              nxlf = nxlfa(i,k)
              nxlc = nxlfa(i,k) + offset_xa(i,k)
@@ -1485,17 +1414,30 @@
              nynf = nynfa(i,k)
              nync = nynfa(i,k) + offset_ya(i,k)
-
+             
+
+             chem_found = .FALSE.
+  
              DO  lsp = 1, nspec
-                spc_name_av  =  TRIM(chem_species(lsp)%name)//'_av'
+
+                spc_name_av  =  TRIM(chem_species(lsp)%name)//'_av'   !< for time-averaged chemical conc.
                 IF (TRIM( field_char ) == TRIM(chem_species(lsp)%name) )  THEN
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d                             !< read data into tmp_3d
-                   chem_species(lsp)%conc = tmp_3d                                !< fill ..%conc in the restart run
-                ELSEIF (TRIM( field_char ) == spc_name_av ) THEN
-                   IF( k == 1 ) READ ( 13 ) tmp_3d
-                   chem_species(lsp)%conc_av = tmp_3d
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =           &
+                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                   chem_found = .TRUE.
+                ELSEIF (TRIM( field_char ) == spc_name_av )  THEN
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
+                                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                   chem_found = .TRUE.
                 ENDIF
 
              ENDDO
-
+             IF (  .NOT.  chem_found )  THEN
+                WRITE( message_string, * ) 'unknown variable named "',         &
+                                        TRIM( field_char ), '" found in',      &
+                                        '&data from prior run on PE ', myid
+                CALL message( 'chem_read_restart_data', 'CM0008', 1, 2, 0, 6, 0 )
+             ENDIF
           ENDDO
 
@@ -1511,13 +1453,16 @@
 
  END SUBROUTINE chem_read_restart_data
-!
-!------------------------------------------------------------------------------!
-!
-! Description:
-! ------------
-!> Subroutine calculating tendencies for chemical species
-!------------------------------------------------------------------------------!
-
- SUBROUTINE chem_tendency ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs,  &
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine calculating prognostic equations for chemical species 
+!> (cache-optimized).
+!> Routine is called separately for each chemical species over a loop from
+!> prognostic_equations.
+!------------------------------------------------------------------------------!
+ SUBROUTINE chem_prognostic_equations_ij ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs,  &
                             i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs, &
                             flux_l_cs, diss_l_cs )
@@ -1530,6 +1475,5 @@
     USE surface_mod,     ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
                                 surf_usm_v
-!> Only one chem_spcs recieved prog_eqns at time.
-!> cem_tendency is called in prog_eqns over loop => nvar
+
 
     IMPLICIT NONE
@@ -1565,5 +1509,5 @@
 !
 
-!-- Diffusion terms ( ie hinteren 3 sind 0 )
+!-- Diffusion terms (the last three arguments are zero)
 
       CALL diffusion_s( i, j, cs_scalar,                                                 &
@@ -1609,6 +1553,135 @@
        ENDIF
     ENDIF
-!
- END SUBROUTINE chem_tendency
+
+ END SUBROUTINE chem_prognostic_equations_ij
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine calculating prognostic equations for chemical species 
+!> (vector-optimized).
+!> Routine is called separately for each chemical species over a loop from
+!> prognostic_equations.
+!------------------------------------------------------------------------------!
+ SUBROUTINE chem_prognostic_equations ( cs_scalar_p, cs_scalar, tcs_scalar_m,  &
+                                        pr_init_cs, ilsp )
+
+    USE advec_s_pw_mod,                                                        &
+        ONLY:  advec_s_pw
+
+    USE advec_s_up_mod,                                                        &
+        ONLY:  advec_s_up
+
+    USE advec_ws,                                                              &
+        ONLY:  advec_s_ws 
+
+    USE diffusion_s_mod,                                                       &
+        ONLY:  diffusion_s
+
+    USE indices,                                                               &
+        ONLY:  nxl, nxr, nyn, nys, wall_flags_0
+
+    USE pegrid
+
+    USE surface_mod,                                                           &
+        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
+               surf_usm_v
+
+    IMPLICIT NONE
+
+    INTEGER ::  i   !< running index
+    INTEGER ::  j   !< running index
+    INTEGER ::  k   !< running index
+
+    INTEGER(iwp),INTENT(IN) ::  ilsp          !<
+
+    REAL(wp), DIMENSION(0:nz+1) ::  pr_init_cs   !<
+
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  cs_scalar      !<
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  cs_scalar_p    !<
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  tcs_scalar_m   !<
+
+
+!
+!-- Tendency terms for chemical species
+    tend = 0.0_wp
+!    
+!-- Advection terms
+    IF ( timestep_scheme(1:5) == 'runge' )  THEN
+       IF ( ws_scheme_sca )  THEN
+          CALL advec_s_ws( cs_scalar, 'kc' )
+       ELSE
+          CALL advec_s_pw( cs_scalar )
+       ENDIF
+    ELSE
+         CALL advec_s_up( cs_scalar )
+    ENDIF
+!
+!-- Diffusion terms  (the last three arguments are zero)
+    CALL diffusion_s( cs_scalar,                                               &
+                      surf_def_h(0)%cssws(ilsp,:),                             &
+                      surf_def_h(1)%cssws(ilsp,:),                             &
+                      surf_def_h(2)%cssws(ilsp,:),                             &
+                      surf_lsm_h%cssws(ilsp,:),                                &
+                      surf_usm_h%cssws(ilsp,:),                                &
+                      surf_def_v(0)%cssws(ilsp,:),                             &
+                      surf_def_v(1)%cssws(ilsp,:),                             &
+                      surf_def_v(2)%cssws(ilsp,:),                             &
+                      surf_def_v(3)%cssws(ilsp,:),                             &
+                      surf_lsm_v(0)%cssws(ilsp,:),                             &
+                      surf_lsm_v(1)%cssws(ilsp,:),                             &
+                      surf_lsm_v(2)%cssws(ilsp,:),                             &
+                      surf_lsm_v(3)%cssws(ilsp,:),                             &
+                      surf_usm_v(0)%cssws(ilsp,:),                             &
+                      surf_usm_v(1)%cssws(ilsp,:),                             &
+                      surf_usm_v(2)%cssws(ilsp,:),                             &
+                      surf_usm_v(3)%cssws(ilsp,:) )
+!    
+!-- Prognostic equation for chemical species
+    DO  i = nxl, nxr
+       DO  j = nys, nyn    
+          DO  k = nzb+1, nzt
+             cs_scalar_p(k,j,i) =   cs_scalar(k,j,i)                           &
+                                  + ( dt_3d  *                                 &
+                                      (   tsc(2) * tend(k,j,i)                 &
+                                        + tsc(3) * tcs_scalar_m(k,j,i)         &
+                                      )                                        & 
+                                    - tsc(5) * rdf_sc(k)                       &
+                                             * ( cs_scalar(k,j,i) - pr_init_cs(k) )    &   
+                                    )                                          &
+                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )       
+
+             IF ( cs_scalar_p(k,j,i) < 0.0_wp )  cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)
+          ENDDO
+       ENDDO
+    ENDDO
+!
+!-- Calculate tendencies for the next Runge-Kutta step
+    IF ( timestep_scheme(1:5) == 'runge' )  THEN
+       IF ( intermediate_timestep_count == 1 )  THEN
+          DO  i = nxl, nxr
+             DO  j = nys, nyn   
+                DO  k = nzb+1, nzt
+                   tcs_scalar_m(k,j,i) = tend(k,j,i)
+                ENDDO
+             ENDDO
+          ENDDO
+       ELSEIF ( intermediate_timestep_count < &
+          intermediate_timestep_count_max )  THEN
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   tcs_scalar_m(k,j,i) = - 9.5625_wp * tend(k,j,i)             &
+                                         + 5.3125_wp * tcs_scalar_m(k,j,i)
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+    ENDIF
+
+ END SUBROUTINE chem_prognostic_equations
+
 
 !------------------------------------------------------------------------------!
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 2814)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 2815)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Rename chem_tendency to chem_prognostic_equations,
+! implement vector version for air chemistry
+! 
+! 2766 2018-01-22 17:17:47Z kanani
 ! Removed preprocessor directive __chem
 ! 
@@ -270,6 +274,6 @@
 
     USE chemistry_model_mod,                                                   &
-        ONLY:  chem_integrate, chem_species,                                   &
-               chem_tendency, nspec, nvar, spc_names
+        ONLY:  chem_integrate, chem_prognostic_equations,                      &
+               chem_species, nspec, nvar, spc_names
            
     USE chem_photolysis_mod,                                                   &
@@ -1301,6 +1305,6 @@
 !
 !--          Loop over chemical species
-             DO  lsp = 1,nvar                         
-                CALL chem_tendency ( chem_species(lsp)%conc_p,                 &
+             DO  lsp = 1, nvar                         
+                CALL chem_prognostic_equations ( chem_species(lsp)%conc_p,     &
                                      chem_species(lsp)%conc,                   &
                                      chem_species(lsp)%tconc_m,                &
@@ -1340,7 +1344,8 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i    !<
-    INTEGER(iwp) ::  j    !<
-    INTEGER(iwp) ::  k    !<
+    INTEGER(iwp) ::  i     !<
+    INTEGER(iwp) ::  j     !<
+    INTEGER(iwp) ::  k     !<
+    INTEGER(iwp) ::  lsp   !< running index for chemical species
 
     REAL(wp)     ::  sbt  !<
@@ -2436,4 +2441,21 @@
     CALL tcm_prognostic()
 
+!
+!-- If required, compute prognostic equation for chemical quantites
+    IF ( air_chemistry )  THEN
+       CALL cpu_log( log_point(83), '(chem advec+diff+prog)', 'start' )
+!
+!--    Loop over chemical species
+       DO  lsp = 1, nvar                         
+          CALL chem_prognostic_equations ( chem_species(lsp)%conc_p,           &
+                                           chem_species(lsp)%conc,             &
+                                           chem_species(lsp)%tconc_m,          &
+                                           chem_species(lsp)%conc_pr_init,     &
+                                           lsp )       
+       ENDDO
+
+       CALL cpu_log( log_point(83), '(chem advec+diff+prog)', 'stop' )             
+    ENDIF   ! Chemicals equations
+
 
  END SUBROUTINE prognostic_equations_vector
